que - posicion del sol durante el dia
Posición del sol dada la hora del día, la latitud y la longitud (6)
Esta pregunta ha sido formulada hace poco más de tres años. Hubo una respuesta, sin embargo, encontré un problema técnico en la solución.
El siguiente código está en R. Lo he portado a otro idioma, sin embargo, he probado el código original directamente en R para asegurarme de que el problema no estuvo relacionado con mi portado.
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 que estoy afectando es que el acimut que devuelve parece estar mal. 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 más o menos las mismas, la tercera un toque más baja y la cuarta mucho más baja. Sin embargo, el primer azimut debe ser más o menos al norte, mientras que el número que da es el opuesto completo. Los tres restantes deben señalar más o menos al sur, sin embargo, solo el último lo hace. Los dos en el punto medio justo al norte, nuevamente a 180º.
Como puedes ver, también hay un par de errores disparados con las latitudes bajas (cerca del ecuador)
Creo que la falla está en esta sección, y el error se desencadena 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 similar de código en C, convertí a R la línea que usa para calcular el azimut sería algo así como
az <- atan(sin(ha) / (cos(ha) * sin(lat) - tan(dec) * cos(lat)))
La salida aquí parece estar yendo en la dirección correcta, pero simplemente no puedo obtener la respuesta correcta todo el tiempo cuando se convierte nuevamente en grados.
Una corrección del código (sospeche que son solo algunas líneas arriba) para hacer que calcule el acimut correcto sería fantástico.
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 azimut calculado usando los algoritmos de Josh y Charlie para comparar. También incluí las simplificaciones al 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 cuales puede escribir pruebas de unidades 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)
)
}
Encontré un pequeño problema con un punto de datos y las funciones de Richie Cotton arriba (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 excitación en coma flotante alrededor de un ángulo de 180 tal que (sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle)) es la cantidad más (sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle)) más de 1, 1.0000000000000004440892098, que genera un NaN ya que la entrada a acos no debe estar por encima de 1 o por debajo de -1.
Sospecho que podría haber casos de borde similares para el cálculo de Josh, donde los efectos de redondeo de coma flotante hacen que la entrada para el paso de entrada esté fuera de -1: 1 pero no los he golpeado en mi conjunto de datos particular.
En la media docena de casos que he encontrado, el "verdadero" (en medio del día o de la noche) es cuando ocurre el problema, así 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 solarAzimuthRadiansJosh y 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 el redondeo a 12 lugares decimales fijó los datos en mi conjunto de datos.
Esta es una actualización sugerida a la excelente respuesta de Josh.
Gran parte del inicio de la función es un código estándar para calcular el número de días desde el mediodía del 1 de enero de 2000. Esto se maneja mucho mejor utilizando 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 R) especificar un objeto de fecha existente o una cadena de fechas + cadenas de formato.
Aquí hay dos funciones de ayuda
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 comienzo de la función ahora se simplifica
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
Como mnlong ha llamado a %% en sus valores, todos deberían ser no negativos, por lo que esta línea es superflua.
Este parece ser un tema importante, así que he publicado una respuesta más larga que la típica: 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 habrás notado, tu código publicado no funciona correctamente para ubicaciones cercanas al 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 ubicación en el mundo.
Discusión
El código de su ejemplo se adapta casi textualmente a partir 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 de polarización en Davis, CA (38 ° 33 ''14 "N, 121 ° 44'' 17" W).
Tanto el código de Michalsky como el de Walraven contienen errores importantes / fatales. En particular, mientras que el algoritmo de Michalsky funciona bien en la mayoría de los Estados Unidos, falla (como lo ha encontrado) en áreas cercanas al ecuador o en el hemisferio sur. En 1989, JW Spencer de Victoria, Australia, notó lo mismo (Solar Energy 42 (4): 353):
Estimado señor:
El método de Michalsky para asignar el azimut 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 al considerar el signo de cos (acimut).
Mis modificaciones a su código se basan en las correcciones sugeridas por Spencer en el comentario publicado. Simplemente los he modificado un poco para asegurar que la función R sunPosition() permanezca ''vectorizada'' (es decir, funcione correctamente en vectores de ubicaciones de puntos, en lugar de tener que pasar un punto por 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 p.m.) en el solsticio de verano del sur (22 de diciembre), 2012. Todos los resultados estuvieron de acuerdo con 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. El primero causa que el 29 de febrero y el 1 de marzo de los años bisiestos se contabilicen como el día 61 del año. El segundo error se 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 las 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ó arriba:
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. Cálculo de la posición del sol. Energía solar. 20: 393-397.
Necesitaba la posición del sol en un proyecto de Python. Adapté el algoritmo de Josh O''Brien.
Gracias Josh.
En caso de que pueda ser útil para cualquier persona, aquí está mi adaptación.
Tenga en cuenta que mi proyecto solo necesitaba una posición de sol instantánea para que el tiempo no sea 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
Usando "NOAA Solar Calculations" de uno de los enlaces anteriores, he cambiado un poco la parte final de la función usando un algoritmo ligeramente diferente que, espero, 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 del azimut en los cuatro casos que mencionaste, tracémoslo en función de 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: