3번째 글을 쓴지 벌써 2년이 넘게 지났습니다. 오랫동안 손을 놓고 있었지만, 올해 겨울이 지나기 전에 이 시리즈를 완결토록 하겠습니다.

4번째 글에서는 태양의 위치 계산에 관한 부분과 장동 계산에 대한 부분을 알아보겠습니다. 여기에는 모두 3개의 함수를 씁니다.

1. 태양의 위치 계산

해의 위치 계산에는 2개의 함수를 씁니다. 하나는 케플러 방정식 풀이, 하나는 해의 황경을 계산합니다.  계산 결과는 모두 황도좌표계를 기준으로 출력합니다. 

ANOMALY: 케플러 방정식 풀이
SUN: 태양의 위치 계산(황도좌표)

(1) ANOMALY(AM, EC, AT, AE)

ANOMALY 함수는 입력받은 2개의 궤도요소를 이용하여 케플러 방정식을 풀고, 그 결과를 출력하는 함수입니다. 입력받는 요소는 AM과 EC, 출력하는 요소는 AT와 AE입니다. 여기에서 AM은 평균근점이각(mean anomaly), EC는 이심률(eccentricity)이고, AT는 진근점이각(true anomaly), AE는 이심률이각(eccentric anomaly)는 입니다.

케플러 방정식은 해석적인 방법으로는 풀 수가 없으므로, 이 함수에서는 계산을 반복해서 오차를 줄여가는 식으로 진근점이각과 이심률이각을 찾아냅니다. 이 함수에서 오차는 0.000001라디안인데, 360도법으로 하면 0.21"입니다. 24절기 입기시각의 위치를 1분 이내의 정밀도로 계산하려면 태양의 위치는 2.5"미만의 오차범위 내로 계산하면 충분하다는 점을 고려한다면, 0.21"면 충분한 정밀도라고 볼 수 있습니다.

계산방법은 Astronomical Algorithms(J.Meeus)의 30장에 나오는 두 번째 방법(199쪽)을 이용하고 있습니다. 케플러 방정식에 관해서는 추후 설명토록 하겠습니다.

Private Sub ANOMALY(ByVal AM As Double, ByVal EC As Double, AT As Double, AE As Double)
  Dim A As Double, E1 As Double, E0 As Double
 
  AM = Rev(AM) * DegtoRad: E1 = AM '계산 초기값을 평균근점이각으로 대입
 
  Do  '이 부분이 케플러 방정식을 실제로 풀이하는 부분
    E0 = E1
    E1 = E0 + (AM + EC * Sin(E0) - E0) / (1 - EC * Cos(E0))
  Loop Until Abs(E1 - E0) < 0.000001
  AE = E1
 
  A = Sqr((1 + EC) / (1 - EC)) * Tan(AE / 2)
  AT = Rev(2 * Atn(A) * RadtoDeg)
  AE = Rev(AE * RadtoDeg)
End Sub

(2) SUN(JD, SL, SB, SR)

SUN 함수는 입력받은 시각의 태양의 위치를 계산하는 함수입니다. JD는 계산하고자 하는 날짜, 시각의 율리우스적일(지구시, 높은 정밀도가 필요하지 않다면 UTC로 입력할 수 있음)이고, SL은 지심황경, SB는 지심황위, SR은 전문단위(AU)로 나타낸 지심거리입니다. 계산은 궤도요소를 이용해서 케플러 방정식을 풀이하는 고전적인 방법을 사용하며, 계산 정밀도를 높이기 위하여 주요 섭동항 역시 포함되어 있습니다. 지심황위는 항상 0으로 출력합니다. 실제로 태양은 황도의 적도를 위도 방향으로 1.2" 이상 벗어나는 일이 없으므로, 이렇게 가정해도 충분한 정밀도를 보여주기 때문입니다.

