이번에는 지난 번에 이야기한 대로 시간 계산에 쓰이는 함수를 알아보겠습니다. 그러나 그 전에 이 모듈에 쓰이는 상수에 관해 설명할 필요가 있는 것 같습니다. 원래 지난 번 글에서 다루어야 했지만 깜빡했네요.

1. 상수

이 모듈은 모두 8개의 상수를 사용합니다. 이 수들은 모듈의 첫 부분에 모여 있습니다.

Const pi As Double = 3.14159265358979 '원주율
Const hpi As Double = 1.5707963267949 '원주율의 반
Const RadtoDeg As Double = 180 / pi '라디안 단위를 360분법으로 바꿀 때 곱하는 수
Const DegtoRad As Double = pi / 180 '360분법으로 나타낸 각도를 라디안 단위로 바꿀 때 곱하는 수
Const MoonMonth As Double = 29.5305882 '달의 공전 주기
Const MoonDay As Double = 12.190749387105 '달이 하루에 이동하는 경도 방향의 평균 각거리
Const OneYear As Double = 365.24219 '1태양년의 길이
Const Oneday As Double = 360 / OneYear '해가 하루 동안 이동하는 경도 방향의 평균 각거리

-pi: 원주율입니다.

-hpi: 원주율을 반으로 나눈 값입니다.

-RadtoDeg, DegtoRad: 라디안 단위에서 360분법으로 또는그 반대로 각도의 단위를 바꾸어줄 때 곱하는 수입니다. 라디안을 360분법으로 바꾸려면 호도법으로 나타낸 각도에 RadtoDeg을, 360분법으로 나타낸 각도를 라디안 단위로 바꾸려면 360분법으로 나타낸 각도에 DegtoRad을 곱해주면 됩니다.

-MoonMonth: 1평균 삭망월입니다. 원래 합삭일을 계산할 때 쓰던 상수인데, 계산 방법을 바꾸면서 더 이상 사용하고 있지는 않습니다.

-MoonDay: 360을 1 평균 삭망월로 나눈 값입니다. 달이 경도 방향으로 하루에 얼마나 이동하는지를 대략 나타내고 있습니다. 역시 이 모듈에서 쓰고 있지는 않습니다.

-OneYear: 일 단위로 나타낸 1 평균 태양년(회귀년)입니다.

-Oneday: 해가 경도 방향으로 하루에 얼마나 이동하는지를 나타낸 값입니다. 360을 1 회귀년으로 나눈 값입니다.

2. 시간 계산에 쓰이는 함수

시간 계산에 쓰이는 함수는 모두 다섯 개가 있습니다. 

JULIANDAY: 입력한 양력 날짜의 율리우스 적일 계산
InvJD: 율리우스 적일에서 양력 날짜 찾기
InvJDYear: 율리우스 적일에서 양력 연도 찾기
GetJD0: 입력한 율리우스 적일의 시간을 0시로 바꾸어 주는 함수
JDtoTDT: 율리우스 적일의 기준 시각을 UT(세계시)에서 TT(지구시)로 바꾸어 주는 함수

(1) JULIANDAY(Year, Month, Day)

JULIANDAY 함수는 상용력으로 표시한 날짜를 율리우스 적일로 바꾸어 주는 함수입니다. 율리우스 적일은 기원전 4713년 1월 1일 12시부터 흘러간 날짜 수로 표시하는 날짜 체계입니다. 이 방법은 일상적으로 흔히 쓰는 방법은 아니지만, 천문 계산에서는 널리 쓰이고 있습니다. 연도와 월, 일이 단절된 형태인 일반 달력과는 다르게 특정 시점을 기준으로 연속된 하나의 수로 특정한 시점을 나타낼 수 있으므로 천문 계산에는 아주 편리한 날짜 표기법입니다.

여기에서 쓰고 있는 함수는 Jean Meeus가 쓴 『Astronomical Algirithms』이라는 책에서 소개하고 있는 방법(제2판 60~62쪽)으로 율리우스 적일을 계산합니다. 입력하는 날짜(상용력으로 표시한 날짜)는 1582년 10월 4일까지는 율리우스력을, 그 이후는 그레고리력을 사용합니다.

