Diese Frage wurde vor etwas mehr als drei Jahren gestellt. Es wurde eine Antwort gegeben, aber ich habe einen Fehler in der Lösung gefunden.

Der folgende Code ist in R. Ich habe ihn in eine andere Sprache portiert, jedoch den Originalcode direkt in R getestet, um sicherzustellen, dass das Problem nicht mit meiner Portierung zusammenhängt.

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

Das Problem, auf das ich stoße, ist, dass der zurückgegebene Azimut falsch erscheint. Wenn ich zum Beispiel die Funktion an der (südlichen) Sommersonnenwende um 12:00 Uhr für die Standorte 0ºE und 41ºS, 3ºS, 3ºN und 41ºN ausführe:

> 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

Diese Zahlen scheinen einfach nicht richtig zu sein. Die Höhe, mit der ich zufrieden bin - die ersten beiden sollten ungefähr gleich sein, die dritte eine Berührung tiefer und die vierte viel niedriger. Der erste Azimut sollte jedoch ungefähr genau nach Norden gehen, während die Zahl, die er angibt, genau das Gegenteil ist. Die restlichen drei sollten ungefähr nach Süden zeigen, aber nur der letzte. Die beiden in der Mitte gleich nördlich, wieder 180º raus.

Wie Sie sehen können, gibt es auch einige Fehler, die bei den niedrigen Breiten ausgelöst werden (schließen Sie den Äquator).

Ich glaube, der Fehler liegt in diesem Abschnitt, wobei der Fehler in der dritten Zeile (beginnend mit elc) ausgelöst wird .

  # 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

Ich googelte herum und fand einen ähnlichen Codeabschnitt in C, der in R konvertiert wurde. Die Zeile, die zur Berechnung des Azimuts verwendet wird, wäre ungefähr so

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

Die Ausgabe hier scheint in die richtige Richtung zu gehen, aber ich kann sie einfach nicht dazu bringen, mir immer die richtige Antwort zu geben, wenn sie wieder in Grad umgewandelt wird.

Eine Korrektur des Codes (vermutlich nur die wenigen Zeilen oben), um den korrekten Azimut zu berechnen, wäre fantastisch.

Dies scheint ein wichtiges Thema zu sein, daher habe ich eine länger als typische Antwort veröffentlicht: Wenn dieser Algorithmus in Zukunft von anderen verwendet werden soll, ist es meiner Meinung nach wichtig, dass er von Verweisen auf die Literatur begleitet wird, aus der er abgeleitet wurde .

Die kurze Antwort

Wie Sie bereits bemerkt haben, funktioniert Ihr Postleitzahl an Orten in der Nähe des Äquators oder auf der südlichen Hemisphäre nicht ordnungsgemäß.

Um dies zu beheben, ersetzen Sie einfach diese Zeilen in Ihrem Originalcode:

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

mit diesen:

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

Es sollte jetzt für jeden Ort auf der Welt funktionieren.

Diskussion

Der Code in Ihrem Beispiel stammt fast wörtlich aus einem Artikel von JJ Michalsky aus dem Jahr 1988 (Solar Energy. 40: 227-235). Dieser Artikel verfeinerte wiederum einen Algorithmus, der 1978 in einem Artikel von R. Walraven (Solar Energy. 20: 393-397) vorgestellt wurde. Walraven berichtete, dass die Methode seit mehreren Jahren erfolgreich eingesetzt wurde, um ein polarisierendes Radiometer in Davis, CA (38 ° 33 '14 "N, 121 ° 44 '17" W) präzise zu positionieren.

Sowohl Michalskys als auch Walravens Code enthält wichtige / schwerwiegende Fehler.Während Michalskys Algorithmus in den meisten USA einwandfrei funktioniert, schlägt er (wie Sie festgestellt haben) in Gebieten in der Nähe des Äquators oder auf der südlichen Hemisphäre fehl. 1989 bemerkte JW Spencer aus Victoria, Australien, dasselbe (Solar Energy. 42 (4): 353):

Lieber Herr:

Michalskys Methode zur Zuordnung des berechneten Azimuts zum richtigen Quadranten, abgeleitet von Walraven, liefert keine korrekten Werte, wenn sie für südliche (negative) Breiten angewendet wird. Ferner schlägt die Berechnung der kritischen Höhe (elc) für einen Breitengrad von Null aufgrund der Division durch Null fehl. Diese beiden Einwände können einfach vermieden werden, indem der Azimut dem richtigen Quadranten zugewiesen wird, indem das Vorzeichen von cos (Azimut) berücksichtigt wird.

Meine Änderungen an Ihrem Code basieren auf den Korrekturen, die Spencer in diesem veröffentlichten Kommentar vorgeschlagen hat. Ich habe sie einfach etwas geändert, um sicherzustellen, dass das R funktioniertsunPosition() "vektorisiert" bleibt (dh ordnungsgemäß an Vektoren von Punktpositionen arbeitet, anstatt Punkt für Punkt übergeben zu werden).

Genauigkeit der Funktion sunPosition()

Um zu testen, ob dies sunPosition()korrekt funktioniert, habe ich die Ergebnisse mit denen verglichen, die mit dem Solarrechner der National Oceanic and Atmospheric Administration berechnet wurden . In beiden Fällen wurden die Sonnenpositionen für den Mittag (12:00 Uhr) an der südlichen Sommersonnenwende (22. Dezember) 2012 berechnet. Alle Ergebnisse stimmten mit 0,02 Grad überein.

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

Andere Fehler im Code

Der veröffentlichte Code enthält mindestens zwei weitere (recht geringfügige) Fehler. Der erste bewirkt, dass der 29. Februar und der 1. März der Schaltjahre beide als Tag 61 des Jahres gewertet werden. Der zweite Fehler ergibt sich aus einem Tippfehler im Originalartikel, der von Michalsky in einer Notiz von 1989 korrigiert wurde (Solar Energy. 43 (5): 323).

Dieser Codeblock zeigt die fehlerhaften Zeilen an, auskommentiert und sofort gefolgt von korrigierten Versionen:

# 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

Korrigierte Version von sunPosition()

Hier ist der korrigierte Code, der oben überprüft wurde:

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

Verweise:

Michalsky, JJ 1988. Der Algorithmus des astronomischen Almanachs für die ungefähre Sonnenposition (1950-2050). Solarenergie. 40 (3): 227 & ndash; 235.

Michalsky, JJ 1989. Errata. Solarenergie. 43 (5): 323.

Spencer, JW 1989. Kommentare zu "Der Algorithmus des astronomischen Almanachs für die ungefähre Sonnenposition (1950-2050)". Solarenergie. 42 (4): 353.

Walraven, R. 1978. Berechnung des Sonnenstandes. Solarenergie. 20: 393 & ndash; 397.

Mit "NOAA Solar Calculations" von einem der obigen Links habe ich den letzten Teil der Funktion ein wenig geändert, indem ich einen etwas anderen Algorithmus verwendet habe, der hoffentlich fehlerfrei übersetzt wurde. Ich habe den jetzt unbrauchbaren Code auskommentiert und den neuen Algorithmus direkt nach der Umrechnung von Breitengrad zu Bogenmaß hinzugefügt:

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

Um den Azimuttrend in den vier genannten Fällen zu überprüfen, zeichnen wir ihn gegen die Tageszeit:

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)

