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

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

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

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

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]

Теперь он должен работать для любого места на земном шаре.

Обсуждение

Код в вашем примере практически дословно адаптирован из статьи J.J. Михальский (Солнечная энергия. 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.

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