我想根据时间、纬度和经度计算太阳的位置。我在这里找到了这个很好的问题和答案:给定一天中的时间、纬度和经度的太阳位置。但是,当我评估函数时,我得到不正确的结果。鉴于答案的质量,我几乎可以肯定我的结果有问题,但我提出这个问题是为了记录试图解决问题。
为方便起见,下面转载了函数的代码:
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)
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)
)
}
如果我运行:
sunPosition(when = Sys.time(),lat = 43, long = -89)
结果是:
elevation azimuthJ azimuthC
1 -24.56604 55.26111 55.26111
Sys.time() 给出:
> Sys.time()
[1] "2016-09-08 09:09:05 CDT"
现在是早上 9 点,太阳已经升起。使用http://www.esrl.noaa.gov/gmd/grad/solcalc/我得到 124 的方位角和 38 的仰角,我认为这是正确的。
我想这可能是代码的问题,但我也从上面的答案中测试了 Josh 的原始 sunPosition 函数并得到了相同的结果。我的下一个想法是我的时间或时区有问题。
如上述问题中所做的那样测试冬至,仍然给出他们发现的相同结果并且看起来是正确的:
testPts <- data.frame(lat = c(-41,-3,3, 41),
long = c(0, 0, 0, 0))
time <- as.POSIXct("2012-12-22 12:00:00")
sunPosition(when = time, lat = testPts$lat, long = testPts$long)
elevation azimuthJ azimuthC
1 72.43112 359.0787 359.0787
2 69.56493 180.7965 180.7965
3 63.56539 180.6247 180.6247
4 25.56642 180.3083 180.3083
当我做同样的测试,但改变经度(-89)时,我在中午得到一个负海拔。
testPts <- data.frame(lat = c(-41,-3,3, 41),
long = c(-89, -89, -89, -89))
time <- as.POSIXct("2012-12-22 12:00:00 CDT")
sunPosition(when = time, lat = testPts$lat, long = testPts$long)
elevation azimuthJ azimuthC
1 16.060136563 107.3420 107.3420
2 2.387033692 113.3522 113.3522
3 0.001378426 113.4671 113.4671
4 -14.190786786 108.8866 108.8866