Gündüz güneş, enlem ve boylam, zaman verilen konum

Bu soru before bir küçük üç yılı aşkın bir süre önce soruldu. Cevap verildi, ancak çözümde bir hata buldum.

Aşağıdaki kodu başka bir dil için taşıdık ettik R. olduğunu, ancak doğrudan orijinal kodu R sorunu benim taşıma ile değildi emin olmak için test ettik.

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

Basıyorum sorun verir Azimut yanlış görünüyor. Eğer kaçarsam örneğin, (Güney) yaz işlevi yerleri 0ºE ve 41ºS, 3ºS, 3ºN ve 41ºN için 12:00 de gündönümü:

> 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

Bu rakamlar doğru değil. Mutluyum yükseklik - ilk ikisi aynı, üçüncüsü bir dokunuş daha düşük kabaca ve çok daha düşük dördüncü olmalıdır. Ancak ilk Azimut verir numarasını tam tersi ise kabaca Kuzeye olmalıdır. Kalan üç kabaca Güneye çevir, ancak sadece sonuncusu. Ortadaki iki Kuzey, 180 paylaştığınız.

Gördüğünüz gibi, hataları da düşük enlemlerde (ekvatora yakın) ile tetiklenen bir çift vardır

Hata hata üçüncü satırı (elc ile başlayan) tetiklenmesi ile bu bölümde, inanıyorum.

  # 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

Etrafta araştırdım ve C, dönüştürülmüş R hesaplamak için kullandığı çizgi Azimut için kod benzer bir yığın gibi bir şey olurdu bulundu

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

Çıkış burada doğru yönde ilerliyor gibi görünüyor, ama ben sadece bana doğru cevabı derece geri dönüştürülmüş olduğunda her zaman vermeye ikna edemez.

Kod (sadece birkaç satır yukarıda şüpheli) bir düzeltme doğru Azimut hesaplamak için harika olurdu.

CEVAP

6 Ocak 2012, Cuma

Bu gibi önemli bir konu, o yüzden yayınladım bir daha tipik bir cevap: Eğer bu algoritma için kullanılan, Diğerleri gelecek, bence önemli olan bu olmalı, beraberinde referanslar için literatürden olan oldu türetilmiş.

Kısa cevap

Belirtildiği gibi, yayınlanan kodunuz doğru, güney yarım kürede ekvatora yakın yerlerde, ya da çalışmıyor.

Bunu düzeltmek için, sadece orijinal kodu: bu satırları değiştir

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

bunlar:

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

Şimdi dünya üzerindeki herhangi bir yer için çalışması gerekir.

Tartışma

Sizin örnek kod neredeyse birebir J. J. Michalsky tarafından 1988 bir madde (Güneş Enerjisi. ' dan uyarlanmıştır 40:227-235). Bu makalede, bir algoritma R. Walraven tarafından 1978 bir makalede sunulan (Güneş Enerjisi. rafine 20:393-397). Walraven yöntemi başarıyla birkaç yıl için tam davis'teki polarize gazları düzeyleri pozisyon için kullanıldığına dair rapor, (38° 33' 14" N, 121° 44' 17" W). CA

Hem Michalsky ve Walraven kodu önemli/önemli hatalar içeriyor.Michalsky algoritma Amerika Birleşik Devletleri en iyi çalışır iken özellikle, buldum gibi ekvatorun, ya da güney yarımkürede yakın alanlar için başarısız olur. 1989, J. W. Victoria, Avustralya, dikkat aynı şey (Güneş Enerjisi. Spencer 42(4):353):

Sevgili Efendim:

Doğru çeyreği için hesaplanan Azimut atama için Michalsky yöntem, Walraven türetilen, (olumsuz) Güney enlemleri için uygulandığında doğru değerler vermez. Daha kritik yükseklik hesaplama (elı) sıfıra bölme yüzünden sıfır bir enlem için başarısız olur. Hem bu itirazlar sadece cos(Azimut) işaretini dikkate alarak doğru çeyreğine Azimut atayarak önlenebilir.

Kodunuzu benim düzenlemeler düzeltmeler yayınlanan Yorum Spencer tarafından önerilen dayanmaktadır. Ben sadece onları biraz değişmiş sağlamak için R fonksiyonu sunPosition() kalır 'vectorized' (yani çalışan düzgün vektörler nokta konumları, yerine ihtiyacı iletilecek bir noktada bir zaman).

Fonksiyonun doğruluk 17**

sunPosition() düzgün çalışıp çalışmadığını test etmek için, bu Ulusal Okyanus tarafından hesaplanan sonuçlar karşılaştırıldığında oldum ve Atmosfer İdaresi Solar Calculator. Her iki durumda da, güneş pozisyonları Güney yaz gündönümü 22 Aralık 2012 tarihinde öğle (12:00) için hesaplandı. Tüm sonuçları içinde 0.02 derece anlaşmaya varıldı.

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

Kodu diğer hatalar

En azından ilan kodu diğer iki (çok) küçük hatalar vardır. İlk sıçrama yıl 29 Şubat ve 1 Mart yılın gün 61 olarak sayımı için iki neden. İkinci hata 1989 bir not (Güneş Enerjisi. içinde Michalsky tarafından düzeltilmiştir orijinal yazıda bir yazım hatası, türemiştir 43(5):323).

Bu kod bloğu soruna neden olan satırları gösterir, diye ve hemen düzeltilmiş sürümleri izledi:

# 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

`sunPosition()` düzeltilmiş versiyonu

İşte yukarıda doğrulandı düzeltilmiş kod:

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