Private Function JULIANDAY(ByVal Year As Double, ByVal Month As Double, ByVal Day As Double) As Double 
    Dim ggg As Double, S As Double, A As Double, J1 As Double, tJD As Double
   
    ggg = 1
    If Year < 1582 Then ggg = 0
    If Year <= 1582 And Month < 10 Then ggg = 0
    If Year <= 1582 And Month = 10 And Day < 5 Then ggg = 0
    tJD = -1 * Int(7 * (Int((Month + 9) / 12) + Year) / 4)
    S = 1
    If (Month - 9) < 0 Then S = -1
    A = Abs(Month - 9)
    J1 = Int(Year + S * Int(A / 7))
    J1 = -1 * Int((Int(J1 / 100) + 1) * 3 / 4)
    tJD = tJD + Int(275 * Month / 9) + Day + (ggg * J1)
    tJD = tJD + 1721027 + 2 * ggg + 367 * Year
    tJD = tJD
    JULIANDAY = tJD
End Function

율리우스 적일 체계에서 하루 미만의 시간은 소수점을 써서 나타내지만, 이 모듈에서는 필요하지 않아서 날짜 단위 까지만 계산하도록 하였습니다.

(2) InvJD(julday, Year, Month, Day)

율리우스 적일로 나타낸 날짜를 상용력으로 바꾸어 주는 함수입니다. 역시 『Astronomical Algirithms』에서 소개하고 있는 방법(63~64쪽)에 따라 계산합니다. 1582년 10월 4일(율리우스 적일로는 2299160.0)까지는 율리우스력으로, 그 이후는 그레고리력으로 출력합니다. 역시 날짜 단위까지만 계산하고 소수점 부분은 생략합니다.

Private Sub InvJD(ByVal julday As Double, Year As Double, Month As Double, Day As Double)
    Dim Z As Double, F As Double, A As Double, i As Double, B As Double
    Dim c As Double, D As Double, T As Double, rj As Double, JJ As Double, RH As Double
    Dim mMe As Double, aae As Double
   
    Z = Int(julday + 0.5)
    F = julday + 0.5 - Z
    If Z < 2299161 Then
       A = Z
    Else
       i = Int((Z - 1867216.25) / 36524.25)
       A = Z + 1 + i - Int(i / 4)
    End If
   
    B = A + 1524
    c = Int((B - 122.1) / 365.25)
    D = Int(365.25 * c)
    T = Int((B - D) / 30.6)
    rj = B - D - Int(30.6001 * T) + F
    JJ = Int(rj)
   
    If T < 14 Then
       mMe = T - 1
    Else
       If T = 14 Or T = 15 Then mMe = T - 13
    End If
   
    If mMe > 2 Then
       aae = c - 4716
    ElseIf mMe = 1 Or mMe = 2 Then
       aae = c - 4715
    End If
   
    Year = aae: Month = mMe: Day = JJ
End Sub

(3) InvJDYear(julday)

InvJD 함수와 동일한 역할을 하지만 연도만 출력해 줍니다.

Private Function InvJDYear(ByVal julday As Double) As Double
    Dim Z As Double, A As Double, i As Double, B As Double
    Dim c As Double, D As Double, T As Double, mMe As Double
   
    Z = Int(julday + 0.5)
   
    If Z < 2299161 Then
       A = Z
    Else
       i = Int((Z - 1867216.25) / 36524.25)
       A = Z + 1 + i - Int(i / 4)
    End If
   
    B = A + 1524
    c = Int((B - 122.1) / 365.25)
    D = Int(365.25 * c)
    T = Int((B - D) / 30.6)
   
    If T < 14 Then
       mMe = T - 1
    Else
       If T = 14 Or T = 15 Then mMe = T - 13
    End If
   
    If mMe > 2 Then
       InvJDYear = c - 4716
    ElseIf mMe = 1 Or mMe = 2 Then
       InvJDYear = c - 4715
    End If