Bild im Anhang:

Hier ist eine Neufassung, die für R idiomatischer und einfacher zu debuggen und zu warten ist. Es ist im Wesentlichen Joshs Antwort, aber mit Azimut, der sowohl mit Josh als auch mit Charlies Algorithmen zum Vergleich berechnet wird. Ich habe auch die Vereinfachungen des Datumscodes aus meiner anderen Antwort aufgenommen. Das Grundprinzip bestand darin, den Code in viele kleinere Funktionen aufzuteilen, für die Sie einfacher Unit-Tests schreiben können.

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

Dies ist ein empfohlenes Update für Joshs ausgezeichnete Antwort.

Ein Großteil des Starts der Funktion besteht aus einem Boilerplate-Code zur Berechnung der Anzahl der Tage seit dem Mittag am 1. Januar 2000. Dies wird viel besser mit der vorhandenen Datums- und Uhrzeitfunktion von R erledigt.

Ich denke auch, dass es einfacher (und konsistenter mit anderen R-Funktionen) ist, ein vorhandenes Datumsobjekt oder eine Datumszeichenfolge + Formatzeichenfolgen anzugeben, anstatt sechs verschiedene Variablen zum Angeben von Datum und Uhrzeit zu haben.

Hier sind zwei Hilfsfunktionen

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

Und der Start der Funktion vereinfacht sich jetzt zu

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

Die andere Kuriosität ist in den Zeilen wie

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

Da mnlonghat sich hatte %%auf seine Werte genannt, sollten sie alle nicht-negativen schon sein, so dass diese Linie überflüssig ist.

Ich brauchte Sonnenstand in einem Python-Projekt. Ich habe den Algorithmus von Josh O'Brien angepasst.

Danke Josh.

Für den Fall, dass es für irgendjemanden nützlich sein könnte, hier ist meine Anpassung.

Beachten Sie, dass mein Projekt nur eine sofortige Sonnenposition benötigte, sodass die Zeit kein Parameter ist.

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

Ich bin auf ein kleines Problem mit einem Datenpunkt und den Funktionen von Richie Cotton oben gestoßen (bei der Implementierung von Charlies Code).

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

weil in der solarAzimuthRadiansCharlie-Funktion eine Gleitkommaerregung um einen Winkel von 180 (sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle))aufgetreten ist, so dass dies der kleinste Betrag über 1 ist, 1.0000000000000004440892098, der ein NaN erzeugt, da die Eingabe in acos nicht über 1 oder unter -1 liegen sollte.

Ich vermute, dass es ähnliche Randfälle für Joshs Berechnung gibt, bei denen Gleitkomma-Rundungseffekte dazu führen, dass die Eingabe für den Asin-Schritt außerhalb von -1: 1 liegt, aber ich habe sie in meinem speziellen Datensatz nicht getroffen.

In den etwa einem halben Dutzend Fällen, in denen ich dies getroffen habe, ist das "wahre" (mitten am Tag oder in der Nacht), wenn das Problem so empirisch auftritt, dass der wahre Wert 1 / -1 sein sollte. Aus diesem Grund würde ich das gerne beheben, indem ich einen Rundungsschritt innerhalb von solarAzimuthRadiansJoshund anwende solarAzimuthRadiansCharlie. Ich bin mir nicht sicher, wie hoch die theoretische Genauigkeit des NOAA-Algorithmus ist (der Punkt, an dem die numerische Genauigkeit ohnehin keine Rolle mehr spielt), aber das Runden auf 12 Dezimalstellen hat die Daten in meinem Datensatz korrigiert.