天体の位置計算 – 黄道傾斜角の計算

「平均黄道傾斜角の計算式」 62, 229頁
「シリーズ現代の天文学 13 天体の位置と運動」 64頁

黄道傾斜角には、平均黄道傾斜角と真黄道傾斜角の２つがあります。主に位置計算に用いられるのは平均黄道傾斜角です。平均黄道傾斜角を求めるには、計算式があるので簡単に求められます。ですが、真黄道傾斜角の方は求めるには大変な計算量が必要です。

今回は、平均黄道傾斜角と真黄道傾斜角を求めてみます。真黄道傾斜角は、黄道章動表から求めています。この表は１０５６ものデータがあり、それらを「天体の位置計算 増補版」６３頁にあるように計算して求めなければなりません。コードにすると少し大変なので、Excelの表で計算したものを紹介します。この表はダウンロード可能となっています。ぜひ、国立天文台天文情報センター暦計算室の結果と較べてみてください。なお、⊿εはおまけです。

　

Fortranによるコード

▶ サブルーチン iauMOE50

    function iauMOE50(t) result(moe)
    ! 平均黄道傾斜角(1950)
    ! t = (JD - 2415020.0) / 36525.0
    
        double precision, intent(in) :: t
        double precision :: moe
    
        moe = 23.45229444d0 - 0.0130125d0 * t - 0.0000016389d0 * t**2 + 0.00000050278d0 * t**3
    end function iauMOE50

▶ サブルーチン iauMOE80

    function iauMOE80(t) result(moe)
    ! 平均黄道傾斜角(1984)
    ! t = (JD - 2451545.0) / 36525.0
    
        double precision, intent(in) :: t
        double precision :: moe
            
        moe = 84381.448d0 - 46.8150d0 * t - 0.00059d0 * t**2 + 0.001813d0 * t**3
        moe = moe / 3600d0
    end function iauMOE80

▶ サブルーチン iauMOE

    function iauMOE(t) result(moe)
    ! 平均黄道傾斜角
    ! t = (JD - 2451545.0) / 36525.0
    
        double precision, intent(in) :: t
        double precision :: moe
    
        moe = 84381.406d0 - 46.836769d0 * t - 0.0001831d0 * t**2 + 0.0020034d0 * t**3 &
            - 0.000000576d0 * t**4 - 0.0000000434d0 * t**5
        moe = moe / 3600d0
    end function iauMOE

▶ メインルーチン pcc_moe.f90

program moe
    use pcclib

    implicit none

    double precision :: julian, t
    double precision :: dy, dt
    double precision :: year, month, day
    double precision :: hour, minute, second
    double precision :: ecl
    double precision :: dtime(3)
    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 - 2415020d0) / 36525d0
    t = (julian - 2451545d0) / 36525d0

    !平均黄道傾斜角の計算
    ! ecl = iauMOE50(t)
    ! ecl = iauMOE80(t)
    ecl = iauMOE(t)

    ! 出力
    write(*, *)
    dfmt = '(" ", i4, " ", i2, " ", f8.5, "       ", f14.5, "               ", f15.12)'
    write(*, *)
    write(*, *) 'Date / UTC               Julian Date           Mean obliquity of ecliptic'
    write(*, *) '-------------------------------------------------------------------------'
    write(*, *) '                                                      。'
    write(*, dfmt) int(year), int(month), day, julian, ecl
    write(*, *)

    stop
end program moe

コンパイル

gfortran pcc_lib.f90 pcc_moe.f90

実行

./a.out

実行結果

DATE AND TIME(JST) ? 20221024,210000


 Date / UTC               Julian Date           Mean obliquity of ecliptic
 -------------------------------------------------------------------------
                                                       。
 2022 10 24.50000        2459877.00000               23.4363116

国立天文台天文情報センター暦計算室の結果と較べても一致していることがわかります。

57