天体の位置計算 – グリニッジ視恒星時の計算（1950）

「章動による赤経・赤緯の変化」 66頁

章動表から黄経の章動と黄道傾斜角の章動を求めることができたら、分点差を求めることができます。分点差を求めることができたなら、グリニッジ平均恒星時に分点差を加算することで、グリニッジ視恒星時を求めることができます。ただし、分点差を求めるには、章動表を使って計算するのですが、「天体の位置計算 増補版」に掲載されている項目数は69項目、1984年で106項目、最新の2021年のデータでは、678項目もあります。また、係数も11もあるために、計算量はかなり多くなります。計算する章動は、それぞれ日月章動と惑星章動を求めなければなりません。方法がわかれば最新のデータを利用することで、国立天文台天文情報センター暦計算室で求められたグリニッジ視恒星時とかなり近い値を求めることができるのではないかと思います。こういった計算は、結構大変なのですね。

計算式
分点差Eqは、黄経の章動⊿φと黄道傾斜角の章動⊿εから下記の計算式で求めることができます。

\[ E_{q} = \Delta\phi + \cos(\epsilon + \Delta\epsilon) \]

最新データ
新しいデータは以下で入手することができます。

「天体の位置計算 増補版」
「天体の位置計算 増補版」 229 – 232頁に新しい計算式と章動表が掲載されています。関数iauNut50のxplと計算式のel、elp、f、d、omの式の差し替えを行えば新しいデータで計算することができます。

Sofa（Standards of Fundamental Astronomy ）
Fortran77とCのソースコードが公開されています。GNU Fortranで移植を考えるのであれば、Fortran77のソースコードを利用するとよいでしょう。アーカイブをダウンロードして、sofa→20210512→f77→srcフォルダのnut00a.forまたはnut80.forがそれにあたります。

▶ メインルーチン pcc_gast.f90

program pcc_gast
    use pcclib

    implicit none
    
    double precision :: dy, dt, julian
    double precision :: year, month, day, hour, minute, second
    double precision :: t0, t1
    double precision :: msdt, sdt_hour, sdt_minute, sdt_second
    double precision :: dp, de, e, eq
    double precision :: sdt(3), eps(2)
    character(255) :: dfmt

    ! 日付の入力
    write(*,*)
    write(*, fmt='(a)', advance='no') 'DATE ? '
    read(*,*) dy

    ! 日付と時刻を抽出
    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)

    ! 時刻引数の計算
    t0 = julian - 2415020d0
    t1 = t0 / 36525d0
    
    ! グリニッジ視恒星時の計算
    e = iauMOE50(t1)                            ! 平均黄道傾斜角
    e = deg2rad(e)
    msdt = MeanSDT50(t1)                        ! グリニッジ平均恒星時（0 UTC）
    eps = iauNut50(t0, t1)                      ! 章動 ⊿Φと⊿εを求める
    dp = eps(1)                                 ! ⊿Φ
    de = eps(2) / 3600d0                        ! ⊿ε
    eq = dp * cos(deg2rad(de) + e)              ! 分点差
    sdt = deg2hms(rad2deg(msdt))
    sdt_hour = sdt(1)
    sdt_minute = sdt(2)
    sdt_second = sdt(3) + (eq / 15d0)           ! 分点差を時刻の秒に加算

    ! 恒星時が範囲を超えたときの処理
    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   Greenwich 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_gast

コンパイル

gfortran pcc_lib.f90 pcc_gast.f90

実行

./a.out

実行結果

DATE ? 19780610

 Date / 0h UTC           Julian Date   Greenwich Sidereal Time
 -------------------------------------------------------------
                                              h   m   s
 1978  6 10.00000       2443669.50000       17  11  58.714

DATE ? 19780620         

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

57