这似乎是一个重要的主题，所以我发布了一个比典型答案更长的答案：如果这个算法将来要被其他人使用，我认为重要的是附上对其衍生文献的引用。
简短的回答
正如您所指出的，您发布的代码对于赤道附近或南半球的位置无法正常工作。
要修复它，只需替换原始代码中的以下行：
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]
它现在应该适用于全球任何地点。
讨论
您示例中的代码几乎是逐字改编自 J.J. 1988 年的文章。 Michalsky（太阳能。40：227-235）。该文章进而改进了 R. Walraven 1978 年文章中提出的算法（太阳能。20：393-397）。 Walraven 报告称，该方法已成功使用多年，在加利福尼亚州戴维斯（北纬 38° 33' 14"，西经 121° 44' 17"）精确定位偏振辐射计。
Michalsky 和 Walraven 的代码都包含重要/致命的错误。特别是，虽然米哈斯基的算法在美国大部分地区运行良好，但对于赤道附近或南半球的地区却失败了（正如您所发现的）。 1989 年，J.W.澳大利亚维多利亚州的 Spencer 也指出了同样的事情（Solar Energy. 42(4):353）：
尊敬的先生：
Michalsky 将计算的方位角分配给正确象限的方法（源自 Walraven）在应用于南纬（负）纬度时不会给出正确的值。此外，由于除以零，对于零纬度，临界高程 (elc) 的计算将失败。只需考虑 cos(方位角) 的符号，将方位角分配给正确的象限，即可避免这两种反对意见。
我对您的代码的编辑是基于斯宾塞在该发表的评论中建议的更正。我只是对它们进行了一些修改以确保 R 函数sunPosition()保持“矢量化”（即在点位置的矢量上正常工作，而不需要一次传递一个点）。
功能的准确性sunPosition()
为了测试这一点sunPosition()工作正常，我将其结果与国家海洋和大气管理局的计算结果进行了比较太阳能计算器 http://www.esrl.noaa.gov/gmd/grad/solcalc/。在这两种情况下，太阳位置都是在 2012 年南方夏至（12 月 22 日）中午（下午 12:00）计算的。所有结果都一致在 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))
}
参考：
米哈尔斯基，J.J. 1988。天文年历的太阳近似位置算法（1950-2050）。太阳能。 40（3）：227-235。
米哈尔斯基，J.J. 1989。勘误表。太阳能。 43（5）：323。
斯宾塞，J.W. 1989年。对《天文年鉴近似太阳位置的算法（1950-2050）》的评论。太阳能。 42（4）：353。
Walraven, R. 1978。计算太阳的位置。太阳能。 20：393-397。
