Esta pregunta se ha hecho antes hace poco más de tres años. Se dio una respuesta, sin embargo, encontré un problema técnico en la solución.

El código a continuación está en R. Lo porté a otro idioma, sin embargo, probé el código original directamente en R para asegurarme de que el problema no estaba en mi portabilidad.

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))
}

El problema al que me golpeo es que el azimut que devuelve parece incorrecto. Por ejemplo, si ejecuto la función en el solsticio de verano (sur) a las 12:00 para las ubicaciones 0ºE y 41ºS, 3ºS, 3ºN y 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

Estos números simplemente no parecen correctos. La elevación con la que estoy contento: las dos primeras deberían ser aproximadamente iguales, la tercera un poco más baja y la cuarta mucho más baja. Sin embargo, el primer acimut debería estar aproximadamente al norte, mientras que el número que da es todo lo contrario. Los tres restantes deben apuntar aproximadamente hacia el sur, sin embargo, solo el último lo hace. Los dos en el punto medio justo al norte, nuevamente 180º hacia afuera.

Como puede ver, también hay un par de errores que se activan con las latitudes bajas (cierre el ecuador)

Creo que la falla está en esta sección, y el error se activa en la tercera línea (comenzando con 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

Busqué en Google y encontré un fragmento de código similar en C, convertido a R, la línea que usa para calcular el acimut sería algo así como

az <- atan(sin(ha) / (cos(ha) * sin(lat) - tan(dec) * cos(lat)))

La salida aquí parece dirigirse en la dirección correcta, pero no puedo hacer que me dé la respuesta correcta todo el tiempo cuando se convierte de nuevo a grados.

Una corrección del código (sospecho que son solo las pocas líneas de arriba) para que calcule el acimut correcto sería fantástico.

Este parece un tema importante, por lo que he publicado una respuesta más larga de lo habitual: si este algoritmo va a ser utilizado por otros en el futuro, creo que es importante que vaya acompañado de referencias a la literatura de la que se ha derivado. .

La respuesta corta

Como ha observado, su código publicado no funciona correctamente para ubicaciones cerca del ecuador o en el hemisferio sur.

Para solucionarlo, simplemente reemplace estas líneas en su código original:

elc <- asin(sin(dec) / sin(lat))
az[el >= elc] <- pi - az[el >= elc]
az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi

con estos:

cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
sinAzNeg <- (sin(az) < 0)
az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
az[!cosAzPos] <- pi - az[!cosAzPos]

Ahora debería funcionar para cualquier lugar del mundo.

Discusión

El código de su ejemplo está adaptado casi literalmente de un artículo de 1988 de JJ Michalsky (Solar Energy. 40: 227-235). Ese artículo a su vez refinó un algoritmo presentado en un artículo de 1978 por R. Walraven (Solar Energy. 20: 393-397). Walraven informó que el método se había utilizado con éxito durante varios años para colocar con precisión un radiómetro polarizador en Davis, California (38 ° 33 '14 "N, 121 ° 44' 17" W).

Tanto el código de Michalsky como el de Walraven contienen errores importantes / fatales. En particular, si bien el algoritmo de Michalsky funciona bien en la mayor parte de los Estados Unidos, falla (como ha descubierto) en áreas cercanas al ecuador o en el hemisferio sur. En 1989, JW Spencer de Victoria, Australia, señaló lo mismo (Solar Energy. 42 (4): 353):

Estimado señor:

El método de Michalsky para asignar el acimut calculado al cuadrante correcto, derivado de Walraven, no da valores correctos cuando se aplica a latitudes meridionales (negativas). Además, el cálculo de la elevación crítica (elc) fallará para una latitud de cero debido a la división por cero. Ambas objeciones pueden evitarse simplemente asignando el acimut al cuadrante correcto considerando el signo de cos (acimut).

Mis ediciones de su código se basan en las correcciones sugeridas por Spencer en ese comentario publicado. Simplemente los he alterado un poco para asegurarme de que la función R sunPosition()permanezca 'vectorizada' (es decir, funcionando correctamente en vectores de ubicaciones de puntos, en lugar de tener que pasar un punto a la vez).

Precisión de la función sunPosition()

Para probar que sunPosition()funciona correctamente, he comparado sus resultados con los calculados por la Calculadora Solar de la Administración Nacional Oceánica y Atmosférica . En ambos casos, las posiciones del sol se calcularon para el mediodía (12:00 PM) en el solsticio de verano del sur (22 de diciembre) de 2012. Todos los resultados estuvieron de acuerdo dentro de 0.02 grados.

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

Otros errores en el código

Hay al menos otros dos errores (bastante menores) en el código publicado. La primera causa que el 29 de febrero y el 1 de marzo de los años bisiestos se cuenten como el día 61 del año. El segundo error deriva de un error tipográfico en el artículo original, que fue corregido por Michalsky en una nota de 1989 (Solar Energy. 43 (5): 323).

Este bloque de código muestra las líneas ofensivas, comentadas y seguidas inmediatamente por versiones corregidas:

# 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

Versión corregida de sunPosition()

Aquí está el código corregido que se verificó anteriormente:

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))
}

