Этот вопрос был задан до того чуть более трех лет назад. Был дан ответ, однако я обнаружил глюк в решении.

Код ниже находится на 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)))

Кажется, что вывод здесь движется в правильном направлении, но я просто не могу заставить его давать мне правильный ответ все время, когда он конвертируется обратно в градусы.

Исправление кода (подозреваю, что это всего несколько строк выше), чтобы он вычислял правильный азимут, было бы фантастическим.

Это кажется важной темой, поэтому я опубликовал более длинный, чем обычно, ответ: если этот алгоритм будет использоваться другими в будущем, я думаю, что важно, чтобы он сопровождался ссылками на литературу, из которой он был получен .

Краткий ответ

Как вы отметили, ваш опубликованный код не работает должным образом для местоположений вблизи экватора или в южном полушарии.

Чтобы исправить это, просто замените эти строки в исходном коде:

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 г. (Solar Energy. 40: 227-235). Эта статья, в свою очередь, усовершенствовала алгоритм, представленный в статье Р. Вальравена 1978 года (Solar Energy. 20: 393-397). Вальравен сообщил, что этот метод успешно использовался в течение нескольких лет для точного позиционирования поляризационного радиометра в Дэвисе, Калифорния (38 ° 33 '14 "с.ш., 121 ° 44' 17" з.д.).

И код Михальского, и код Валравена содержат важные / фатальные ошибки. В частности, хотя алгоритм Михальского отлично работает на большей части территории Соединенных Штатов, он не работает (как вы обнаружили) для областей вблизи экватора или в южном полушарии. В 1989 году Дж. В. Спенсер из Виктории, Австралия, отметил то же самое (Solar Energy. 42 (4): 353):

Уважаемый господин:

Метод Михальского для отнесения вычисленного азимута к правильному квадранту, полученный из Walraven, не дает правильных значений при применении для южных (отрицательных) широт. Кроме того, вычисление критической высоты (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))
}

Рекомендации:

Михальский, Дж. Дж. 1988. Алгоритм Астрономического Альманаха для приблизительного положения Солнца (1950-2050). Солнечная энергия. 40 (3): 227-235.

Михальский, JJ 1989. Errata. Солнечная энергия. 43 (5): 323.

Спенсер, Дж. У. 1989. Комментарии к "Алгоритму приблизительного положения Солнца в астрономическом альманахе (1950-2050)". Солнечная энергия. 42 (4): 353.

Вальравен Р. 1978. Расчет положения солнца. Солнечная энергия. 20: 393-397.

Используя "NOAA Solar Calculations" по одной из приведенных выше ссылок, я немного изменил заключительную часть функции, используя немного другой алгоритм, который, я надеюсь, был переведен без ошибок. Я закомментировал бесполезный код и добавил новый алгоритм сразу после преобразования широты в радианы:

# -----------------------------------------------
# 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 знаков после запятой зафиксировало данные в моем наборе данных.