Это кажется важной темой, поэтому я написал более длинный, чем обычно, ответ: если этот алгоритм будет использоваться другими в будущем, я думаю, что важно, чтобы он сопровождался ссылками на литературу, из которой он был получен.
Краткий ответ
Как вы заметили, ваш опубликованный код не работает должным образом для мест вблизи экватора или в южном полушарии.
Чтобы исправить это, просто замените эти строки в исходном коде:
elc <- asin(sin(dec) / sin(lat))
az[el >= elc] <- pi - az[el >= elc]
az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi
с этими:
cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
sinAzNeg <- (sin(az) < 0)
az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
az[!cosAzPos] <- pi - az[!cosAzPos]
Теперь он должен работать для любой точки земного шара.
Обсуждение
Код в вашем примере адаптирован почти дословно из статьи Дж.Дж. 1988 года. Михальский (Солнечная энергия. 40: 227-235). Эта статья, в свою очередь, усовершенствовала алгоритм, представленный в статье Р. Уолравена в 1978 году (Solar Energy. 20: 393-397). Уолравен сообщил, что этот метод успешно использовался в течение нескольких лет для точного позиционирования поляризационного радиометра в Дэвисе, Калифорния (38 ° 33 '14 "с.ш., 121 ° 44' 17" з.д.).
Код как Михальского, так и Вальравена содержит важные / фатальные ошибки. В частности, хотя алгоритм Михальского прекрасно работает в большинстве Соединенных Штатов, он выходит из строя (как вы обнаружили) для областей вблизи экватора или в южном полушарии. В 1989 году Ю.В. Спенсер из Виктории, Австралия, отметил то же самое (Solar Energy. 42 (4): 353):
Уважаемый сэр:
Метод Михальского для назначения вычисленного азимута правильному квадранту, полученному из Уолравена, не дает правильных значений при применении для южных (отрицательных) широт. Кроме того, расчет критической высоты (elc) потерпит неудачу на нулевой широте из-за деления на ноль. Оба эти возражения можно избежать, просто назначив азимут правильному квадранту, учитывая знак cos (азимут).
Мои изменения в вашем коде основаны на исправлениях, предложенных Спенсером в опубликованном комментарии. Я просто несколько изменил их, чтобы гарантировать, что функция R sunPosition() остается «векторизованной» (то есть работает должным образом на векторах местоположений точек, вместо того, чтобы проходить по одной точке за раз).
Точность функции sunPosition()
Чтобы проверить, что sunPosition() работает правильно, я сравнил его результаты с результатами, рассчитанными солнечным калькулятором Национального управления океанических и атмосферных исследований . В обоих случаях положения солнца были рассчитаны для полудня (12:00 вечера) южного летнего солнцестояния (22 декабря) 2012 года. Все результаты были в пределах 0,02 градуса.
testPts <- data.frame(lat = c(-41,-3,3, 41),
long = c(0, 0, 0, 0))
# Sun's position as returned by the NOAA Solar Calculator,
NOAA <- data.frame(elevNOAA = c(72.44, 69.57, 63.57, 25.6),
azNOAA = c(359.09, 180.79, 180.62, 180.3))
# Sun's position as returned by sunPosition()
sunPos <- sunPosition(year = 2012,
month = 12,
day = 22,
hour = 12,
min = 0,
sec = 0,
lat = testPts$lat,
long = testPts$long)
cbind(testPts, NOAA, sunPos)
# lat long elevNOAA azNOAA elevation azimuth
# 1 -41 0 72.44 359.09 72.43112 359.0787
# 2 -3 0 69.57 180.79 69.56493 180.7965
# 3 3 0 63.57 180.62 63.56539 180.6247
# 4 41 0 25.60 180.30 25.56642 180.3083
Другие ошибки в коде
В опубликованном коде есть как минимум две другие (довольно незначительные) ошибки. Первое приводит к тому, что 29 февраля и 1 марта високосных лет считаются 61-м днем года. Вторая ошибка происходит от опечатки в оригинальной статье, которая была исправлена Михальским в заметке 1989 года (Solar Energy. 43 (5): 323).
Этот блок кода показывает ошибочные строки, закомментированные и сопровождаемые исправленными версиями:
# leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & day >= 60
leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) &
day >= 60 & !(month==2 & day==60)
# oblqec <- 23.429 - 0.0000004 * time
oblqec <- 23.439 - 0.0000004 * time
Исправленная версия sunPosition()
Вот исправленный код, который был проверен выше:
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 & !(month==2 & 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.439 - 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))
# For logic and names, see Spencer, J.W. 1989. Solar Energy. 42(4):353
cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
sinAzNeg <- (sin(az) < 0)
az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
az[!cosAzPos] <- pi - az[!cosAzPos]
# if (0 < sin(dec) - sin(el) * sin(lat)) {
# if(sin(az) < 0) az <- az + twopi
# } else {
# az <- pi - az
# }
el <- el / deg2rad
az <- az / deg2rad
lat <- lat / deg2rad
return(list(elevation=el, azimuth=az))
}
Ссылки:
Михальский, J.J. 1988. Алгоритм астрономического альманаха для приблизительного положения Солнца (1950-2050). Солнечная энергия. 40 (3):. 227-235 * +1056 *
Михальский, J.J. 1989. Ошибки. Солнечная энергия. 43 (5): 323.
Спенсер, J.W. 1989. Комментарии к «Алгоритму астрономического альманаха для приблизительного положения Солнца (1950-2050)». Солнечная энергия. 42 (4): 353.
Walraven, R. 1978. Расчет положения солнца. Солнечная энергия. 20: 393-397
.