天体の位置計算 - ユリウス日から年月日の計算

天体の位置計算 – ユリウス日から年月日の計算

「マイコン宇宙講座」には、jdateというサブルーチンがあります。このサブルーチンはユリウス日から年月日を求めるものですが、ここでは、「新装改訂版天文計算入門」（恒星社 1996 新装版第1刷長谷川一郎著）62頁「JDからグレゴリオ暦法への換算」に掲載されている方法を使って求めたいと思います。

計算手順（16.6）が掲載されていますので、そのとおりに計算していけば簡単に求めることができます。計算は電卓を使っても求めることができるようになっており、ややプログラミング的ではないですが、簡便ですのでコード化はしやすいと思いますし、Excelなどの表計算ソフトを使ってもいいと思います。

計算式
ユリウス日は世界時になりますので、求められる年月日も世界時になります。たとえば、日本時間で2022年７月７日09時00分00秒のユリウス日は2459767.5 UTCとなり、そこから求められる年月日は2022年7月7日00時00分00秒UTCとなります。

\[ \begin{eqnarray} A &=& [JD + 68569.5]\\ B &=& [A \div 36524.25]\\ C &=& A – [36524.25B + 0.75]\\ E &=& [(C + 1) \div 365.25025]\\ F &=& C – [365.25E] + 31\\ G &=& [F \div 30.59]\\ D &=& F – [30.59G] + (JD + 0.5の少数部分)\\ H &=& [G + 11]\\ M &=& G – 12H + 2\\ Y &=& 100(B – 49) + E + H\\ \end{eqnarray} \] \[ \begin{array}{l} ※[]は整数部を取る\\ ※引用元:「新装改訂版天文計算入門」63頁式16.6より \end{array} \]

計算式のY、M、Dが年月日になります。もし、12月32日になった場合は、翌年の1月1日とします。また、Dはユリウス日によっては小数点がつくことがあり、小数点部分は時刻となります。時刻の計算は以下のようになります。

\[ \begin{eqnarray} A &=& 24(D – [D])\\ H &=& [A]\\ B &=& 60(A – H)\\ M &=& [B]\\ S &=& [60(B- M)] \end{eqnarray} \] \[ \begin{array}{l} ※[]は整数部を取る\\ \end{array} \]

Fortranによるコード

▶ サブルーチン Julian2date_2

    function Julian2date_2(julian) result(dtime)
    ! ユリウス日から年月日を求める
    !「天文計算入門」63頁 式(16.6)
        
        double precision, intent(in) :: julian
        double precision :: year, month, day
        double precision :: dtime(3)
        double precision :: a, b, c, e, f, g, h

        a = int(julian + 68569.5d0)
        b = int(a / 36524.25d0)
        c = a - int(36524.25d0 * b + 0.75d0)
        e = int((c + 1) / 365.25025d0)
        f = c - int(365.25d0 * e) + 31
        g = int(f / 30.59d0)
        day = f - int(30.59d0 * g) + (julian + 0.5d0) - int(julian + 0.5d0)
        h = int(g / 11d0)
        month = int(g - 12d0 * h + 2d0)
        year = int(100 * (b - 49d0) + e + h)
        
        ! 12月32日になったときの処理
        if (int(month) == 12d0) then
            if (int(day) > 31d0) then
                year = year + 1d0
                month = (month + 1d0) - 12d0
                day = 1d0
            end if
        end if

        dtime(1) = year
        dtime(2) = month
        dtime(3) = day
    end function Julian2date_2

▶ メインルーチン pcc_jdate.f90

32行目にコメントをつけていますが、コメントアウトすることで、UTCの年月日および時刻を返します。実行結果はコメントアウトしています。

program jdate
    use pcclib

    implicit none

    double precision :: julian
    double precision :: dy, dt
    double precision :: year, month, day, day_temp
    double precision :: hour, minute, second
    double precision :: dtime(3)
    character(128) :: 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)

    ! 入力した年月日を表示
    write(*, *)
    dfmt = '(" INPUT DATE = ", i4, "年 ", i2, "月 ", i2, "日 ", i2, "時 ", i2, "分 ", i2, "秒")'
    write(*, dfmt) int(year), int(month), int(day), int(hour), int(minute), int(second)

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

    ! ユリウス日から年月日を求める
    dtime = Julian2date_2(julian)
    year = dtime(1)
    month = dtime(2)
    day = dtime(3)

    ! 日の小数点から時刻を求める
    day_temp = 24.0d0 * (day - int(day))
    hour = int(day_temp)
    second = 60d0 * (day_temp - hour)
    minute = int(second)
    second = int(60d0 * (second - minute))

    ! 出力
    write(*, *)
    dfmt = '("RETURN DATE = ", i4, "年 ", i2, "月 ", i2, "日 ", i2, "時 ", i2, "分 ", i2, "秒")'
    write(*, dfmt) int(year), int(month), int(day), int(hour), int(minute), int(second)
    write(*, *)

    stop
end program jdate

コンパイル

gfortran pcc_lib.f90 pcc_jdate.f90

実行

./a.out

実行結果
うるう年が考慮されているかをチェックした実行結果です。きちんと考慮されていますね。

DATE AND TIME(JST) ? 20000301,063045

 INPUT DATE = 2000年  3月  1日  6時 30分 45秒

RETURN DATE = 2000年  2月 29日 21時 30分 45秒

DATE AND TIME(JST) ? 20220301,063045

 INPUT DATE = 2022年  3月  1日  6時 30分 45秒

RETURN DATE = 2022年  2月 28日 21時 30分 45秒

老猫こてつ