주어진 시간, 위도 및 경도가 지정된 태양의 위치
이 질문은 요청 된 전 이상 3 년 전에 조금. 주어진 답변이 있었지만 솔루션에서 결함을 발견했습니다.
아래 코드는 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))
}
내가 치는 문제는 그것이 반환하는 방위각이 잘못 보인다는 것입니다. 예를 들어, 0ºE 및 41ºS, 3ºS, 3ºN 및 41ºN 위치에 대해 12:00에 (남부)하지에서 함수를 실행하면 :
> 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
이 숫자는 옳지 않은 것 같습니다. 내가 만족하는 고도-처음 두 개는 거의 같고 세 번째는 터치가 더 낮아야하며 네 번째는 훨씬 더 낮아야합니다. 그러나 첫 번째 방위각은 대략 정북이어야하며, 제공하는 숫자는 완전히 반대입니다. 나머지 3 개는 대략 남쪽을 가리켜 야하지만 마지막 하나만 그렇습니다. 중간에있는 두 개는 북쪽에서 약간 떨어져 다시 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
나는 주위를 훑어 보았고 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]
이제 지구상의 모든 위치에서 작동합니다.
토론
이 예제의 코드는 JJ Michalsky (Solar Energy. 40 : 227-235)의 1988 년 기사에서 거의 그대로 적용되었습니다. 그 기사는 R. Walraven (Solar Energy. 20 : 393-397)의 1978 년 기사에 제시된 알고리즘을 차례로 개선했습니다. Walraven은이 방법이 캘리포니아 데이비스 (38 ° 33 '14 "N, 121 ° 44'17"W)에 편광 복사계를 정확하게 배치하는 데 몇 년 동안 성공적으로 사용되었다고보고했습니다.
Michalsky와 Walraven의 코드에는 모두 중요 / 치명적인 오류가 포함되어 있습니다. 특히 Michalsky의 알고리즘은 대부분의 미국에서 잘 작동하지만 적도 근처 또는 남반구 지역에서는 실패합니다. 1989 년에 호주 빅토리아에있는 JW Spencer는 똑같은 사실을 언급했습니다 (Solar Energy. 42 (4) : 353) :
근계:
계산 된 방위각을 Walraven에서 파생 된 올바른 사분면에 할당하는 Michalsky의 방법은 남부 (음수) 위도에 적용될 때 올바른 값을 제공하지 않습니다. 또한 0으로 나누기 때문에 위도가 0 인 경우 임계 고도 (elc) 계산이 실패합니다. cos (azimuth)의 부호를 고려하여 방위각을 올바른 사분면에 할당함으로써이 두 가지 반대를 피할 수 있습니다.
코드에 대한 수정 사항은 게시 된 주석에서 Spencer가 제안한 수정 사항을 기반으로합니다. R 함수 sunPosition()가 '벡터화'상태를 유지 하도록 약간 변경했습니다 (즉, 한 번에 한 지점을 전달하는 대신 지점 위치의 벡터에서 제대로 작동 함).
기능의 정확성 sunPosition()
sunPosition()제대로 작동하는지 테스트하기 위해 그 결과를 국립 해양 대기 청의 태양열 계산기 에서 계산 한 결과와 비교했습니다 . 두 경우 모두 2012 년 12 월 22 일 남부하지 (12 월 22 일)의 정오 (오후 12시) 동안 태양 위치를 계산했습니다. 모든 결과는 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))
}
참조 :
Michalsky, JJ 1988. 대략적인 태양 위치 (1950-2050)에 대한 천문 연감의 알고리즘. 태양 에너지. 40 (3) : 227-235.
Michalsky, JJ 1989. 정오표. 태양 에너지. 43 (5) : 323.
Spencer, JW 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의 기존 날짜 및 시간 함수를 사용하는 것이 훨씬 더 좋습니다.
또한 날짜와 시간을 지정하기 위해 6 개의 다른 변수를 사용하는 것보다 기존 날짜 개체 또는 날짜 문자열 + 형식 문자열을 지정하는 것이 더 쉽고 (다른 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했다 %%그것의 값이라는 라인이 불필요하므로, 그들은 모두, 이미 음수가되어야한다.
파이썬 프로젝트에서 태양 위치가 필요했습니다. 저는 Josh O'Brien의 알고리즘을 적용했습니다.
고마워요 조쉬.
누구에게나 유용 할 수 있다면 여기에 내 적응이 있습니다.
내 프로젝트에는 즉각적인 태양 위치 만 필요하므로 시간은 매개 변수가 아닙니다.
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의 기능에 약간의 문제가 발생했습니다 (찰리의 코드 구현에서).
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, 1.0000000000000004440892098보다 작은 양이어서 acos에 대한 입력이 1보다 크거나 -1보다 작아서는 안되므로 NaN을 생성합니다.
Josh의 계산과 유사한 경우가있을 수 있습니다. 부동 소수점 반올림 효과로 인해 asin 단계에 대한 입력이 -1 : 1을 벗어나게되지만 특정 데이터 세트에서 해당 항목에 도달하지 않았습니다.
6 개 정도의 경우에이 문제가 발생했을 때 "참"(낮이나 밤의 중간)은 문제가 발생하는 경우이므로 경험적으로 실제 값은 1 / -1이어야합니다. 따라서 solarAzimuthRadiansJosh및 내에서 반올림 단계를 적용하여 문제를 해결할 수 solarAzimuthRadiansCharlie있습니다. NOAA 알고리즘의 이론적 정확도가 무엇인지 (수치 정확도가 문제가되지 않는 지점) 확실하지 않지만 소수점 이하 12 자리로 반올림하면 데이터 세트의 데이터가 수정되었습니다.
'Programing' 카테고리의 다른 글
| asp.net-mvc의 캐싱 (0) | 2020.09.24 |
|---|---|
| 목록에서 id 삭제 (0) | 2020.09.24 |
| 별칭 함수를 입력하고 캐스팅하지 않고 사용할 수있는 이유는 무엇입니까? (0) | 2020.09.24 |
| 안전하게 홈 디렉토리 변경 (0) | 2020.09.24 |
| iOS 4.0 iTunes 백업에서 Manifest.mbdb 파일을 구문 분석하는 방법 (0) | 2020.09.24 |