给定时间、纬度和经度的太阳位置 [英] Position of the sun given time of day, latitude and longitude
问题描述
此问题已被问到 之前 三年前.给出了一个答案，但是我发现解决方案中存在一个小故障.
This question has been asked before a little over three years ago. There was an answer given, however I've found a glitch in the solution.
下面的代码是用 R 编写的.我已经将它移植到了另一种语言，但是已经直接在 R 中测试了原始代码以确保问题不在于我的移植.
Code below is in R. I've ported it to another language, however have tested the original code directly in R to ensure the issue wasn't with my porting.
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 对位置 0ºE 和 41ºS、3ºS、3ºN 和 41ºN 运行函数:
The problem I'm hitting is that the azimuth it returns seems wrong. For example, if I run the function on the (southern) summer solstice at 12:00 for locations 0ºE and 41ºS, 3ºS, 3ºN and 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 度.
These numbers just don't seem right. The elevation I'm happy with - the first two should be roughly the same, the third a touch lower, and the fourth much lower. However the first azimuth should be roughly due North, whereas the number it gives is the complete opposite. The remaining three should point roughly due South, however only the last one does. The two in the middle point just off North, again 180º out.
正如你所看到的，低纬度(靠近赤道)也触发了一些错误
As you can see there are also a couple of errors triggered with the low latitudes (close the equator)
我认为问题出在这一部分，错误在第三行(以 elc 开头)触发.
I believe the fault is in this section, with the error being triggered at the third line (starting with 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
我用谷歌搜索并在 C 中找到了类似的代码块，转换为 R 用于计算方位角的行将类似于
I googled around and found a similar chunk of code in C, converted to R the line it uses to calculate the azimuth would be something like
az <- atan(sin(ha) / (cos(ha) * sin(lat) - tan(dec) * cos(lat)))
这里的输出似乎正朝着正确的方向发展，但是当它转换回度数时，我无法一直给我正确的答案.
The output here seems to be heading in the right direction, but I just can't get it to give me the right answer all the time when it's converted back to degrees.
更正代码(怀疑只是上面的几行)以使其计算正确的方位角会很棒.
A correction of the code (suspect it's just the few lines above) to make it calculate the correct azimuth would be fantastic.
推荐答案
这似乎是一个重要的话题，所以我发布了一个比典型答案更长的答案:如果这个算法将来要被其他人使用，我想重要的是它伴随着对它的来源文献的引用.
This seems like an important topic, so I've posted a longer than typical answer: if this algorithm is to be used by others in the future, I think it's important that it be accompanied by references to the literature from which it has been derived.
正如您所指出的，您发布的代码在赤道附近或南半球的位置无法正常工作.
As you've noted, your posted code does not work properly for locations near the equator, or in the southern hemisphere.
要修复它，只需替换原始代码中的这些行:
To fix it, simply replace these lines in your original code:
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]
它现在应该适用于地球上的任何位置.
It should now work for any location on the globe.
您示例中的代码几乎逐字改编自 J.J. 1988 年的一篇文章.Michalsky(太阳能.40:227-235).那篇文章又改进了 R. Walraven 1978 年的一篇文章中提出的算法(Solar Energy. 20:393-397).Walraven 报告说，该方法已成功用于精确定位位于加利福尼亚州戴维斯的偏振辐射计(38° 33' 14" N，121° 44' 17" W).
The code in your example is adapted almost verbatim from a 1988 article by J.J. Michalsky (Solar Energy. 40:227-235). That article in turn refined an algorithm presented in a 1978 article by R. Walraven (Solar Energy. 20:393-397). Walraven reported that the method had been used successfully for several years to precisely position a polarizing radiometer in Davis, CA (38° 33' 14" N, 121° 44' 17" W).
Michalsky 和 Walraven 的代码都包含重要/致命的错误.特别是，虽然 Michalsky 的算法在美国大部分地区运行良好，但它在靠近赤道或南半球.1989 年，J.W.澳大利亚维多利亚州的 Spencer 也注意到了同样的事情(太阳能.42(4):353):
Both Michalsky's and Walraven's code contains important/fatal errors. In particular, while Michalsky's algorithm works just fine in most of the United States, it fails (as you've found) for areas near the equator, or in the southern hemisphere. In 1989, J.W. Spencer of Victoria, Australia, noted the same thing (Solar Energy. 42(4):353):
尊敬的先生:
Michalsky 用于将计算的方位角分配给正确的象限的方法(源自 Walraven)在应用于南纬(负)纬度时没有给出正确的值.此外，由于被零除，临界高程 (elc) 的计算将失败为零纬度.通过考虑 cos(azimuth) 的符号将方位角分配给正确的象限，可以避免这两种反对意见.
Michalsky's method for assigning the calculated azimuth to the correct quadrant, derived from Walraven, does not give correct values when applied for Southern (negative) latitudes. Further the calculation of the critical elevation (elc) will fail for a latitude of zero because of division by zero. Both these objections can be avoided simply by assigning the azimuth to the correct quadrant by considering the sign of cos(azimuth).
我对您的代码所做的编辑基于 Spencer 在该已发布评论中建议的更正.我只是稍微改变了它们以确保 R 函数 sunPosition() 保持矢量化"(即在点位置的矢量上正常工作，而不是一次传递一个点).
My edits to your code are based on the corrections suggested by Spencer in that published Comment. I have simply altered them somewhat to ensure that the R function sunPosition() remains 'vectorized' (i.e. working properly on vectors of point locations, rather than needing to be passed one point at a time).
为了测试 sunPosition() 是否正常工作，我将其结果与美国国家海洋和大气管理局的 太阳能计算器.在这两种情况下，太阳位置都是在 2012 年南夏至(12 月 22 日)正午(12:00 PM)计算的.所有结果都在 0.02 度以内.
To test that sunPosition() works correctly, I've compared its results with those calculated by the National Oceanic and Atmospheric Administration's Solar Calculator. In both cases, sun positions were calculated for midday (12:00 PM) on the southern summer solstice (December 22nd), 2012. All results were in agreement to within 0.02 degrees.
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).
Other errors in the code
There are at least two other (quite minor) errors in the posted code. The first causes February 29th and March 1st of leap years to both be tallied as day 61 of the year. The second error derives from a typo in the original article, which was corrected by Michalsky in a 1989 note (Solar Energy. 43(5):323).
此代码块显示了有问题的行，已注释掉并紧随其后的是更正的版本:
This code block shows the offending lines, commented out and followed immediately by corrected versions:
# 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()
这是上面验证过的更正代码:
Corrected version of sunPosition()
Here is the corrected code that was verified above:
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. The Astronomical Almanac 的近似太阳位置算法(1950-2050).太阳能.40(3):227-235.
References:
Michalsky, J.J. 1988. The Astronomical Almanac's algorithm for approximate solar position (1950-2050). Solar Energy. 40(3):227-235.
迈克尔斯基，J.J.1989. 勘误表.太阳能.43(5):323.
Michalsky, J.J. 1989. Errata. Solar Energy. 43(5):323.
斯宾塞，J.W.1989. 对近似太阳位置的天文年历算法(1950-2050)"的评论.太阳能.42(4):353.
Spencer, J.W. 1989. Comments on "The Astronomical Almanac's Algorithm for Approximate Solar Position (1950-2050)". Solar Energy. 42(4):353.
Walraven, R. 1978.计算太阳的位置.太阳能.20:393-397.
Walraven, R. 1978. Calculating the position of the sun. Solar Energy. 20:393-397.
这篇关于给定时间、纬度和经度的太阳位置的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
