太阳的位置，一天中的时间，经度和纬度 [英] Position of the sun given time of day, latitude and longitude
问题描述
这个问题已被问及之前的三年前。有一个答案，但我发现在解决方案中有一个小故障。
下面的代码在R.我已经将它移植到另一种语言，但是已经测试过原始代码直接在R中以确保问题不在我的移植中。
sunPosition< - 函数（年，月，日，小时= 12，最小= 0，秒= 0，
lat = 46.5，长= 6.5）{
twopi < - 2 * pi
deg2rad <-pi / 180
＃获取一年中的某一天，例如2月1日= 32，3月1日= 61闰年
month.days < - c（0,31,28,31,30,31,30,31,31,30,31,30）
日< - 日+ cumsum（month.days）[月]
过期日 - - 年%% 4 == 0& （year %% 400 == 0 | year %% 100！= 0）&日期> = 60
日期[leapdays]< - day [leapdays] + 1
＃获取儒略日期 - 2400000
小时< - 小时+分钟/ 60 +秒/ 3600＃小时加分数
delta < - year-1949
leap < - trunc（delta / 4）＃前闰年
jd < - 32916.5 + delta * 365 +跃点+日+小时/ 24
＃Atronomer的账单的输入是
＃朱利安日期和JD 2451545.0（2000年1月1日中午）
时间之间的差额< - jd - 51545.
＃黄道坐标
＃平均经度
mnlong < - 280.460 + .9856474 *时间
mnlong < - mnlong %% 360
mnlong [mnlong< 0]< - mnlong [mnlong< 0] + 360
＃平均异常
mnanom < - 357.528 + .9856003 *时间
mnanom < - mnanom %% 360
mnanom [mnanom < 0]< - mnanom [mnanom< 0] + 360
mnanom < - mnanom * deg2rad
＃黄道经度和黄道倾角
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 *时间
eclong < - eclong * deg2rad
oblqec < - oblqec * deg2rad
＃天体坐标
＃赤经和赤纬
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））
＃本地坐标
＃格林威治平均恒星时间
gmst< ; - 6.697375 + .0657098242 *时间+小时
gmst< - gmst %% 24
gmst [gmst< 0]< - gmst [gmst< 0] + 24.
＃本地平均恒星时间
lmst < - gmst + long / 15.
lmst < - lmst %% 24.
lmst [lmst< 0]< - lmst [lmst< 0] + 24.
lmst < - lmst * 15. * deg2rad
＃小时角度
ha < - lmst - ra
ha [ha < -pi] - ha [ha< -pi] + twopi
ha [ha> [ha> pi - twopi
＃纬度到弧度
lat < - lat * deg2rad
＃方位和高程
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的位置运行函数：
> sunPosition（2012,12,22,12,0,0，-41,0）
$海拔
[1] 72.42113
$方位角
[1] 180.9211
> sunPosition（2012,12,22,12,0,0，-3,0）
$海拔
[1] 69.57493
$方位角
[1] -0.79713
警告信息：
asin（sin（dec）/ sin（lat））：NaNs产生
> sunPosition（2012,12,22,12,0,0,3,0）
$海拔
[1] 63.57538
$方位角
[1] - 0.6250971
警告信息：
在asin（sin（dec）/ sin（lat））中：NaNs产生
> sunPosition（2012,12,22,12,0,0,41,0）
$海拔
[1] 25.57642
$方位角
[1] 180.3084
这些数字看起来不对。我很满意的海拔 - 前两个应该大致相同，第三个降低，第四个降低。然而，第一个方位角应该大致为北，而它给出的数字是完全相反的。剩下的三个应该大致指向南方，但只有最后一个。这两个在中间位置，刚刚离开北部，再次是180度。
正如您所看到的，还有一些低纬度（赤道附近）引发的错误，
我认为错误在本节中，错误在第三行触发（从 elc ）。
＃方位角和高程
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它用来计算方位角的线将会是类似于
az < - atan（sin（ha）/（cos（ha ）* sin（lat） - tan（dec）* cos（lat）））
似乎正在朝着正确的方向前进，但我无法一直让它回到正确的答案，当它转化为度数时。
更正的代码（怀疑它只是上面几行），使它计算出正确的方位角会很棒。
就像一个重要的话题一样，所以我发布了比典型的答案更长的时间：如果这个算法将来会被其他人使用，我认为重要的是它附有对其衍生出来的文献的参考。
简短回答
正如您所注意到的，您发布的代码无法正常工作赤道线或南半球。
要解决这个问题，只需在原始代码中替换这些行即可：
<$ （sin（dec）/ sin（lat））
az [el> = elc] <-pi-az [el> = elc ]
az [el< = elc& ha> 0]< - az [el< = elc& ha> 0] + twopi
包含以下内容：
<$ （sin（az））< 0（sin（az））< 0（sin ）
az [cosAzPos& sinAzNeg]< - az [cosAzPos& sinAzNeg] + twopi
az [！cosAzPos] <-pi-az [！cosAzPos]
<它现在应该适用于全球任何地点。
讨论
你的例子中的代码几乎可以逐字从1988年的J.J. Michalsky（太阳能40：227-235）。该文章反过来改进了R. Walraven在1978年的一篇文章中提出的算法（Solar Energy 20：393-397）。 Walraven报道说，这种方法已经成功使用了好几年，可以在加利福尼亚州戴维斯（北纬38°33'14，东经121°44'17）精确定位偏振辐射计。
Michalsky's和Walraven的代码都包含重要的/致命的错误。特别是，虽然Michalsky的算法在美国大部分地区都可以正常工作，但它并没有（如您所见）靠近赤道，或在南半球。 1989年，J.W.澳大利亚维多利亚州的斯宾塞也注意到了这一点（太阳能42（4）：353）：
亲爱的主席先生：
Michalsky用于将计算得到的方位角分配给来自Walraven的正确象限的方法在应用于南方（负）纬度时不会给出正确的值。此外，临界高程（elc）的计算由于被零除而在纬度为零时将失败。通过考虑cos（方位角）的符号将方位角分配给正确的象限，可以避免这两种反对。
我的编辑到您的代码是基于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 ））
＃NOAA太阳能计算器返回的Sun的位置，
NOAA < - data.frame（elevNOAA = c（72.44,69.57,63.57,25.6），
azNOAA = c（359.09,180.79,180.62,180.3））
＃sun由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仰角方位角
＃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）&
天> = 60& ！（month == 2& day == 60）
＃oblqec< - 23.429 - 0.0000004 * time
oblqec< - 23.439 - 0.0000004 * time
更正后的版本 sunPosition（）
这是上面已验证的已更正的代码：
sunPosition< - 函数（年，月，日，小时= 12，最小= 0，秒= 0，
lat = 46.5，长= 6.5）{
twopi < - 2 * pi
deg2rad< - pi / 180
＃获取一年中的某一天，例如2月1日= 32，3月1日= 61闰年
month.days < - c（0,31,28,31,30,31,30,31,31,30,31,30）
日< - 日+ cumsum（month.days）[月]
过期日 - - 年%% 4 == 0& （year %% 400 == 0 | year %% 100！= 0）&
天> = 60& ！（month == 2& day == 60）
day [leapdays]< - day [leapdays] + 1
＃获取Julian日期 - 2400000
hour< ; - 小时+分钟/ 60 +秒/ 3600＃小时加分数
delta < - year-1949
leap < - trunc（delta / 4）＃前闰年
jd< - 32916.5 +增量* 365 +飞跃+日+小时/ 24
＃Atronomer年报的输入是
＃Julian日期和JD 2451545.0之间的差额（2000年1月1日中午）
时间< - jd - 51545.
＃黄道坐标
＃平均经度
mnlong< - 280.460 + .9856474 * time
mnlong< - mnlong％> 360
mnlong [mnlong< 0]< - mnlong [mnlong< 0] + 360
＃平均异常
mnanom < - 357.528 + .9856003 *时间
mnanom < - mnanom %% 360
mnanom [mnanom < 0]< - mnanom [mnanom< 0] + 360
mnanom < - mnanom * deg2rad
＃黄道经度和黄道倾角
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 *时间
eclong < - eclong * deg2rad
oblqec < - oblqec * deg2rad
＃天体坐标
＃赤经和赤纬
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））
＃本地坐标
＃格林威治平均恒星时间
gmst< ; - 6.697375 + .0657098242 *时间+小时
gmst< - gmst %% 24
gmst [gmst< 0]< - gmst [gmst< 0] + 24.
＃本地平均恒星时间
lmst < - gmst + long / 15.
lmst < - lmst %% 24.
lmst [lmst< 0]< - lmst [lmst< 0] + 24.
lmst < - lmst * 15. * deg2rad
＃小时角度
ha < - lmst - ra
ha [ha < -pi] - ha [ha< -pi] + twopi
ha [ha> [ha> pi - twopi
＃纬度到弧度
lat < - lat * deg2rad
＃方位和高程
el < - asin（sin （dec）* sin（lat）+ cos（dec）* cos（lat）* cos（ha））
az <-asin（-cos（dec）* sin（ha）/ cos（el））
＃逻辑和名称请参见Spencer，JW 1989年。太阳能。 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
＃} 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，J.J. 1989年。勘误。太阳能。 43（5）：323。
Spencer，J.W. 1989.关于天文年历的近似太阳位置算法（1950-2050）的评论。太阳能。 42（4）：353。
Walraven，R. 1978.计算太阳的位置。太阳能。 20：393-397。
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.
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))
}
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
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)
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
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.
The short answer
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
with these:
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.
Discussion
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).
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):
Dear Sir:
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).
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).
Accuracy of the function sunPosition()
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
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
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))
}
References:
Michalsky, J.J. 1988. The Astronomical Almanac's algorithm for approximate solar position (1950-2050). Solar Energy. 40(3):227-235.
Michalsky, J.J. 1989. Errata. Solar Energy. 43(5):323.
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. Calculating the position of the sun. Solar Energy. 20:393-397.
这篇关于太阳的位置，一天中的时间，经度和纬度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
