To pytanie zostało zadane przed nieco ponad trzy lata temu. Dano odpowiedź, ale znalazłem usterkę w rozwiązaniu.

Poniższy kod jest w R. Przenieśliłem go na inny język, jednak przetestowałem oryginalny kod bezpośrednio w R, aby upewnić się, że problem nie dotyczył mojego przenoszenia.

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

Problem polega na tym, że azymut, który zwraca, wydaje się zły. Na przykład, jeśli uruchomię funkcję podczas przesilenia letniego (południowego) o godzinie 12:00 dla lokalizacji 0ºE i 41ºS, 3ºS, 3ºN i 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

Te liczby po prostu nie wydają się prawidłowe. Wysokość, z której jestem zadowolony - pierwsze dwie powinny być mniej więcej takie same, trzecia odrobinę niższa, a czwarta znacznie niższa. Jednak pierwszy azymut powinien być mniej więcej na północ, podczas gdy liczba, którą podaje, jest całkowitym przeciwieństwem. Pozostałe trzy powinny wskazywać mniej więcej dokładnie na południe, ale tylko ostatnia tak robi. Dwie w środkowym punkcie tuż przy północy, znowu 180º na zewnątrz.

Jak widać, istnieje również kilka błędów wywoływanych przy niskich szerokościach geograficznych (blisko równika)

Uważam, że usterka jest w tej sekcji, a błąd jest wyzwalany w trzeciej linii (zaczynając od 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

Przeszukałem go i znalazłem podobny fragment kodu w C, przekonwertowany na R, linia, której używa do obliczenia azymutu, wyglądałaby tak

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

Wydaje się, że dane wyjściowe zmierzają we właściwym kierunku, ale po prostu nie mogę ich uzyskać, aby dać mi właściwą odpowiedź przez cały czas, gdy jest konwertowany z powrotem na stopnie.

Korekta kodu (podejrzewam, że to tylko kilka linii powyżej), aby obliczyć poprawny azymut, byłaby fantastyczna.

Wydaje się, że to ważny temat, więc zamieściłem dłuższą niż typowa odpowiedź: jeśli ten algorytm ma być używany przez innych w przyszłości, myślę, że ważne jest, aby towarzyszyły mu odniesienia do literatury, z której został wyprowadzony .

Krótka odpowiedź

Jak zauważyłeś, Twój opublikowany kod nie działa poprawnie w przypadku lokalizacji w pobliżu równika lub na półkuli południowej.

Aby to naprawić, po prostu zamień te wiersze w oryginalnym kodzie:

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

z nimi:

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

Powinien teraz działać w dowolnym miejscu na świecie.

Dyskusja

Kod w twoim przykładzie jest prawie dosłownie zaadaptowany z artykułu JJ Michalsky'ego z 1988 roku (Solar Energy. 40: 227-235). Ten artykuł z kolei udoskonalił algorytm przedstawiony w artykule R. Walravena z 1978 roku (Solar Energy. 20: 393-397). Walraven poinformował, że metoda była z powodzeniem stosowana od kilku lat do precyzyjnego pozycjonowania radiometru polaryzacyjnego w Davis w Kalifornii (38 ° 33 '14 "N, 121 ° 44' 17" W).

Kod Michalsky'ego i Walravena zawiera ważne / krytyczne błędy. W szczególności, chociaż algorytm Michalsky'ego działa dobrze w większości Stanów Zjednoczonych, zawodzi (jak odkryłeś) w obszarach w pobliżu równika lub na półkuli południowej. W 1989 roku JW Spencer z Wiktorii w Australii zauważył to samo (Solar Energy. 42 (4): 353):

Szanowny Panie:

Metoda Michalsky'ego przypisywania obliczonego azymutu do właściwego kwadrantu, pochodząca z Walraven, nie daje poprawnych wartości, gdy jest stosowana dla południowych (ujemnych) szerokości geograficznych. Ponadto obliczenie krytycznej wysokości (elc) nie powiedzie się dla zerowej szerokości geograficznej z powodu dzielenia przez zero. Obu tych zastrzeżeń można uniknąć, po prostu przypisując azymut do właściwego kwadrantu, biorąc pod uwagę znak cos (azymut).

Moje zmiany w Twoim kodzie są oparte na poprawkach sugerowanych przez Spencera w opublikowanym Komentarzu. Po prostu nieco je zmieniłem, aby upewnić się, że funkcja R sunPosition()pozostaje „wektoryzowana” (tj. Działa poprawnie na wektorach lokalizacji punktów, a nie musi być przekazywana po jednym punkcie na raz).

Dokładność funkcji sunPosition()

Aby sprawdzić, czy to sunPosition()działa poprawnie, porównałem jego wyniki z wynikami obliczonymi przez kalkulator słoneczny National Oceanic and Atmospheric Administration . W obu przypadkach pozycje słońca obliczono dla południa (12:00) południowego przesilenia letniego (22 grudnia) 2012. Wszystkie wyniki były zgodne z dokładnością do 0,02 stopnia.

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