End Function

이 함수는 아래처럼 간결하게 바꾸어 줄 수 있습니다.

Private Function InvJDYear(ByVal julday As Double) As Double
    Dim A As Double, B As Double, C As Double
   

    InvJD julday, A, B, C  'InvJD 함수를 이용해 계산(A: 연도, B: 월, C: 일)
    InvJDYear = A '위 값에서 연도만 출력
End Function

(4) GetJD0(DateJD)

율리우스 적일은 기원전 4713년 1월 1일 정오(12시)를 기점으로 삼고 있습니다. 따라서 자정(00시)을 기준으로 시간을 표시하면 xxx.5 형태가 됩니다. 이 함수는 소수점이 포함된 율리우스 적일을 입력 받아 그 날의 자정(00시)으로 바꾸어 주는 함수입니다. 예를 들어 2425780.25라는 값을 입력 받으면 해당 일의 자정 때인 2425779.5로 바꾸어 출력해 줍니다. 정수 부분이 바뀌었다는 점에 주의를 해야 하는데 이유는 이렇습니다. 2425780.25일은 율리우스 적일 체계의 원기를 기준으로 2425780.25일이 지난 시각입니다. 0.25일은 6시간이고, 따라서 2425780.25일은 2425780번째 날의 정오(12시)에서 6시간 후가 됩니다. 그래서 2425780.25일은 2425780일의 18시(12+6=18)가 되며 이날의 자정은 2425780.25에서 18시간(0.75일)을 뺀 2425779.5가 됩니다.

Private Function GetJD0(ByVal DateJD As Double) As Double
  If DateJD - Int(DateJD) >= 0.5 Then '만일 정오 이후라면 소수점 부분을 없앤 다음 0.5를 더해 주고
    GetJD0 = Int(DateJD) + 0.5
  Else '정오가 되기 전이라면 소수점 부분을 없앤 다음 0.5를 빼 줌.
    GetJD0 = Int(DateJD) - 0.5
  End If
End Function

이 함수는 이 모듈에서 상당히 많이 쓰이고 있습니다. 

(5) JDtoTDT(julday)

JDtoTDT 함수는 율리우스 적일로 나타낸 시각을 UT(세계시)에서 TT(지구시)로 바꾸어 주는 함수입니다. 이 모듈에서는 달의 정삭 시각과 해의 절입 시각을 계산하는 부분에서 쓰이고 있습니다.

