天体の位置計算 – グリニッジ平均恒星時の計算

「天球座標-恒星時」26-28頁

グリニッジ平均恒星時は春分点を基準とする恒星時で歳差の影響を考慮したものになります。本来は歳差と章動の影響を含めて計算するのですが、計算量の多さと複雑さで大変な計算となります。天文計算で使用する恒星時は、グリニッジ視恒星時との差は１″しかなく、よほど精密な計算をしない限り、グリニッジ平均恒星時でも充分です。

今回は、グリニッジ平均恒星時を求めます。「天体の位置計算 増補版」では、新しい計算システムとして、新しい計算式も掲載されていますので、それも併せて紹介して計算していきたいと思います。

　

グリニッジ平均恒星時を求める計算式
26頁の式2-5を使用して求めることができます。サブルーチン化するにあたり、これを角度に換算しています。

\[ \theta_{0} = 6^{h}38^{m}45^{s}.836 + 8640184^{s}.542T_{u} + 0^{s}.0929T_{u}^{2}\\ \]

\(T_{u}\) は1900年１月０日正午から計算日までの経過日数を36525日を単位として求めた値になります。

\[ T_{u} = \frac{(JD\,-\,2415020)}{36525} \]

▶ サブルーチン MeanSDT50

    function MeanSDT50(t) result(msdt)
    ! グリニッジ平均恒星時(1950)
    ! t = (JD - 2415020) / 36525
    
        double precision, intent(in) :: t
        double precision :: msdt
    
        msdt = Mla(99.690983333d0 + 36000.768925d0 * t + 0.0003870833d0 * t**2)
    end function MeanSDT50

新しいグリニッジ平均恒星時を求める計算式
235頁式9-6、この式もサブルーチン化するにあたり、角度に換算しています。

\[ \theta_{0} = 6^{h}41^{m}50^{s}.54841 + 8640184^{s}.812866T + 0^{s}.093104T^{2}\,-\,0^{s}.0000062T^{3} \]

Tは2000年１月１日正午から計算日の世界時０時のユリウス通日（JD)を36525日を単位として求めた値になります。こちらは「天体の位置計算 増補版」に掲載されているような経過日数を求める計算式はなく、必ず世界０時ユリウス通日を求めなければなりません

\[ T = \frac{(JD – 2451545)}{36525} \]

▶ サブルーチン MeanSDT80

    function MeanSDT80(t) result(msdt)
    ! グリニッジ平均恒星時(1984)
    ! t = (JD - 2451545) / 36525
    
            double precision, intent(in) :: t
            double precision :: msdt
    
            msdt = Mla(100.460618375d0 + 36000.770053608d0 * t + 0.000387933d0 * t**2 - 0.00000002583d0 * t**3)
    end function MeanSDT80

Fortranによるコード

▶ メインルーチン pcc_msdt.f90

program pcc_msdt
    use pcclib
    
    implicit none
    
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: RAD = 180.0d0 / PI

    double precision :: julian, dy, dt
    double precision :: year, month, day, hour, minute, second
    double precision :: sdt_hour, sdt_minute, sdt_second
    double precision :: msdt(3), dtime(3)
    double precision :: t
    character(255) :: dfmt

    write(*,*)
    write(*, fmt='(a)', advance='no') 'DATE AND TIME(JST) ? '
    read(*,*) dy, dt

    year = int(dy / 10000.0d0)
    month = mod(int(dy / 100.0d0), 100)
    day = mod(dy, 100.0d0)
    hour = int(dt / 10000.0d0)
    minute = mod(int(dt / 100.0d0), 100)
    second = mod(dt, 100.0d0)

    ! ユリウス日の計算
    julian = julianDay(year, month, day)
    julian = julian + (hour / 24d0 + minute / 1440d0 + second / 86400d0)
    julian = julian - 0.375

    ! ユリウス日から年月日を求める
    dtime = Julian2date(julian)
    year = dtime(1)
    month = dtime(2)
    day = dtime(3)
    
    ! 時刻引数の計算
    ! t = (julian - 2451545d0) / 36525d0
    t = (julian - 2415020d0) / 36525d0
    
    ! グリニッジ平均恒星時の計算
    ! msdt = deg2hms(rad2deg(MeanSDT80(t)))
    msdt = deg2hms(rad2deg(MeanSDT50(t)))
    sdt_hour = msdt(1)
    sdt_minute = msdt(2)
    sdt_second = msdt(3)

    ! 恒星時が範囲を超えたときの処理
    if (sdt_hour > 24d0) then
        sdt_hour = sdt_hour - 24d0
    end if
    if (sdt_second >= 60d0) then
        sdt_second = sdt_second - 60d0
        sdt_minute = sdt_minute + 1d0
    end if
    if (sdt_minute == 60d0) then
        sdt_minute = 0d0
        sdt_hour = sdt_hour + 1d0
    end if
    if (sdt_hour == 24d0) then
        sdt_hour = 0d0
    end if

    ! 出力
    dfmt = '(" ", i4, " ", i2, " ", f8.5, "      ", f14.5, "         ", i2, "  ", i2, "  ", f6.3)'
    write(*, *)
    write(*, *) 'Date / 0h UTC           Julian Date        Mean Sidereal Time'
    write(*, *) '-------------------------------------------------------------'
    write(*, *) '                                               h   m   s'
    write(*, dfmt) int(year), int(month), day, julian, int(sdt_hour), int(sdt_minute), sdt_second
    write(*, *)

    stop
end program pcc_msdt

コンパイル

gfortran pcc_lib.f90 pcc_msdt.f90

実行

./a.out

実行結果

DATE AND TIME(JST) ? 19780620,090000

 Date / 0h UTC           Julian Date        Mean Sidereal Time
 -------------------------------------------------------------
                                                h   m   s
 1978  6 20.00000       2443679.50000         17  51  24.247

57