Inne błędy w kodzie

W przesłanym kodzie są co najmniej dwa inne (dość drobne) błędy. Pierwsza z nich powoduje, że 29 lutego i 1 marca lat przestępnych są liczone jako 61 dzień roku. Drugi błąd wynika z literówki w oryginalnym artykule, którą poprawił Michalsky w notatce z 1989 r. (Solar Energy. 43 (5): 323).

Ten blok kodu pokazuje obraźliwe wiersze, zakomentowane i zaraz po nich poprawione wersje:

# 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

Poprawiona wersja sunPosition()

Oto poprawiony kod, który został zweryfikowany powyżej:

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

Bibliografia:

Michalsky, JJ 1988. Algorytm Astronomical Almanac's dla przybliżonej pozycji Słońca (1950-2050). Energia słoneczna. 40 (3): 227-235.

Michalsky, JJ 1989. Errata. Energia słoneczna. 43 (5): 323.

Spencer, JW 1989. Komentarze na temat "Algorytm Astronomical Almanac's for Approximate Solar Position (1950-2050)". Energia słoneczna. 42 (4): 353.

Walraven, R. 1978. Obliczanie położenia słońca. Energia słoneczna. 20: 393-397.

Używając „NOAA Solar Calculations” z jednego z linków powyżej, zmieniłem nieco ostatnią część funkcji, używając nieco innego algorytmu, który, mam nadzieję, przetłumaczył bez błędów. Skomentowałem teraz bezużyteczny kod i dodałem nowy algorytm zaraz po konwersji szerokości geograficznej na radiany:

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

Aby zweryfikować trend azymutu we wspomnianych czterech przypadkach, wykreślmy go względem pory dnia:

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)

Obraz załączony:

Oto przepis, który jest bardziej idiomatyczny do R i łatwiejszy do debugowania i utrzymania. Zasadniczo jest to odpowiedź Josha, ale z azymutem obliczonym przy użyciu algorytmów Josha i Charliego dla porównania. Dodałem również uproszczenia do kodu daty z mojej drugiej odpowiedzi. Podstawową zasadą było podzielenie kodu na wiele mniejszych funkcji, dla których można łatwiej napisać testy jednostkowe.

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

To sugerowana aktualizacja doskonałej odpowiedzi Josha.

Znaczna część początku funkcji jest kodem standardowym do obliczania liczby dni od południa 1 stycznia 2000 r. O wiele lepiej radzi sobie z tym przy użyciu istniejącej funkcji daty i czasu R.

Myślę też, że zamiast mieć sześć różnych zmiennych do określenia daty i godziny, łatwiej (i bardziej spójnie z innymi funkcjami R) jest określić istniejący obiekt daty lub ciągi daty + ciągi formatujące.

Oto dwie funkcje pomocnicze

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

A początek funkcji upraszcza się teraz do

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

Inną osobliwością są takie linie

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

Ponieważ mnlongnie miał %%wezwał swoich wartości, wszystkie one powinny być już nieujemna, więc ta linia jest zbędna.

Potrzebowałem pozycji słońca w projekcie Pythona. Dostosowałem algorytm Josha O'Briena.

Dziękuję Josh.

Na wypadek, gdyby przydała się komukolwiek, oto moja adaptacja.

Zauważ, że mój projekt wymagał tylko natychmiastowej pozycji słońca, więc czas nie jest parametrem.

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

Napotkałem niewielki problem z punktem danych i powyższymi funkcjami Richiego Cottona (w implementacji kodu Charliego)

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

ponieważ w funkcji solarAzimuthRadiansCharlie podekscytowanie zmiennoprzecinkowe występowało wokół kąta 180, (sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle))czyli najmniejszej wartości powyżej 1, 1,0000000000000004440892098, co generuje NaN, ponieważ wartość wejściowa acos nie powinna przekraczać 1 ani poniżej -1.

Podejrzewam, że mogą istnieć podobne przypadki krawędzi w obliczeniach Josha, w których efekty zaokrąglania zmiennoprzecinkowego powodują, że dane wejściowe dla kroku asin są poza -1: 1, ale nie trafiłem ich w moim konkretnym zbiorze danych.

W kilku lub kilkunastu przypadkach, w których trafiłem, „prawda” (w środku dnia lub nocy) występuje wtedy, gdy problem występuje, więc z empirycznego punktu widzenia prawdziwa wartość powinna wynosić 1 / -1. Z tego powodu czułbym się komfortowo, stosując krok zaokrąglania wewnątrz solarAzimuthRadiansJoshi solarAzimuthRadiansCharlie. Nie jestem pewien, jaka jest teoretyczna dokładność algorytmu NOAA (punkt, w którym dokładność numeryczna i tak przestaje mieć znaczenie), ale zaokrąglenie do 12 miejsc po przecinku naprawiło dane w moim zestawie danych.