Положение солнца с учетом времени суток, широты и долготы

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

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

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

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