우리가 일상생활에서 쓰고 있는 시각 체계는 세계협정시(UTC)를 기반으로 하고 있습니다. 각 지역의 표준시는 세계협정시에서 경도 차이로 인한 시차를 더한 것입니다. 그런데 세계협정시는 원자시계를 기준으로 측정하고 있으므로 일정할 것 같지만, 지구의 운동과 어긋나지 않도록 윤초(閏秒, leap second)를 넣어서 조절하고 있어서 하루의 길이가 일정하지 않습니다. 따라서 세계협정시로 나타낸 날짜는 천문 계산을 위해서는 적합하지 않은데, 이를 보완하기 위해 천문학에서는 지구시(Terrestiral Time)라는 균일한 시간 체계를 사용합니다(이 문제에 대한 자세한 설명은 http://blog.jidolstar.com/486, http://blog.jidolstar.com/482에서 볼 수 있습니다).

Private Function JDtoTDT(ByVal julday As Double) As Double
  Dim Y As Long, T As Double, dt As Double
 
  Y = CInt(InvJDYear(julday))
  Select Case Y
   Case Is < 949
     T = (Y - 2000) / 100
     dt = (2715.6 + 573.36 * T + 46.5 * T * T) / 3600
   Case 949 To 1619
     T = (Y - 1850) / 100
     dt = (22.5 * T * T) / 3600
   Case 1620 To 1621: dt = 124 / 3600
   Case 1622 To 1623: dt = 115 / 3600
   Case 1624 To 1625: dt = 106 / 3600
   Case 1626 To 1627: dt = 98 / 3600
   Case 1628 To 1629: dt = 91 / 3600
   Case 1630 To 1631: dt = 85 / 3600
   Case 1632 To 1633: dt = 79 / 3600
   Case 1634 To 1635: dt = 74 / 3600
   Case 1636 To 1637: dt = 70 / 3600
   Case 1638 To 1639: dt = 65 / 3600
   Case 1640 To 1645: dt = 60 / 3600
   Case 1646 To 1653: dt = 50 / 3600
   Case 1654 To 1661: dt = 40 / 3600
   Case 1662 To 1671: dt = 30 / 3600
   Case 1672 To 1681: dt = 20 / 3600
   Case 1682 To 1691: dt = 10 / 3600
   Case 1692 To 1707: dt = 9 / 3600
   Case 1708 To 1717: dt = 10 / 3600
   Case 1718 To 1733: dt = 11 / 3600
   Case 1734 To 1743: dt = 12 / 3600
   Case 1744 To 1751: dt = 13 / 3600
   Case 1752 To 1757: dt = 14 / 3600
   Case 1758 To 1765: dt = 15 / 3600
   Case 1766 To 1775: dt = 16 / 3600
   Case 1776 To 1791: dt = 17 / 3600
   Case 1792 To 1795: dt = 16 / 3600
   Case 1796 To 1797: dt = 15 / 3600
   Case 1798 To 1799: dt = 14 / 3600
   Case 1800 To 1899
     T = (Y - 1900) / 100
     dt = (-0.067471 + T * (-0.058091))
     dt = (0.161416 + T * (0.145932 + T * dt))
     dt = (-0.14696 + T * (-0.146279 + T * dt))
     dt = (0.062971 + T * (0.079441 + T * dt))
     dt = (-0.012462 + T * (-0.022542 + T * dt))
     dt = (-0.000014 + T * (0.001148 + T * (0.003357 + T * dt)))
     dt = dt * 24
   Case 1900 To 1987
     T = (Y - 1900) / 100
     dt = (-0.861938 + T * (0.677066 + T * -0.212591))
     dt = (0.025184 + T * (-0.181133 + T * (0.55304 + T * dt)))
     dt = (-0.00002 + T * (0.000297 + T * dt))
     dt = dt * 24
   Case 1988 To 1996
     T = (Y - 2000) / 100
     dt = (67 + 123.5 * T + 32.5 * T * T) / 3600
   Case 1997: dt = 62 / 3600
   Case 1998 To 1999: dt = 63 / 3600
   Case 2000 To 2001: dt = 64 / 3600
   Case 2002 To 2020
     T = (Y - 2000) / 100
     dt = (63 + 123.5 * T + 32.5 * T * T) / 3600
   Case Is > 2020
     T = (Y - 1875.1) / 100
     dt = 45.39 * T * T / 3600
   Case Else
     dt = 0
  End Select
 
  JDtoTDT = dt / 24 + julday
End Function

이 함수 내부에서 쓰고 있는 dt는 천문학에서는 delta-T(ΔT)라고 부르는 값입니다. ΔT는 TT(지구시)에서 UT(세계시)를 뺀 값인데, 관측을 통해서만 알 수 있습니다. 이 함수에서 과거의 값은 역사 자료와 관측 자료를 바탕으로 추정한 값을 쓰고 있고, 미래의 값(2002년 이후)은 과거의 변화 추세를 바탕으로 예상한 값을 근사적으로 계산해서 얻은 값을 쓰고 있습니다. ΔT의 값은 연구자에 따라 계산하는 방법과 추정치가 조금씩 다릅니다. NASA의 일식 사이트에는 이 함수에서 쓴 것과는 다른 계산법을 소개하고 있는데, 이 방법으로 바꾸어 주어도 됩니다.

Polynomial Expressions for Delta T (ΔT)

http://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html

참고로 성도 프로그램과 달력 프로그램은 NASA의 일식 사이트에서 소개해 놓은 방법으로 ΔT를 계산합니다.