이번 편에서는 해의 절입시각 계산, 달의 정삭시각 계산에 사용하는 함수에 대해 알아보겠습니다. 

입기시각은 중기 또는 절기에 이르는 순간의 시점, 정삭시각은 달이 삭에 이르는 순간의 시점입니다. 계산정확도가 1분 이내라면 대부분 정확하게 음양력을 계산할 수 있습니다. 다만 극히 드물게 입기시각이나 정삭시각이 자정에 극히 가까울 때라면 음력 초하루 날짜의 계산이나 치윤에 문제가 될 수 있지만, 1900~2100년 사이에는 이런 사례가 없으므로 안심해도 됩니다.   

입기시각과 정삭시각 계산은 음력 계산 과정에서 가장 중요한 과정 가운데 하나입니다. 두 시각을 얼마나 정확하게 계산하느냐에 따라 음양력 변환의 정확도가 결정되기 때문입니다. 음력-양력 상호 변환 함수 모듈에서는 정삭시각과 입기시각 계산을 위해 4개의 함수를 사용합니다. 

1. 달의 정삭시각 계산

달의 정삭시각 계산은 단 하나의 함수를 이용하여 처리합니다.

NewMoon: 달의 정삭 시각 계산

(1) NewMoon(JDN)

정삭시각 계산은 기본적으로 해와 달의 황경 차이를 이용해서 계산합니다. 해와 달의 겉보기 황경의 차이가 0도가 시점이 정삭시각이 되는데요, 달의 위치를 높은 정밀도로 구하려면 계산식이 많이 길어지고 복잡해지므로, 여기에서는 그 방법 대신 직접 달의 정삭시각을 구하는 방식을 택했습니다. 

이 함수에서 사용하고 있는 계산법은 J.Meeus의 『Astronomical Algorithms』(Second Edition), 제 49장 Phases of the Moon에서 설명하고 있는 방법을 사용했습니다. 이 계산식은 VSOP87과 ELP-2000에 기반을 두고 만들었으며, 음양력 변환에는 충분한 정도의 정확도를 보여즙니다.

입력인 JDN 변수에는 계산하고자하는 정삭일 근처의 날짜(율리우스 적일)를 입력해주면 되며, 출력은 율리우스적일로 나타낸 지구시(TT)로 됩니다.

