Posición del sol dada la hora del día, latitud y longitud

83

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.

CucharaNZ
fuente
2
Es posible que tenga más suerte en el intercambio de pila de matemáticas
abcde123483
1
Hay un código para hacer esto en el paquete maptools, ver? Solarpos
mdsumner
Gracias @ulvund, podría intentarlo a continuación.
SpoonNZ
4
Bien, entonces creo que deberías copiar el Javascript del sitio de la NOAA, esa es la fuente de muchas versiones. El código que escribimos colapsó todo esto en lo que necesitábamos en dos funciones más pequeñas, pero eso fue solo para elevación y sintonizado con una aplicación en particular. Solo vea la fuente de srrb.noaa.gov/highlights/sunrise/azel.html
mdsumner
1
¿Has probado mi respuesta a la pregunta anterior ? ephemincluso podría tener en cuenta la refracción de la atmósfera (influenciada por la temperatura, la presión) y la elevación de un observador.
jfs

Respuestas:

110

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.

Josh O'Brien
fuente
¡Gracias por la fantástica respuesta! No he estado aquí durante el fin de semana, así que lo perdí, lo siento. No tendré la oportunidad de probar esto hasta al menos esta noche, pero parece que funcionará. ¡Salud!
SpoonNZ
1
@SpoonNZ - Es un placer. Si necesita copias en pdf de cualquiera de esas referencias citadas, avíseme a mi dirección de correo electrónico y se las puedo enviar.
Josh O'Brien
1
@ JoshO'Brien: Acabo de agregar algunas sugerencias en una respuesta separada. Es posible que desee echar un vistazo e incorporarlos a los suyos.
Richie Cotton
@RichieCotton - Gracias por publicar sus sugerencias. No los voy a agregar aquí, pero solo porque son Respecíficos y el OP estaba usando el código R para intentar depurarlo antes de portarlo a otro idioma. (De hecho, acabo de editar mi publicación para corregir un error de procesamiento de fecha en el código original, y es exactamente el tipo de error que justifica el uso de código de nivel superior como el que propuso). ¡Salud!
Josh O'Brien
También se podrían combinar las fechas julianas para: hora = 365 * (año - 2000) + piso ((año - 1949) / 4) + día + hora - 13.5
Hawk
19

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:

ingrese la descripción de la imagen aquí

mbask
fuente
@Josh O'Brien: Su respuesta muy detallada es una lectura excelente. Como nota relacionada, nuestras funciones SunPosition producen exactamente los mismos resultados.
mbask
Adjunto el archivo de imagen, si lo desea.
mdsumner
1
@Charlie - Gran respuesta, y las parcelas son una adición especialmente agradable. Antes de verlos, no me había dado cuenta de lo diferentes que serían las coordenadas azimutales nocturnas del sol en ubicaciones 'ecuatoriales' frente a ubicaciones más 'templadas'. Realmente genial.
Josh O'Brien
12

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)
  )
}
4 revoluciones
fuente
Tenga en cuenta que cuando realice pruebas en el sitio web de la NOAA aquí: esrl.noaa.gov/gmd/grad/solcalc/azel.html que NOAA usa Longitude West como + ve. Este algoritmo usa Longitude West como -ve.
Neon22
Cuando ejecuto "sunPosition (lat = 43, long = -89)", obtengo una elevación de 52 y un acimut de 175. Pero al usar la aplicación web de la NOAA esrl.noaa.gov/gmd/grad/solcalc , obtengo una elevación de alrededor de 5 y acimut de 272. ¿Me falta algo? NOAA es correcta, pero no puedo obtener sunPosition para dar resultados precisos.
Tedward
@Tedward utiliza de forma sunPositionpredeterminada la fecha y la hora actuales. ¿Es eso lo que querías?
Richie Cotton
Si. También probé con diferentes tiempos. Esto fue al final del día, voy a intentarlo de nuevo hoy con un nuevo comienzo. Estoy bastante seguro de que estoy haciendo algo mal, pero no sé qué. Seguiré trabajando en eso.
Tedward
Necesitaba convertir "cuándo" a UTC para obtener resultados precisos. Consulte stackoverflow.com/questions/39393514/… . @aichao sugiere código para la conversión.
Tedward
10

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.

Richie algodón
fuente
¡Muchas gracias! Como se mencionó, he transferido esto a PHP (y probablemente iré a Javascript, solo tengo que decidir dónde quiero qué funciones se manejan) para que el código no me ayude mucho, pero debería poder ser transferido (aunque con un poco ¡Más pensamiento involucrado que con el código original!). Necesito modificar un poco el código que maneja las zonas horarias, por lo que podría integrar este cambio al mismo tiempo.
SpoonNZ
2
Cambios ingeniosos @Richie Cotton. Tenga en cuenta que la hora de asignación <- hora_de_día debería ser realmente la hora <- hora_de_día (cuándo) y que el tiempo variable debería contener el número de días, no un objeto de clase "difftime". La segunda línea de la función astronomers_almanac_time debe cambiarse a algo como as.numeric (difftime (x, origin, units = "days"), units = "days").
mbask
1
Gracias por las excelentes sugerencias. Podría ser bueno (si está interesado) incluir en su publicación una versión editada de toda la sunPosition()función que sea más R-ish en su construcción.
Josh O'Brien
@ JoshO'Brien: Hecho. Hice la wiki de la comunidad de respuestas, ya que es una combinación de todas nuestras respuestas. Da la misma respuesta que la suya para la hora actual y las coordenadas predeterminadas (¿suizas?), Pero se necesitan muchas más pruebas.
Richie Cotton
@RichieCotton - Qué buena idea. Analizaré más a fondo lo que has hecho tan pronto como tenga la oportunidad.
Josh O'Brien
4

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
Jérôme
fuente
Esto fue realmente útil para mí. Gracias. Una cosa que hice fue agregar un ajuste para el horario de verano. En caso de que sea útil, fue simplemente: if (time.localtime (). Tm_isdst == 1): hora + = 1
Mark Ireland
1

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.

David Hood
fuente