拝啓。

計算された方位を正しい象限に割り当てるMichalskyの方法は、Walravenから派生したものですが、南（負）緯度に適用すると正しい値を与えません。また、臨界高度(elc)の計算も緯度がゼロの場合、ゼロで割ってしまうので失敗します。これらの反論は両方とも、cos(方位)の符号を考慮することによって正しい象限に方位を割り当てることで簡単に回避することができます。

あなたのコードへの私の編集は、その公開されたコメントでSpencerによって提案された修正に基づきます。私は単にそれらを多少変更し、R 関数 sunPosition() が「ベクトル化」されたままであること（つまり、一度に1つの点を渡される必要があるのではなく、点の位置のベクトル上で適切に動作すること）を保証するために、それらを多少変更しただけです。

関数の精度 sunPosition()

をテストするために sunPosition() が正しく動作することを確認するために、その結果を米国海洋大気庁（National Oceanic and Atmospheric Administration）の ソーラー計算機 . いずれも、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 つの他の (非常に小さな) エラーがあります。1 つは、うるう年の 2 月 29 日と 3 月 1 日が両方ともその年の 61 日目として集計されることです。2つ目の誤りは、原著論文のタイプミスに由来するもので、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, J.J. 1988. 天文年鑑の太陽位置の近似値アルゴリズム（1950-2050）。太陽エネルギー。40(3):227-235.

Michalsky, J.J. 1989. Errata. 太陽エネルギー。43(5):323.

スペンサー、J.W.1989。The Astronomical Almanac's Algorithm for Approximate Solar Position (1950-2050) "へのコメント。太陽エネルギー。42(4):353.

Walraven, R. 1978. 太陽の位置の計算. 太陽エネルギー。20:393-397.