Private Function NewMoon(ByVal JDN As Double) As Double

  Dim jd As Double, cor As Double, M As Double, Mp As Double, F As Double, T As Double, OMG As Double

  Dim YR As Double, k As Double, nph As Double, dt As Double, T2 As Double, T3 As Double, T4 As Double

  Dim A1 As Double, A2 As Double, A3 As Double, A4 As Double, A5 As Double, A6 As Double, A7 As Double

  Dim A8 As Double, A9 As Double, A10 As Double, E As Double, E2 As Double, TX As Double

  Dim A11 As Double, A12 As Double, A13 As Double, A14 As Double, cor2 As Double

  YR = InvJDYear(JDN)

  k = JDN - JULIANDAY(YR, 1, 1)

  YR = YR + k / 365.25

  nph = (YR - 2000) * 12.3685

  nph = Int(nph)

  k = nph: T = k / 1236.85

  T2 = T * T: T3 = T2 * T: T4 = T2 * T2

  jd = 2451550.09766 + 29.530588861 * k + 0.00015437 * T2 - 0.00000015 * T3 + 0.00033 * 0.00000000073 * T4

       

  TX = (jd - 2451545) / 36525

  E = 1 - 0.002516 * TX - 0.0000074 * TX * TX: E2 = E * E

  M = Rev(2.5534 + 29.1053567 * k - 0.0000014 * T2 - 0.00000011 * T3)

  Mp = Rev(201.5643 + 385.81693528 * k + 0.0107582 * T2 + 0.00001238 * T3 - 0.000000058 * T4)

  F = Rev(160.7108 + 390.67050284 * k - 0.0016118 * T2 - 0.00000227 * T3 + 0.000000011 * T4)

  OMG = Rev(124.7746 - 1.56375588 * k + 0.0020672 * T2 + 0.00000215 * T3)

  

  A1 = 299.77 + 0.107408 * k - 0.009173 * T2

  A2 = 251.88 + 0.016321 * k

  A3 = 251.83 + 26.651886 * k

  A4 = 349.42 + 36.412478 * k

  A5 = 84.66 + 18.206239 * k

  A6 = 141.74 + 53.303771 * k

  A7 = 207.14 + 2.453732 * k

  A8 = 154.84 + 7.30686 * k

  A9 = 34.52 + 27.261239 * k

  A10 = 207.19 + 0.121824 * k

  A11 = 291.34 + 1.844379 * k

  A12 = 161.72 + 24.198154 * k

  A13 = 239.56 + 25.513099 * k

  A14 = 331.55 + 3.592518 * k

  

  cor = -0.4072 * Sind(Mp) + 0.17241 * E * Sind(M) + 0.01608 * Sind(2 * Mp) _

        + 0.01039 * Sind(2 * F) + 0.00739 * E * Sind(Mp - M) - 0.00514 * E * Sind(Mp + M) _

        + 0.00208 * E2 * Sind(2 * M) - 0.00111 * Sind(Mp - 2 * F) - 0.00057 * Sind(Mp + 2 * F) _

        + 0.00056 * E * Sind(2 * Mp + M) - 0.00042 * Sind(3 * Mp) + 0.00042 * E * Sind(M + 2 * F) _

        + 0.00038 * E * Sind(M - 2 * F) - 0.00024 * E * Sind(2 * Mp - M) - 0.00017 * Sind(OMG) _

        - 0.00007 * Sind(Mp + 2 * M) + 0.00004 * Sind(2 * Mp - 2 * F) + 0.00004 * Sind(3 * M) _

        + 0.00003 * Sind(Mp + M - 2 * F) + 0.00003 * Sind(2 * Mp + 2 * F) - 0.00003 * Sind(Mp + M + 2 * F) _

        + 0.00003 * Sind(Mp - M + 2 * F) - 0.00002 * Sind(Mp - M - 2 * F) - 0.00002 * Sind(3 * Mp + M) + 0.00002 * Sind(4 * Mp)

  

  cor2 = 0.000325 * Sind(A1) + 0.000165 * Sind(A2) + 0.000164 * Sind(A3) + 0.000126 * Sind(A4) _

       + 0.00011 * Sind(A5) + 0.000062 * Sind(A6) + 0.00006 * Sind(A7) + 0.000056 * Sind(A8) _

       + 0.000047 * Sind(A9) + 0.000042 * Sind(A10) + 0.00004 * Sind(A11) + 0.000037 * Sind(A12) _

       + 0.000035 * Sind(A13) + 0.000023 * Sind(A14)

      

  jd = jd + cor + cor2

  dt = JDtoTDT(jd) - jd

  NewMoon = jd - dt + timezone / 24

End Function

해와 달의 위치를 이용해서 정삭시각을 구하려면 해와 달은 모두 지심좌표를 취하고, 해는 광행차, 달은 광차(light-time correction)를 적용하여 계산하면 됩니다. 장동은 해와 달 모두에 똑같이 적용되는 양이므로 따로 계산에 포함시키지 않아도 됩니다.

2. 해의 입기시각 계산

해의 입기시각 계산은 4편에서 다룬 해의 위치계산함수를 이용하여 대략 해의 위치를 구한 다음, 점차 오차를 줄여나가는 방식으로 계산합니다. 계산에는 모두 3개의 함수를 씁니다. 계산순서는 Set24Julgi, cJunggi, CalcJulGi 입니다(정확히는 CalcJulGi에서 cJunggi를 보조함수로 사용).

CalcJulGi: 12중기 입기 시각 저장
cJunggi: 12중기 입기 시각 계산
Set24Julgi: 12중기 계산에 필요한 기초 자료 입력