Private Sub SUN(ByVal jd As Double, SL As Double, SB As Double, SR As Double)
  Dim T As Double, T2 As Double, A As Double, B As Double, L As Double, M1 As Double, EC As Double
  Dim A1 As Double, B1 As Double, C1 As Double, D1 As Double, E1 As Double, H1 As Double
  Dim D2 As Double, D3 As Double, AT As Double, AE As Double
 
  T = (jd - 2415020) / 36525: T2 = T * T
 
  A = 100.0021359 * T: B = 360 * (A - Int(A)) '궤도요소 계산
  L = 279.69668 + 0.0003025 * T2 + B
  A = 99.99736042 * T: B = 360 * (A - Int(A))
  M1 = 358.47583 - (0.00015 + 0.0000033 * T) * T2 + B
  EC = 0.01675104 - 0.0000418 * T - 0.000000126 * T2
 
  Call ANOMALY(M1, EC, AT, AE) '케플러 방정식 풀이
 
  A = 62.55209472 * T: B = 360 * (A - Int(A)) '섭동항 계산을 위한 인수 계산
  A1 = 153.23 + B
  A = 125.1041894 * T: B = 360 * (A - Int(A))
  B1 = 216.57 + B
  A = 91.56766028 * T: B = 360 * (A - Int(A))
  C1 = 312.69 + B
  A = 1236.853095 * T: B = 360 * (A - Int(A))
  D1 = 350.74 - 0.00144 * T2 + B
  E1 = 231.19 + 20.2 * T
  A = 183.1353208 * T: B = 360 * (A - Int(A))
  H1 = 353.4 + B
 
  D2 = 0.00134 * Cosd(A1) + 0.00154 * Cosd(B1) + 0.002 * Cosd(C1) '섭동(황위)
  D2 = D2 + 0.00179 * Sind(D1) + 0.00178 * Sind(E1)
 
  D3 = 0.00000543 * Sind(A1) + 0.00001575 * Sind(B1) '섭동(지심거리)
  D3 = D3 + 0.00001627 * Sind(C1) + 0.00003076 * Cosd(D1)
  D3 = D3 + 0.00000927 * Sind(H1)
 
  SL = Rev(AT + L - M1 + D2) '해의 황경 출력
  SR = 1.0000002 * (1 - EC * Cosd(AE)) + D3 '해의 지심거리 계산
  SB = 0 '황위는 0으로 가정
End Sub

함수 내부를 살펴보면, L은 태양의 평균황경, M1은 평균근점이각, EC는 이심률입니다. 섭동은 황경과 지심거리만 계산하는데, D2가 황경, D3가 지심거리에 대한 섭동입니다. 계산방법은 Astronomy with your personal computer(Peter Duffett-Smith, 제2판)에 있는 방법을 사용했습니다. 계산정밀도를 살펴보면 VSOP87과 비교했을 때 1900~2100년 사이에 황경의 평균오차는 3.46", 평균편차는 4.62", 최대오차는 20.20"이고, 거리의 최대오차는 0.000021AU입니다. 최대오차량으로 보아 입기시각 계산에 9분 정도의 오차가 발생할 수는 있지만, 실제로 이 기간동안의 음양력 계산에는 영향을 미치지 않습니다.

2. 장동

장동은 달이 지구 주위를 공전하면서 지구의 자전축이 조금씩 기울어지는 현상으로, 황경 방향으로는 ±17", 황도경사각은 ±9" 범위 내에서 영향을 줍니다. 입기시각 예측에 최대 7분 정도의 영향을 주는 요소이므로, 반드시 계산에 반영해야합니다. 장동 계산에는 단 하나의 함수만 사용하는데, 황경 방향의 요소만 계산합니다.

Nutation: 황경 방향의 장동 계산

(1) Nutation(JDE, dPsi)

Nutation 함수는 입력받은 시각의 황경방향 장동운동의 양을 계산합니다. 입력값인 JDE는 지구시로 나타낸 시각이고 출력값인 dPsi는 황경방향의 장동운동의 양입니다. 출력값의 단위는 360분법으로 나타낸 초(acsecond)입니다.

Private Sub Nutation(ByVal JDE As Double, dPsi As Double)
  Dim eps0 As Double
  Dim OM As Double, T As Double, L1 As Double, L2 As Double, T2 As Double
 
  T = (JDE - 2451545) / 36525
  T2 = T * T
  L1 = Rev(280.466457 + 36000.7698278 * T + 0.00030322 * T2 + 0.00000002 * T2 * T)
  L2 = Rev(218.3164477 + 481267.88123421 * T - 0.0015786 * T2 + T2 * T / 538841 - T2 * T2 / 65194000)
  OM = Rev(125.04452 - 1934.136261 * T + 0.0020708 * T2 + (T2 * T) / 450000) 
 
  dPsi = -17.2 * Sind(OM) - 1.32 * Sind(2 * L1) - 0.23 * Sind(2 * L2) + 0.21 * Sind(2 * OM)
End Sub

이 함수에 쓰인 계산식은 IAU1980에서 정의되어 있는 계산법에서 가장 중요한 몇 개의 항만 추려서 간단히 줄인 계산식으로 계산 정밀도는 황경방향으로 0.5" 이내입니다(입기시각 계산에서는 10초 이내의 오차요인으로 작용). 함수 내의 L1(태양의 평균황경), L2(달의 평균황경), OM(달의 승교점 경도)는 ELP-2000/82B에서 사용하는 식으로 고쳤습니다. 정밀도를 더 높이기 위한 조치이지만, 단순히 각 상수의 유효자릿수만 늘인 것으로 큰 의미는 없습니다.

다음 글에서는 해와 달의 절입시각/정삭시각 계산에 관한 내용을 다루도록 하겠습니다.