天体の位置計算 – 続・グリニッジ平均恒星時の計算
前記事「天体の位置計算-グリニッジ平均恒星時の計算」では、「天体の位置計算 増補版」235頁にあるJ2000.0の式を紹介しましたが、国立天文台天文情報センター暦計算室の歴表年表の値と較べると、0s.002〜3ほどの差がでます。この値でも構わないのですが、もう少し精度を上げるために、「シリーズ現代の天文学13 天体の位置と運動 第2版」（日本評論社 2017 福島登志夫[編]）66頁に紹介されている式を使って求めることにします。解説については本書を参照してください。
計算式
\[ \begin{array}{l} ERA = 2\pi(0.7790572732640 + 1.00273781191135448(JD\,-\,2451545.0)\\ GMST = ERA + 0^{″}.014506 + 4612^{″}.156534T+1^{″}.3915817T^{2}\\ \hspace{32pt}\,-\,0^{″}.00000044T^{3}\,-\,0^{″}.000029956T^{4}\,-\,0^{″}.0000000368T^{5}\\ \\ JD : ユリウス日\\ T : (JD\,-\,2451545.0)/36525.0 \end{array} \]
計算は、ERAを計算したら、GMSTを計算しますが、このときERAの値は加算しません。別々に求めた後にERAとGMSTを加算します。最初、この式を見たときにどうすればよいかわからず、長い間放置していましたが、計算方法が「マイコン天文学Ⅰ」（恒星社 1983 中野主一著）の49頁にある地方恒星時の計算にありました。コードはN88-BASICで書かれたものですが、これが解決方法になりました。
ERAとGMSTの2つの式の違いは単位がERAはラジアン値、GMSTは秒だということです。ですからGMSTの計算結果を角度に変換し、加算するときにERAを角度に変換しています。
▶ サブルーチン MeanSDT
function MeanSDT(julian, t) result(gmst)
! グリニッジ平均恒星時
! t = (JD - 2451545.0) / 36525.0
double precision, parameter :: PI = 3.141592653589793238462643d0
double precision, parameter :: PI2 = 6.283185307179586476925287d0
double precision, parameter :: RAD = 180d0 / PI
double precision, intent(in) :: julian, t
double precision :: gmst, era
era = 0.7790572732640d0 + 1.00273781191135448d0 * (julian - 2451545.0d0)
era = PI2 * (era - int(era))
gmst = 0.014506d0 + 4612.156534d0 * t + 1.3915817d0 * t**2 - 0.00000044d0 * t**3 &
- 0.000029956d0 * t**4 - 0.0000000368d0 * t**5
gmst = gmst / 3600d0
gmst = ((era * RAD) + gmst) / 15d0
end function MeanSDTFortranによるコード
▶ メインルーチン pcc_msdt_2.f90
program pcc_msdt_2
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
! グリニッジ平均恒星時の計算
msdt = deg2dms(MeanSDT(julian, 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_2コンパイル
gfortran pcc_lib.f90 pcc_msdt_2.f90実行
./a.out実行結果
DATE AND TIME(JST) ? 20221023,090000
Date / 0h UTC Julian Date Mean Sidereal Time
-------------------------------------------------------------
h m s
2022 10 23.00000 2459875.50000 2 5 35.042国立天文台天文情報センター暦計算室の歴表年表と比較する場合は、時刻系を「+09時間(日本)」にセットしてください。これは、時刻の入力が日本標準時刻（JST）として入力しているからです。ユリウス日およびグリニッジ平均恒星時は日本標準時刻 – 9時間として計算されています。ちなみに、国立天文台天文情報センター暦計算室の歴表年表の時刻を０時、時刻系を「+00時間(世界時)」に変更して計算しても同じ結果になります。
