Position du soleil donnée heure, latitude et longitude
cette question a été posée avant il y a un peu plus de trois ans. Il y avait une réponse donnée, mais j'ai trouvé un bug dans la solution.
Code ci-dessous est en R. Je l'ai porté dans une autre langue, cependant j'ai testé le code original directement en R pour m'assurer que le problème ne se posait pas avec mon portage.
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))
}
le problème que je rencontre est que l'Azimut qu'il renvoie semble erroné. Par exemple, si je lance le fonction sur le solstice d'été (Sud) à 12h00 pour les emplacements 0ºE et 41ºS, 3ºS, 3ºN et 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
ces chiffres ne semblent pas justes. L'élévation dont je suis heureux - les deux premiers devraient être à peu près les mêmes, le troisième un peu plus bas, et le quatrième beaucoup plus bas. Cependant le premier Azimut doit être approximativement plein nord, tandis que le nombre qu'il donne est le contraire complet. Les trois autres devraient pointer vers le Sud, mais seul le dernier le fait. Les deux dans la pointe du milieu juste au nord, encore 180º dehors.
comme vous pouvez le voir Il ya aussi un couple d'erreurs déclenchées avec les basses latitudes (fermer l'Équateur)
je crois que la faute est dans cette section, l'erreur étant déclenchée à la troisième ligne (commençant par 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
j'ai cherché sur Google autour et trouvé un morceau de code similaire en C, converti en R la ligne qu'il utilise pour calculer l'Azimut serait être quelque chose comme
az <- atan(sin(ha) / (cos(ha) * sin(lat) - tan(dec) * cos(lat)))
la sortie ici semble se diriger dans la bonne direction, mais je ne peux pas l'obtenir pour me donner la bonne réponse tout le temps quand il est converti en degrés.
une correction du code (suspecte que ce ne sont que les quelques lignes ci-dessus) pour le faire calculer l'Azimut correct serait fantastique.
6 réponses
cela semble être un sujet important, donc j'ai posté une réponse plus longue que typique: si cet algorithme doit être utilisé par d'autres à l'avenir, je pense qu'il est important qu'il soit accompagné de références à la littérature à partir de laquelle il a été dérivé.
La réponse courte
Comme vous l'avez remarqué, votre posté code ne fonctionne pas correctement pour les emplacements près de l'équateur, ou dans l'hémisphère sud.
pour le fixer, il suffit de remplacer ces lignes dans votre code d'origine:
elc <- asin(sin(dec) / sin(lat))
az[el >= elc] <- pi - az[el >= elc]
az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi
avec ceci:
cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
sinAzNeg <- (sin(az) < 0)
az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
az[!cosAzPos] <- pi - az[!cosAzPos]
Il devrait maintenant fonctionner pour n'importe quel endroit sur le globe.
Discussion
le code dans votre exemple est adapté presque mot pour mot d'un article de 1988 par J. J. Michalsky (énergie solaire. 40: 227-235). Cet article affine à son tour un algorithme présenté dans un article de 1978 par R. Walraven (Solar Energy. 20: 393-397). Walraven a indiqué que la méthode avait été utilisée avec succès pendant plusieurs années pour positionner avec précision un radiomètre polarisant à Davis, CA (38° 33' 14" N, 121° 44' 17" O).
le code de Michalsky et de Walraven contient des erreurs importantes/fatales. en particulier, alors que L'algorithme de Michalsky fonctionne très bien dans la plupart des États-Unis, il échoue (comme vous l'avez trouvé) pour les zones proches de l'Équateur, ou dans l'hémisphère sud. En 1989, J. W. Spencer de Victoria, Australie, a noté la même chose (L'énergie solaire. 42(4):353):
Cher Monsieur:
la méthode de Michalsky pour attribuer l'Azimut calculé au quadrant correct, dérivé de Walraven, ne donne pas de valeurs correctes lorsqu'elle est appliquée aux latitudes méridionales (négatives). De plus, le calcul de l'élévation critique (elc) échouera pour une latitude de zéro à cause de la division par zéro. Ces deux objections peuvent être évitées simplement en assignant l'Azimut au quadrant correct en considérant le signe de cos (Azimut).
mes modifications à votre code sont basées sur les corrections suggérées par Spencer dans ce commentaire publié. Je les ai simplement quelque peu modifiés pour m'assurer que la fonction r sunPosition()
reste "vectorisée" (c'est-à-dire qu'elle fonctionne correctement sur les vecteurs de localisation des points, plutôt que d'avoir à passer un point à la fois).
précision de la fonction sunPosition()
pour tester que sunPosition()
fonctionne correctement, j'ai comparé ses résultats avec ceux calculés par la calculatrice solaire de la National Oceanic and Atmospheric Administration . Dans les deux cas, les positions du soleil ont été calculées pour Midi (12:00 PM) au solstice d'été du Sud (22 décembre 2012). Tous les résultats concordaient à 0,02 degré près.
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
autres erreurs dans le code
il y a au moins deux autres erreurs (mineures) dans le code affiché. Les premières causes le 29 février et le 1er mars des années bissextiles à tous les deux être comptés comme le jour 61 de l'année. La deuxième erreur provient d'une faute de frappe dans l'article original, qui a été corrigée par Michalsky dans une note de 1989 (énergie solaire. 43 (5): 323).
ce bloc de code affiche les lignes offensantes, commentées et suivies immédiatement des versions corrigées:
# 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
Version corrigée de sunPosition()
voici le code corrigé qui a été vérifié ci-dessus:
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))
}
, les Références:
Michalsky, J. J. 1988. L'Almanach Astronomique de l'algorithme approximatif de la position solaire (1950-2050). L'Énergie Solaire. 40 (3):227-235.
Michalsky, J. J. 1989. Errata. L'Énergie Solaire. 43 (5): 323.
Spencer, J. W. 1989. Commentaires sur " The Astronomical Almanach de l'Algorithme Approximatif de la Position Solaire (1950-2050)". L'Énergie Solaire. 42 (4): 353.
Walraven, R. 1978. Calcul de la position du soleil. L'Énergie Solaire. 20: 393-397.
en utilisant "NOAA Solar Calculations" de l'un des liens ci-dessus, j'ai changé un peu la partie finale de la fonction en utilisant un algorithme légèrement différent qui, je l'espère, ont traduit sans erreurs. J'ai commenté le code maintenant inutile et j'ai ajouté le nouvel algorithme juste après la conversion de latitude à radians:
# -----------------------------------------------
# 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))
pour vérifier la tendance de l'azimut dans les quatre cas que vous avez mentionnés, comparons-le à l'heure du jour:
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)
l'Image ci-jointe:
Voici une réécriture en plus idiomatique R, et plus facile à déboguer et à entretenir. C'est essentiellement la réponse de Josh, mais avec de l'Azimut calculé en utilisant les algorithmes de Josh et de Charlie pour la comparaison. J'ai également inclus les simplifications au code de date de mon autre réponse. Le principe de base est de diviser le code en beaucoup de petites fonctions que vous pouvez plus facilement écrire des tests unitaires pour.
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)
)
}
il s'agit d'une mise à jour suggérée à L'excellente réponse de Josh.
une grande partie du début de la fonction est code boilerplate pour calculer le nombre de jours depuis midi le 1er janvier 2000. Cela est beaucoup mieux traité en utilisant la date et la fonction de temps existantes de R.
je pense aussi que plutôt que d'avoir six variables différentes pour spécifier la date et l'heure, il est plus facile (et plus cohérent avec d'autres fonctions de R) de spécifier une date existante objet ou chaîne de date + chaîne de format.
voici deux fonctions d'aide
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)
}
Et le début de la fonction simplifie
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)
#...
L'autre bizarrerie est dans les lignes comme
mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360
depuis mnlong
a eu %%
appelé sur ses valeurs, ils devraient tous être non-négatif déjà, de sorte que cette ligne est superflue.
j'avais besoin de la position du soleil dans un projet Python. J'ai adapté L'algorithme de Josh O'Brien.
Merci Josh.
au cas où il pourrait être utile à n'importe qui, Voici mon adaptation.
notez que mon projet ne nécessitait que la position instantanée du soleil, donc le temps n'est pas un paramètre.
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'ai rencontré un léger problème avec un data point et les fonctions de Richie Cotton ci-dessus (dans l'implémentation du code 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
parce que dans la fonction solarAzimuthRadiansCharlie il y a eu excitation flottante autour d'un angle de 180 tel que (sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle))
est la plus petite quantité sur 1, 1.00000000000000000004440892098, qui génère un NaN que l'entrée à acos ne devrait pas être au-dessus de 1 ou en dessous de -1.
je soupçonne il pourrait y avoir des cas de bord semblables pour le calcul de Josh, où les effets d'arrondi à la virgule flottante font que l'entrée pour l'étape asin est à l'extérieur de -1:1, mais je ne les ai pas touchés dans mon ensemble de données particulier.
dans la demi-douzaine de cas que j'ai cités, le" vrai " (au milieu du jour ou de la nuit) est le moment où la question se produit de façon empirique la vraie valeur devrait être 1/-1. Pour cette raison, je serais à l'aise de corriger cela en appliquant une étape d'arrondissement dans solarAzimuthRadiansJosh
et solarAzimuthRadiansCharlie
. Je ne sais pas quelle est la précision théorique de L'algorithme NOAA (le point auquel la précision numérique cesse de compter de toute façon) mais l'arrondissement à 12 décimales a fixé les données dans mon ensemble de données.