또한 절입시각을 저장하기 위해 Junggi라는 변수배열을 사용합니다. 이 변수는 Julgi12라는 구조체를 사용하는데, 각각의 요소마다 4개의 자료가 저장됩니다. 각각의 요소에 대한 내용은 아래 코의의 주석을 참고하기 바랍니다. 계산 결과는 RealDay라는 부분에 저장하는데, 계산 연도 전년의 소설(小雪)부터 계산 연도 다음 해의 우수(雨水)까지 입기시각을 율리우스적일 형태로 저장합니다. timezone 변수에는 시간 단위로 나타낸 시간권을 기억하는 변수로, 시간권이 UTC보다 빠르면 양의 값(+), 느리면 음의 값(-)을 가집니다. 한국은 +9입니다.

Type Julgi12

  MonNumber As Byte '중기에 해당되는 월 번호. 예를 들어 동지는 11.

  Ref_Day As Double '양력 1월 1일부터 해당 중기까지의 대략적인 일수. 춘분은 80.

  Longitude As Double '해당 중기의 태양 황경. 동지는 270도

  RealDay As Double '해당 중기의 입기시각(계산결과 저장)

End Type

Dim Junggi(15) As Julgi12, timezone As Double

여기에 쓰인 계산법은 비효율적인 부분이 있으며, 조만간 개선된 알고리즘으로 대체할 예정입니다(현재 다른 곳에는 개선된 방법을 쓰고 있음). 

(1) Set24Julgi()

이 함수는 Junggi 변수배열에 계산을 위한 기초자료를 저장줍니다. 

Private Sub Set24Julgi()

  Junggi(0).MonNumber = 10: Junggi(0).Longitude = 240: Junggi(0).Ref_Day = -43 

  Junggi(1).MonNumber = 11: Junggi(1).Longitude = 270: Junggi(1).Ref_Day = -13     

  Junggi(2).MonNumber = 12: Junggi(2).Longitude = 300: Junggi(2).Ref_Day = 20    

  Junggi(3).MonNumber = 1: Junggi(3).Longitude = 330: Junggi(3).Ref_Day = 50     

  Junggi(4).MonNumber = 2: Junggi(4).Longitude = 0: Junggi(4).Ref_Day = 80     

  Junggi(5).MonNumber = 3: Junggi(5).Longitude = 30: Junggi(5).Ref_Day = 110     

  Junggi(6).MonNumber = 4: Junggi(6).Longitude = 60: Junggi(6).Ref_Day = 140     

  Junggi(7).MonNumber = 5: Junggi(7).Longitude = 90: Junggi(7).Ref_Day = 170     

  Junggi(8).MonNumber = 6: Junggi(8).Longitude = 120: Junggi(8).Ref_Day = 200      

  Junggi(9).MonNumber = 7: Junggi(9).Longitude = 150: Junggi(9).Ref_Day = 230      

  Junggi(10).MonNumber = 8: Junggi(10).Longitude = 180: Junggi(10).Ref_Day = 260    

  Junggi(11).MonNumber = 9: Junggi(11).Longitude = 210: Junggi(11).Ref_Day = 290     

  Junggi(12).MonNumber = 10: Junggi(12).Longitude = 240: Junggi(12).Ref_Day = 320      

  Junggi(13).MonNumber = 11: Junggi(13).Longitude = 270: Junggi(13).Ref_Day = 350      

  Junggi(14).MonNumber = 12: Junggi(14).Longitude = 300: Junggi(14).Ref_Day = 385     

  Junggi(15).MonNumber = 1: Junggi(15).Longitude = 330: Junggi(15).Ref_Day = 415     

End Sub

Ref_Day는 양력 계산연도 1월 1일로부터의 대략적인 날짜 수입니다. 입기시각을 찾는 참고자료로 cJunggi함수에서 이용합니다. 추후 사용할 개선된 함수에서는 불필요한 부분입니다.

(2)  cJunggi(cYear, LonSun, RefDay)

