To pytanie zostało zadane przed nieco ponad trzy lata temu. Dano odpowiedź, ale znalazłem usterkę w rozwiązaniu.
Poniższy kod jest w R. Przenieśliłem go na inny język, jednak przetestowałem oryginalny kod bezpośrednio w R, aby upewnić się, że problem nie dotyczył mojego przenoszenia.
sunPosition <- function(year, month, day, hour=12, min=0, sec=0,
lat=46.5, long=6.5) {
twopi <- 2 * pi
deg2rad <- pi / 180
# Get day of the year, e.g. Feb 1 = 32, Mar 1 = 61 on leap years
month.days <- c(0,31,28,31,30,31,30,31,31,30,31,30)
day <- day + cumsum(month.days)[month]
leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & day >= 60
day[leapdays] <- day[leapdays] + 1
# Get Julian date - 2400000
hour <- hour + min / 60 + sec / 3600 # hour plus fraction
delta <- year - 1949
leap <- trunc(delta / 4) # former leapyears
jd <- 32916.5 + delta * 365 + leap + day + hour / 24
# The input to the Atronomer's almanach is the difference between
# the Julian date and JD 2451545.0 (noon, 1 January 2000)
time <- jd - 51545.
# Ecliptic coordinates
# Mean longitude
mnlong <- 280.460 + .9856474 * time
mnlong <- mnlong %% 360
mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360
# Mean anomaly
mnanom <- 357.528 + .9856003 * time
mnanom <- mnanom %% 360
mnanom[mnanom < 0] <- mnanom[mnanom < 0] + 360
mnanom <- mnanom * deg2rad
# Ecliptic longitude and obliquity of ecliptic
eclong <- mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)
eclong <- eclong %% 360
eclong[eclong < 0] <- eclong[eclong < 0] + 360
oblqec <- 23.429 - 0.0000004 * time
eclong <- eclong * deg2rad
oblqec <- oblqec * deg2rad
# Celestial coordinates
# Right ascension and declination
num <- cos(oblqec) * sin(eclong)
den <- cos(eclong)
ra <- atan(num / den)
ra[den < 0] <- ra[den < 0] + pi
ra[den >= 0 & num < 0] <- ra[den >= 0 & num < 0] + twopi
dec <- asin(sin(oblqec) * sin(eclong))
# Local coordinates
# Greenwich mean sidereal time
gmst <- 6.697375 + .0657098242 * time + hour
gmst <- gmst %% 24
gmst[gmst < 0] <- gmst[gmst < 0] + 24.
# Local mean sidereal time
lmst <- gmst + long / 15.
lmst <- lmst %% 24.
lmst[lmst < 0] <- lmst[lmst < 0] + 24.
lmst <- lmst * 15. * deg2rad
# Hour angle
ha <- lmst - ra
ha[ha < -pi] <- ha[ha < -pi] + twopi
ha[ha > pi] <- ha[ha > pi] - twopi
# Latitude to radians
lat <- lat * deg2rad
# Azimuth and elevation
el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
az <- asin(-cos(dec) * sin(ha) / cos(el))
elc <- asin(sin(dec) / sin(lat))
az[el >= elc] <- pi - az[el >= elc]
az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi
el <- el / deg2rad
az <- az / deg2rad
lat <- lat / deg2rad
return(list(elevation=el, azimuth=az))
}
Problem polega na tym, że azymut, który zwraca, wydaje się zły. Na przykład, jeśli uruchomię funkcję podczas przesilenia letniego (południowego) o godzinie 12:00 dla lokalizacji 0ºE i 41ºS, 3ºS, 3ºN i 41ºN:
> sunPosition(2012,12,22,12,0,0,-41,0)
$elevation
[1] 72.42113
$azimuth
[1] 180.9211
> sunPosition(2012,12,22,12,0,0,-3,0)
$elevation
[1] 69.57493
$azimuth
[1] -0.79713
Warning message:
In asin(sin(dec)/sin(lat)) : NaNs produced
> sunPosition(2012,12,22,12,0,0,3,0)
$elevation
[1] 63.57538
$azimuth
[1] -0.6250971
Warning message:
In asin(sin(dec)/sin(lat)) : NaNs produced
> sunPosition(2012,12,22,12,0,0,41,0)
$elevation
[1] 25.57642
$azimuth
[1] 180.3084
Te liczby po prostu nie wydają się prawidłowe. Wysokość, z której jestem zadowolony - pierwsze dwie powinny być mniej więcej takie same, trzecia odrobinę niższa, a czwarta znacznie niższa. Jednak pierwszy azymut powinien być mniej więcej na północ, podczas gdy liczba, którą podaje, jest całkowitym przeciwieństwem. Pozostałe trzy powinny wskazywać mniej więcej dokładnie na południe, ale tylko ostatnia tak robi. Dwie w środkowym punkcie tuż przy północy, znowu 180º na zewnątrz.
Jak widać, istnieje również kilka błędów wywoływanych przy niskich szerokościach geograficznych (blisko równika)
Uważam, że usterka jest w tej sekcji, a błąd jest wyzwalany w trzeciej linii (zaczynając od elc).
# Azimuth and elevation
el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
az <- asin(-cos(dec) * sin(ha) / cos(el))
elc <- asin(sin(dec) / sin(lat))
az[el >= elc] <- pi - az[el >= elc]
az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi
Przeszukałem go i znalazłem podobny fragment kodu w C, przekonwertowany na R, linia, której używa do obliczenia azymutu, wyglądałaby tak
az <- atan(sin(ha) / (cos(ha) * sin(lat) - tan(dec) * cos(lat)))
Wydaje się, że dane wyjściowe zmierzają we właściwym kierunku, ale po prostu nie mogę ich uzyskać, aby dać mi właściwą odpowiedź przez cały czas, gdy jest konwertowany z powrotem na stopnie.
Korekta kodu (podejrzewam, że to tylko kilka linii powyżej), aby obliczyć poprawny azymut, byłaby fantastyczna.
