Положение солнца с учетом времени суток, широты и долготы
Этот вопрос был задан чуть более трех лет назад. Был дан ответ, однако я нашел ошибку в решении.
Код ниже написан на R. Я перенес его на другой язык, однако я проверил исходный код непосредственно на R, чтобы убедиться, что проблема не с моим портированием.
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))
}
Проблема, которую я поражаю, состоит в том, что азимут, который это возвращает, кажется неправильным. Например, если я запускаю функцию в (южное) летнее солнцестояние в 12:00 для мест 0ºE и 41ºS, 3ºS, 3ºN и 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
Эти цифры просто не кажутся правильными. Возвышение, которым я доволен - первые два должны быть примерно одинаковыми, третье - чуть ниже, а четвертое - намного ниже. Однако первый азимут должен быть примерно на север, тогда как число, которое он дает, является полной противоположностью. Остальные три должны указывать примерно на юг, однако только последний делает. Два в средней точке рядом с севером, снова на 180º.
Как вы можете видеть, есть также несколько ошибок, вызванных низкими широтами (закройте экватор)
Я полагаю, что ошибка в этом разделе, с ошибкой, вызванной в третьей строке (начиная с 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
Я гуглил и нашел похожий кусок кода в C, преобразовал в R строку, которую он использует для вычисления азимута, что-то вроде
az <- atan(sin(ha) / (cos(ha) * sin(lat) - tan(dec) * cos(lat)))
Вывод здесь, кажется, идет в правильном направлении, но я просто не могу получить его, чтобы дать мне правильный ответ все время, когда он переводится обратно в градусы.
Исправление кода (подозреваю, что это всего лишь несколько строк выше), чтобы заставить его вычислять правильный азимут, было бы фантастическим.
6 ответов
Это кажется важной темой, поэтому я разместил более длинный, чем обычно, ответ: если этот алгоритм будет использоваться другими в будущем, я думаю, важно, чтобы он сопровождался ссылками на литературу, из которой он был получен,
Краткий ответ
Как вы уже заметили, ваш опубликованный код не работает должным образом для мест вблизи экватора или в южном полушарии.
Чтобы это исправить, просто замените эти строки в исходном коде:
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 года JJ Michalsky (Solar Energy. 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))
}
Рекомендации:
Михальский, JJ 1988. Алгоритм астрономического альманаха для приближенного положения Солнца (1950-2050). Солнечная энергия. 40(3):227-235.
Михальский, JJ 1989. Ошибки. Солнечная энергия. 43(5):323.
Спенсер, JW 1989. Комментарии к "Алгоритму астрономического альманаха для приблизительного положения Солнца (1950-2050)". Солнечная энергия. 42 (4): 353.
Walraven, R. 1978. Расчет положения солнца. Солнечная энергия. 20:393-397.
Используя "Солнечные вычисления NOAA" по одной из ссылок выше, я немного изменил финальную часть функции, используя немного другой алгоритм, который, я надеюсь, был переведен без ошибок. Я закомментировал бесполезный теперь код и добавил новый алгоритм сразу после преобразования широты в радианы:
# -----------------------------------------------
# New code
# Solar zenith angle
zenithAngle <- acos(sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(ha))
# Solar azimuth
az <- acos(((sin(lat) * cos(zenithAngle)) - sin(dec)) / (cos(lat) * sin(zenithAngle)))
rm(zenithAngle)
# -----------------------------------------------
# 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
# -----------------------------------------------
# New code
if (ha > 0) az <- az + 180 else az <- 540 - az
az <- az %% 360
# -----------------------------------------------
return(list(elevation=el, azimuth=az))
Чтобы проверить азимутальный тренд в четырех упомянутых вами случаях, построим график в зависимости от времени суток:
hour <- seq(from = 0, to = 23, by = 0.5)
azimuth <- data.frame(hour = hour)
az41S <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,-41,0)$azimuth)
az03S <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,-03,0)$azimuth)
az03N <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,03,0)$azimuth)
az41N <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,41,0)$azimuth)
azimuth <- cbind(azimuth, az41S, az03S, az41N, az03N)
rm(az41S, az03S, az41N, az03N)
library(ggplot2)
azimuth.plot <- melt(data = azimuth, id.vars = "hour")
ggplot(aes(x = hour, y = value, color = variable), data = azimuth.plot) +
geom_line(size = 2) +
geom_vline(xintercept = 12) +
facet_wrap(~ variable)
Изображение прикреплено:
Вот переписывание, которое более идиоматично для R, и его легче отлаживать и поддерживать. По сути, это ответ Джоша, но с азимутом, рассчитанным с использованием алгоритмов Джоша и Чарли для сравнения. Я также включил упрощения в код даты из моего другого ответа. Основной принцип состоял в том, чтобы разбить код на множество мелких функций, для которых вам будет проще написать модульные тесты.
astronomersAlmanacTime <- function(x)
{
# Astronomer's almanach time is the number of
# days since (noon, 1 January 2000)
origin <- as.POSIXct("2000-01-01 12:00:00")
as.numeric(difftime(x, origin, units = "days"))
}
hourOfDay <- function(x)
{
x <- as.POSIXlt(x)
with(x, hour + min / 60 + sec / 3600)
}
degreesToRadians <- function(degrees)
{
degrees * pi / 180
}
radiansToDegrees <- function(radians)
{
radians * 180 / pi
}
meanLongitudeDegrees <- function(time)
{
(280.460 + 0.9856474 * time) %% 360
}
meanAnomalyRadians <- function(time)
{
degreesToRadians((357.528 + 0.9856003 * time) %% 360)
}
eclipticLongitudeRadians <- function(mnlong, mnanom)
{
degreesToRadians(
(mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)) %% 360
)
}
eclipticObliquityRadians <- function(time)
{
degreesToRadians(23.439 - 0.0000004 * time)
}
rightAscensionRadians <- function(oblqec, eclong)
{
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] + 2 * pi
ra
}
rightDeclinationRadians <- function(oblqec, eclong)
{
asin(sin(oblqec) * sin(eclong))
}
greenwichMeanSiderealTimeHours <- function(time, hour)
{
(6.697375 + 0.0657098242 * time + hour) %% 24
}
localMeanSiderealTimeRadians <- function(gmst, long)
{
degreesToRadians(15 * ((gmst + long / 15) %% 24))
}
hourAngleRadians <- function(lmst, ra)
{
((lmst - ra + pi) %% (2 * pi)) - pi
}
elevationRadians <- function(lat, dec, ha)
{
asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
}
solarAzimuthRadiansJosh <- function(lat, dec, ha, el)
{
az <- asin(-cos(dec) * sin(ha) / cos(el))
cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
sinAzNeg <- (sin(az) < 0)
az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + 2 * pi
az[!cosAzPos] <- pi - az[!cosAzPos]
az
}
solarAzimuthRadiansCharlie <- function(lat, dec, ha)
{
zenithAngle <- acos(sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(ha))
az <- acos((sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle)))
ifelse(ha > 0, az + pi, 3 * pi - az) %% (2 * pi)
}
sunPosition <- function(when = Sys.time(), format, lat = 46.5, long = 6.5)
{
if(is.character(when)) when <- strptime(when, format)
when <- lubridate::with_tz(when, "UTC")
time <- astronomersAlmanacTime(when)
hour <- hourOfDay(when)
# Ecliptic coordinates
mnlong <- meanLongitudeDegrees(time)
mnanom <- meanAnomalyRadians(time)
eclong <- eclipticLongitudeRadians(mnlong, mnanom)
oblqec <- eclipticObliquityRadians(time)
# Celestial coordinates
ra <- rightAscensionRadians(oblqec, eclong)
dec <- rightDeclinationRadians(oblqec, eclong)
# Local coordinates
gmst <- greenwichMeanSiderealTimeHours(time, hour)
lmst <- localMeanSiderealTimeRadians(gmst, long)
# Hour angle
ha <- hourAngleRadians(lmst, ra)
# Latitude to radians
lat <- degreesToRadians(lat)
# Azimuth and elevation
el <- elevationRadians(lat, dec, ha)
azJ <- solarAzimuthRadiansJosh(lat, dec, ha, el)
azC <- solarAzimuthRadiansCharlie(lat, dec, ha)
data.frame(
elevation = radiansToDegrees(el),
azimuthJ = radiansToDegrees(azJ),
azimuthC = radiansToDegrees(azC)
)
}
Это рекомендуемое обновление отличного ответа Джоша.
Большая часть запуска функции - это стандартный код для расчета количества дней с полудня 1 января 2000 года. Это гораздо лучше, если использовать существующую функцию даты и времени R.
Я также думаю, что вместо шести различных переменных для указания даты и времени проще (и более согласованно с другими функциями R) указать существующий объект даты или строки даты + строки формата.
Вот две вспомогательные функции
astronomers_almanac_time <- function(x)
{
origin <- as.POSIXct("2000-01-01 12:00:00")
as.numeric(difftime(x, origin, units = "days"))
}
hour_of_day <- function(x)
{
x <- as.POSIXlt(x)
with(x, hour + min / 60 + sec / 3600)
}
И запуск функции теперь упрощается до
sunPosition <- function(when = Sys.time(), format, lat=46.5, long=6.5) {
twopi <- 2 * pi
deg2rad <- pi / 180
if(is.character(when)) when <- strptime(when, format)
time <- astronomers_almanac_time(when)
hour <- hour_of_day(when)
#...
Другая странность в линиях как
mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360
поскольку mnlong
имел %%
Обращаясь к его значениям, они все уже должны быть неотрицательными, так что эта строка лишняя.
Мне нужно было положение солнца в проекте Python. Я адаптировал алгоритм Джоша О'Брайена.
Спасибо, Джош.
На случай, если это кому-нибудь пригодится, вот моя адаптация.
Обратите внимание, что моему проекту требовалось только мгновенное положение солнца, поэтому время не является параметром.
def sunPosition(lat=46.5, long=6.5):
# Latitude [rad]
lat_rad = math.radians(lat)
# Get Julian date - 2400000
day = time.gmtime().tm_yday
hour = time.gmtime().tm_hour + \
time.gmtime().tm_min/60.0 + \
time.gmtime().tm_sec/3600.0
delta = time.gmtime().tm_year - 1949
leap = delta / 4
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)
t = jd - 51545
# Ecliptic coordinates
# Mean longitude
mnlong_deg = (280.460 + .9856474 * t) % 360
# Mean anomaly
mnanom_rad = math.radians((357.528 + .9856003 * t) % 360)
# Ecliptic longitude and obliquity of ecliptic
eclong = math.radians((mnlong_deg +
1.915 * math.sin(mnanom_rad) +
0.020 * math.sin(2 * mnanom_rad)
) % 360)
oblqec_rad = math.radians(23.439 - 0.0000004 * t)
# Celestial coordinates
# Right ascension and declination
num = math.cos(oblqec_rad) * math.sin(eclong)
den = math.cos(eclong)
ra_rad = math.atan(num / den)
if den < 0:
ra_rad = ra_rad + math.pi
elif num < 0:
ra_rad = ra_rad + 2 * math.pi
dec_rad = math.asin(math.sin(oblqec_rad) * math.sin(eclong))
# Local coordinates
# Greenwich mean sidereal time
gmst = (6.697375 + .0657098242 * t + hour) % 24
# Local mean sidereal time
lmst = (gmst + long / 15) % 24
lmst_rad = math.radians(15 * lmst)
# Hour angle (rad)
ha_rad = (lmst_rad - ra_rad) % (2 * math.pi)
# Elevation
el_rad = math.asin(
math.sin(dec_rad) * math.sin(lat_rad) + \
math.cos(dec_rad) * math.cos(lat_rad) * math.cos(ha_rad))
# Azimuth
az_rad = math.asin(
- math.cos(dec_rad) * math.sin(ha_rad) / math.cos(el_rad))
if (math.sin(dec_rad) - math.sin(el_rad) * math.sin(lat_rad) < 0):
az_rad = math.pi - az_rad
elif (math.sin(az_rad) < 0):
az_rad += 2 * math.pi
return el_rad, az_rad
Я столкнулся с небольшой проблемой с точкой данных и функциями Ричи Коттона выше (в реализации кода Чарли)
longitude= 176.0433687000000020361767383292317390441894531250
latitude= -39.173830619999996827118593500927090644836425781250
event_time = as.POSIXct("2013-10-24 12:00:00", format="%Y-%m-%d %H:%M:%S", tz = "UTC")
sunPosition(when=event_time, lat = latitude, long = longitude)
elevation azimuthJ azimuthC
1 -38.92275 180 NaN
Warning message:
In acos((sin(lat) * cos(zenithAngle) - sin(dec))/(cos(lat) * sin(zenithAngle))) : NaNs produced
потому что в функции solarAzimuthRadiansCharlie было возбуждение с плавающей запятой вокруг угла 180, так что (sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle))
самое маленькое количество больше 1, 1.0000000000000004440892098, которое генерирует NaN, так как входные данные для acos не должны быть выше 1 или ниже -1.
Я подозреваю, что могут быть аналогичные крайние случаи для вычисления Джоша, где эффекты округления с плавающей запятой приводят к тому, что вход для шага asin выходит за пределы -1: 1, но я не поразил их в моем конкретном наборе данных.
В полудюжине или около того случаев, которые я затронул, "истина" (середина дня или ночи) - это когда проблема возникает так эмпирически, истинное значение должно быть 1/-1. По этой причине мне было бы удобно исправить это, применив шаг округления внутри solarAzimuthRadiansJosh
а также solarAzimuthRadiansCharlie
, Я не уверен, какова теоретическая точность алгоритма NOAA (точка, в которой численная точность все равно перестает иметь значение), но округление до 12 десятичных знаков зафиксировало данные в моем наборе данных.