Referencias:

Michalsky, JJ 1988. El algoritmo del Almanaque Astronómico para la posición solar aproximada (1950-2050). Energía solar. 40 (3): 227-235.

Michalsky, JJ 1989. Errata. Energía solar. 43 (5): 323.

Spencer, JW 1989. Comentarios sobre "El algoritmo del almanaque astronómico para la posición solar aproximada (1950-2050)". Energía solar. 42 (4): 353.

Walraven, R. 1978. Calculando la posición del sol. Energía solar. 20: 393-397.

Usando "Cálculos solares NOAA" de uno de los enlaces de arriba he cambiado un poco la parte final de la función usando un algoritmo ligeramente diferente que, espero, se haya traducido sin errores. He comentado el código ahora inútil y agregué el nuevo algoritmo justo después de la conversión de latitud a radianes:

# -----------------------------------------------
# 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))

Para verificar la tendencia de acimut en los cuatro casos que mencionó, grafiquemos la hora del día:

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)

Imagen adjunta:

Aquí hay una reescritura que es más idiomática para R y más fácil de depurar y mantener. Es esencialmente la respuesta de Josh, pero con el acimut calculado usando los algoritmos de Josh y Charlie para comparar. También incluí las simplificaciones del código de fecha de mi otra respuesta. El principio básico era dividir el código en muchas funciones más pequeñas para las que se pueden escribir pruebas unitarias más fácilmente.

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)
  )
}

Esta es una actualización sugerida de la excelente respuesta de Josh.

Gran parte del inicio de la función es un código repetitivo para calcular el número de días desde el mediodía del 1 de enero de 2000. Esto se trata mucho mejor con la función de fecha y hora existente de R.

También creo que en lugar de tener seis variables diferentes para especificar la fecha y la hora, es más fácil (y más consistente con otras funciones de R) especificar un objeto de fecha existente o cadenas de fecha + cadenas de formato.

Aquí hay dos funciones auxiliares

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)
}

Y el inicio de la función ahora se simplifica a

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)
  #...

La otra rareza está en las líneas como

mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360

Dado que mnlongha %%invocado sus valores, todos deberían ser no negativos, por lo que esta línea es superflua.

Necesitaba la posición del sol en un proyecto de Python. Adapté el algoritmo de Josh O'Brien.

Gracias Josh.

Por si pudiera ser de utilidad para alguien, aquí está mi adaptación.

Tenga en cuenta que mi proyecto solo necesitaba la posición instantánea del sol, por lo que el tiempo no es un parámetro.

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

Encontré un pequeño problema con un punto de datos y las funciones de Richie Cotton anteriores (en la implementación del código de Charlie)

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

porque en la función solarAzimuthRadiansCharlie ha habido una excitación de punto flotante alrededor de un ángulo de 180, que (sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle))es la cantidad más pequeña sobre 1, 1.0000000000000004440892098, lo que genera un NaN ya que la entrada de acos no debe estar por encima de 1 ni por debajo de -1.

Sospecho que puede haber casos extremos similares para el cálculo de Josh, donde los efectos de redondeo de punto flotante hacen que la entrada para el paso asin esté fuera de -1: 1 pero no los he alcanzado en mi conjunto de datos particular.

En la media docena de casos en los que he llegado a este punto, el "verdadero" (en medio del día o de la noche) es cuando el problema ocurre, de manera que empíricamente el valor verdadero debería ser 1 / -1. Por esa razón, me sentiría cómodo arreglando eso aplicando un paso de redondeo dentro de solarAzimuthRadiansJoshy solarAzimuthRadiansCharlie. No estoy seguro de cuál es la precisión teórica del algoritmo NOAA (el punto en el que la precisión numérica deja de importar de todos modos), pero redondear a 12 lugares decimales fijó los datos en mi conjunto de datos.