이것은 중요한 주제처럼 보이므로 일반적인 답변보다 더 긴 답변을 게시했습니다.이 알고리즘이 미래에 다른 사람들이 사용할 예정이라면 파생 된 문헌에 대한 참조를 첨부하는 것이 중요하다고 생각합니다. .

짧은 대답

앞서 언급했듯이 게시 된 코드는 적도 근처 또는 남반구의 위치에서 제대로 작동하지 않습니다.

이 문제를 해결하려면 원래 코드에서 다음 줄을 바꾸면됩니다.

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

이들과 함께 :

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

이제 지구상의 모든 위치에서 작동합니다.

토론

이 예제의 코드는 JJ Michalsky (Solar Energy. 40 : 227-235)의 1988 년 기사에서 거의 그대로 적용되었습니다. 이 기사는 R. Walraven (Solar Energy. 20 : 393-397)의 1978 년 기사에 제시된 알고리즘을 차례로 개선했습니다. Walraven은이 방법이 캘리포니아 데이비스 (38 ° 33 ’14 “N, 121 ° 44’17″W)에 편광 복사계를 정확하게 배치하는 데 몇 년 동안 성공적으로 사용되었다고보고했습니다.

Michalsky와 Walraven의 코드에는 모두 중요 / 치명적인 오류가 포함되어 있습니다. 특히, Michalsky의 알고리즘은 대부분의 미국에서 잘 작동하지만 적도 근처 또는 남반구 지역에서는 실패합니다. 1989 년에 호주 빅토리아에있는 JW Spencer는 똑같은 사실을 언급했습니다 (Solar Energy. 42 (4) : 353) :

근계:

계산 된 방위각을 Walraven에서 파생 된 올바른 사분면에 할당하는 Michalsky의 방법은 남부 (음수) 위도에 적용될 때 올바른 값을 제공하지 않습니다. 또한 0으로 나누기 때문에 위도가 0 인 경우 임계 고도 (elc) 계산이 실패합니다. cos (azimuth)의 부호를 고려하여 방위각을 올바른 사분면에 할당함으로써이 두 가지 반대를 피할 수 있습니다.

코드에 대한 편집 내용은 게시 된 주석에서 Spencer가 제안한 수정 사항을 기반으로합니다. R 함수 sunPosition()가 ‘벡터화’상태를 유지 하도록 약간 변경했습니다 (즉, 한 번에 한 지점을 전달하는 대신 지점 위치의 벡터에서 제대로 작동 함).

기능의 정확성 sunPosition()

sunPosition()제대로 작동하는지 테스트하기 위해 그 결과를 국립 해양 대기 청의 태양열 계산기 에서 계산 한 결과와 비교했습니다 . 두 경우 모두 2012 년 12 월 22 일 남하의 정오 (오후 12시)에 태양 위치를 계산했습니다. 모든 결과는 0.02도 이내로 일치했습니다.

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

코드의 기타 오류

게시 된 코드에 적어도 두 개의 다른 (매우 사소한) 오류가 있습니다. 첫 번째는 윤년의 2 월 29 일과 3 월 1 일이 모두 해당 연도의 61 일로 집계되도록합니다. 두 번째 오류는 원래 기사의 오타에서 비롯된 것으로 Michalsky가 1989 년 메모에서 수정했습니다 (Solar Energy. 43 (5) : 323).

이 코드 블록은 문제가되는 줄을 표시하고 주석 처리 한 후 바로 수정 된 버전을 표시합니다.

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

위에서 확인한 수정 된 코드는 다음과 같습니다.

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

참조 :

Michalsky, JJ 1988. 대략적인 태양 위치에 대한 천문 연감의 알고리즘 (1950-2050). 태양 에너지. 40 (3) : 227-235.

Michalsky, JJ 1989. 정오표. 태양 에너지. 43 (5) : 323.

Spencer, JW 1989. “근사적인 태양 위치에 대한 천문 연감의 알고리즘 (1950-2050)”에 대한 주석. 태양 에너지. 42 (4) : 353.

Walraven, R. 1978 년. 태양의 위치 계산. 태양 에너지. 20 : 393-397.