이 함수는 Set24Julgi에서 cYear, LonSun, RefDay에서 입력된 값을 바탕으로 입기시각을 계산해서 출력합니다. 출력결과는 율리우스적일 형태이며, 시간권을 적용한 지방시(LST)입니다. 시간권은 timezone 변수에 저장된 값을 이용합니다.  

입력받는 인수는 3가지인데, cYear는 계산하려는 연도, LonSun는 계산하려는 중기의 태양황경, RefDay는 1월 1일부터 계산하려는 중기까지의 대략적인 날짜수입니다. 초기값으로 1월 1일의 위치를 입력한 다음, 태양의 평균운동속도를 이용하여 점점 오차를 줄여나가는 식으로 작동합니다. 

입기시각의 계산에는 해의 정확한 겉보기 지심황경을 구하기 위해 장동과 광행차 보정을 포함하고 있으며, 계산오차는 시간으로 1초 이내로 하고 있지만, 태양의 위치계산 함수의 오차로 인해 실제의 입기시각과 통상 2~3분, 최대 10분까지의 차이를 보일 수 있습니다. 

Private Function cJunggi(ByVal cYear As Double, ByVal LonSun As Double, ByVal RefDay As Double) As Double

  Dim JDyear As Double, aDay As Double

  Dim LamSun As Double, dt As Double, dLam As Double, tJD As Double, i As Long

  Dim dl As Double, SL As Double, SB As Double, SR As Double, JDE As Double

  

  dt = 0: i = 0

  JDyear = JULIANDAY(cYear, 1, 0) - 0.5 '연초의 율리우스 적일 계산

  

  tJD = JDyear + RefDay

  JDE = JDtoTDT(tJD)

  SUN JDE, SL, SB, SR

  LamSun = Rev(SL - 0.005691611 / SR)

  dLam = AngDistLon(LamSun, LonSun)

  

  Do

    dt = dLam / Oneday '보정항은 해의 평균 운동속도를 이용하여 추정

    

    If LonSun > 357 Or LonSun < 3 Then

      If LamSun > 180 Then LamSun = LamSun - 360

    End If

    

    If LonSun > LamSun Then

      tJD = tJD + dt

    Else

      tJD = tJD - dt

    End If

    

    JDE = JDtoTDT(tJD)

    SUN JDE, SL, SB, SR

    Nutation JDE, dl '장동 계산

    LamSun = Rev(SL + dl / 3600 - 0.005691611 / SR) ''(- 0.005691611 / SR)이 광행차 보정 부분

    

    dLam = AngDistLon(LamSun, LonSun) '해와 달의 경도 차이

    i = i + 1

  Loop Until (dLam / Oneday * 86400) < 1 Or i > 50  '계산오차가 1초 이내 또는 반복횟수가 50회 미만인 동안 계산 반복

  cJunggi = tJD + timezone / 24

End Function

이 함수도 조만간 더 나은 방법으로 바꿀 예정입니다. 

(3) CalcJulGi(JDt)

CalcJulGi 함수는 입력받은 율리우스적일이 포함된 해(양력)의 입기시각을 한 번에 계산하여 결과를 Junggi 변수배열에 저장하는 역할을 합니다. 입력값인 JDt는 계산하고자하는 연도가 포함된 날의 율리우스적일입니다. 

Private Sub CalcJulGi(ByVal JDt As Double)

  Dim i As Integer, nYear As Double

  

  nYear = InvJDYear(JDt) '입력받은 날의 양력 연도 얻기

  

  For i = 0 To 15 'cJunggi 함수를 이용해 절입시간을 한 번에 계산한 다음, 입기일의 0시로 바꾸어 

Junggi 변수에 저장.

    Junggi(i).RealDay = GetJD0(cJunggi(nYear, Junggi(i).Longitude, Junggi(i).Ref_Day)) + 0.5

  Next i

End Sub

이번 글에서 설명한 함수들은 조만간 논리적 구조가 더 명료한 함수로 대체할 예정입니다.

다음번 글에서는 음양력 변환에 사용하는 함수에 대한 설명을 하겠습니다. 아마도 2편 정도 더 쓰면 마무리가 될 것 같습니다.