这个问题在三年多前就被问到了before。给出了答案，但是我发现解决方案存在问题。
下面的代码在R.我已经将它移植到另一种语言，但是直接在R中测试了原始代码，以确保问题不在于我的移植。
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))
}
我遇到的问题是它返回的方位角似乎是错误的。例如，如果我在12:00（夏令时）和40ºS，3ºS，3ºN和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
这些数字似乎不对。我很满意的高度 - 前两个应该大致相同，第三个触摸较低，第四个低得多。然而，第一个方位角大致应该是北方，而它给出的数字完全相反。其余三个应该大致指向南方，但只有最后一个。中间的两个位于北方，再次180º。
正如您所看到的那样，低纬度（接近赤道）也会引发一些错误
我认为故障在本节中，错误在第三行触发（从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
我用Google搜索并在C中找到了类似的代码块，转换为R，它用来计算方位角的行就像是
az <- atan(sin(ha) / (cos(ha) * sin(lat) - tan(dec) * cos(lat)))
这里的输出似乎朝着正确的方向前进，但是当它被转换回度数时，我无法让它给我正确的答案。
对代码进行修正（怀疑它只是上面几行）以使其计算正确的方位角将是非常棒的。
这似乎是一个重要的主题，所以我发布了一个比典型答案更长的时间：如果这个算法将来被其他人使用，我认为重要的是它必须附带参考文献的来源。 。
如您所知，您发布的代码无法在赤道附近或南半球的位置正常工作。
要修复它，只需在原始代码中替换这些行：
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年的一篇文章中提出的算法（Solar Energy。20：393-397）。 Walraven报道说，该方法已经成功使用了几年，以便在加利福尼亚州戴维斯（38°33'14“N，121°44'17”W）精确定位偏振辐射计。
Michalsky和Walraven的代码都包含重要/致命的错误。特别是，虽然Michalsky的算法在美国的大部分地区运作得很好，但是对于赤道附近或南半球的区域来说，它失败了（正如你所发现的那样）。 1989年，J.W。澳大利亚维多利亚州斯宾塞也注意到了同样的事情（太阳能.42（4）：353）：
亲爱的先生：
Michalsky将计算出的方位角分配给正确象限的方法（源自Walraven）在应用于南方（负）纬度时没有给出正确的值。此外，临界高程（elc）的计算将因纬度为零而失败，因为除以零。通过考虑cos（方位角）的符号，可以简单地通过将方位角分配给正确的象限来避免这两个异议。
我对您的代码的编辑是基于Spencer在已发布的评论中提出的更正。我只是稍微修改它们以确保R函数sunPosition()保持“矢量化”（即在点位置的矢量上正常工作，而不是一次需要传递一个点）。
sunPosition()为了测试sunPosition()是否正常工作，我将其结果与国家海洋和大气管理局的Solar Calculator计算结果进行了比较。在这两种情况下，太阳位置都计算在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年的一篇文章（太阳能.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.勘误表。太阳能。 43（5）：323。
斯宾塞，J.W。 1989年。关于“天文年历近似太阳位置算法（1950-2050）”的评论。太阳能。 42（4）：353。
Walraven，R。1978.计算太阳的位置。太阳能。 20：393-397。
使用上面某个链接中的“NOAA Solar Calculations”，我通过使用一种略有不同的算法改变了函数的最后部分，我希望这种算法能够无误地翻译。我已经注释掉了现在没用的代码，并在纬度到弧度转换之后添加了新算法：
# -----------------------------------------------
# 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))
为了验证你提到的四种情况中的方位角趋势，我们将其与一天中的时间进行对比：
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)
附图：
这是一个重写，它更加惯用于R，更容易调试和维护。这基本上是Josh的答案，但是使用Josh和Charlie的算法计算方位角进行比较。我还从其他答案中包含了对日期代码的简化。基本原则是将代码分成许多较小的函数，您可以更轻松地编写单元测试。
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)
)
}
这是Josh优秀答案的建议更新。
该函数的大部分开始是用于计算自2000年1月1日中午以来的天数的样板代码。使用R的现有日期和时间函数可以更好地处理。
我还认为，不是指定六个不同的变量来指定日期和时间，而是指定现有日期对象或日期字符串+格式字符串更容易（并且与其他R函数更一致）。
这是两个辅助函数
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)
}
现在，该功能的开始简化为
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)
#...
另一个奇怪的是像
mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360
由于mnlong已经让%%调用了它的值，所以它们都应该是非负的，因此这条线是多余的。
我在Python项目中需要太阳位置。我改编了Josh O'Brien的算法。
谢谢Josh。
如果它对任何人都有用，这是我的改编。
请注意，我的项目只需要即时太阳位置，因此时间不是参数。
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
我遇到了一个数据点和Richie Cotton上面的函数的一个小问题（在Charlie的代码的实现中）
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
因为在solarAzimuthRadiansCharlie函数中，在180度角附近存在浮点兴奋，使得(sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle))是超过1,10000000000000004440892098的最小量，其产生NaN作为acos的输入不应高于1或低于-1。
我怀疑Josh的计算可能存在类似的边缘情况，其中浮点舍入效应导致asin步的输入在-1：1之外，但我没有在我的特定数据集中击中它们。
在我打了六个左右的情况下，“真实”（白天或晚上）是问题发生的时候，因此根据经验，真值应为1 / -1。出于这个原因，我可以通过在solarAzimuthRadiansJosh和solarAzimuthRadiansCharlie中应用舍入步骤来解决这个问题。我不确定NOAA算法的理论精度是什么（数值精度无论如何都会停止重要），但是舍入到12位小数将数据固定在我的数据集中。