海上保安庁海洋情報部-天体の位置計算式をコード化してみた
海上保安庁海洋情報部のサイトには、天文・暦情報のページに天測・天測略歴の付録に掲載されているコンピュータによる天体の位置計算式が紹介されています。解説には計算方法と計算例、PDFおよびTXTには計算に使用する係数が掲載されています。これらは毎年発行されており、天体の位置計算をするならばその年に発行された係数を使用しなければなりません。
計算対象となる天体は、太陽、月、金星、火星、木星、土星、41個の明るい恒星となります。天体は限定的ですが、もともとは天測のためのものですのでしかたがありません。ですが、天体はいずれも肉眼でも見えるものですし、いずれも天文現象で観測または観望のニュースになる天体です。
今回は、紹介されているコンピュータによる天体の位置計算式の解説に紹介されている計算方法を、Fortranでコード化してみました。ほかのプログラミング言語を使用している方は、コードを書き直してみるのもいいかもしれません。
なお、プログラム中で使用している係数は太陽のものを使用しています。ほかの天体の位置を計算する場合は係数のデータを置き換えてください。
▶ サブルーチン pcob_sub.f90
subroutine PlanetaryPosition(year, month, day, hour, minute, second, pos)
! 視赤経・視赤緯・地心距離および視半径の計算
! 海上保安庁 海洋情報部による略算式(2022版)
! https:!www1.kaiho.mlit.go.jp/KOHO/index.html
!
! 年度が替わる場合は、x()およびr()の値を変更する。
! 太陽以外の天体の位置を計算する場合は、x()およびSDの値を変更する。
!
! 引 数: 観測年月日＋時分秒(日本標準時)
! year : 年
! month : 月
! day : 日
! hour : 時
! minute : 分
! second : 秒
! 戻り値: 視赤経・視赤緯・地心距離・視半径を返す
! pos(1): 視赤経
! pos(2): 視赤緯
! pos(3): 地心距離
! pos(4): 視半径
!
! 計算精度について:
! 略算式とはいえ、精度は国立天文台天文情報センター暦計算室の値と較べても
! 非常に高い一致をみる。ただ、使用にあたっては、係数の桁数の関係で視赤経
! および視赤緯、地心距離は小数点以下5-6桁までにとどめる。
implicit none
double precision, parameter :: PI = 3.141592653589793d0
double precision, parameter :: RAD = 180.0d0 / PI
double precision, parameter :: DELTAT = 70.0d0
double precision, parameter :: SD = 16.02d0 ! ベースとなる太陽の視半径
double precision, intent(in) :: year, month, day
double precision, intent(in) :: hour, minute, second
double precision, intent(out) :: pos(5)
double precision, dimension(0:29, 1:3) :: x
double precision, dimension(0:9) :: r
double precision :: f, tp, theta, rr
double precision :: ra ! 赤経
double precision :: dc ! 赤緯
double precision :: ds ! 地心距離
integer :: i, a, b, p, q, y, s, t
! 各種引数の計算
p = int(month) - 1
q = int((month + 7) / 10)
y = int((year / 4) - int(year / 4) + 0.77)
s = int(p * 0.55 - 0.33)
t = int(30 * p + q * (s - y) + p * (1 - q) + day)
f = (hour - 9.0d0) / 24.0d0 + minute / 1440.0d0 + second / 86400.0d0
tp = t + f + DELTAT / 86400.0d0
! 係数 - 太陽
if ((0 <= tp) .and. (tp <= 121)) then
! 1 / 0 - 5 / 1
a = 0
b = 121
x( 0, 1:3) = (/ 0.000000d0, 0.00000d0, 0.000000d0/)
x( 1, 1:3) = (/22.714111d0, -5.77971d0, 0.993148d0/)
x( 2, 1:3) = (/ 3.895588d0, 20.02854d0, 0.012659d0/)
x( 3, 1:3) = (/-0.102703d0, 1.72530d0, 0.002294d0/)
x( 4, 1:3) = (/ 0.034943d0, -0.98556d0, -0.000642d0/)
x( 5, 1:3) = (/ 0.006578d0, 0.01260d0, -0.000037d0/)
x( 6, 1:3) = (/-0.002269d0, 0.00796d0, 0.000013d0/)
x( 7, 1:3) = (/ 0.000090d0, -0.00351d0, 0.000009d0/)
x( 8, 1:3) = (/ 0.000090d0, 0.00081d0, 0.000005d0/)
x( 9, 1:3) = (/-0.000024d0, 0.00022d0, -0.000007d0/)
x(10, 1:3) = (/ 0.000024d0, 0.00010d0, 0.000002d0/)
x(11, 1:3) = (/ 0.000022d0, 0.00009d0, -0.000014d0/)
x(12, 1:3) = (/-0.000066d0, -0.00033d0, -0.000006d0/)
x(13, 1:3) = (/-0.000019d0, -0.00024d0, 0.000014d0/)
x(14, 1:3) = (/ 0.000040d0, 0.00022d0, 0.000004d0/)
x(15, 1:3) = (/ 0.000006d0, 0.00014d0, -0.000006d0/)
x(16, 1:3) = (/-0.000014d0, -0.00008d0, -0.000001d0/)
x(17, 1:3) = (/ 0.000004d0, -0.00003d0, 0.000002d0/)
x(18, 1:3) = (/ 0.000002d0, 0.00003d0, 0.000000d0/)
r(0) = 0.000000d0
r(1) = 0.618166d0
r(2) = 3.975436d0
r(3) = -0.000012d0
r(4) = 0.000005d0
r(5) = 0.000003d0
r(6) = 0.000000d0
r(7) = -0.000002d0
r(8) = 0.000000d0
else
if ((120 <= tp) .and. (tp <= 245)) then
! 4 / 30 - 9 / 1
a = 120
b = 244
x( 0, 1:3) = (/ 0.000000d0, 0.00000d0, 0.000000d0/)
x( 1, 1:3) = (/ 6.617773d0, 17.16569d0, 1.012365d0/)
x( 2, 1:3) = (/ 4.135644d0, -3.37392d0, 0.001158d0/)
x( 3, 1:3) = (/-0.041829d0, -5.79007d0, -0.004190d0/)
x( 4, 1:3) = (/-0.040668d0, 0.21379d0, -0.000055d0/)
x( 5, 1:3) = (/ 0.005360d0, 0.15980d0, 0.000107d0/)
x( 6, 1:3) = (/ 0.003267d0, -0.01083d0, -0.000003d0/)
x( 7, 1:3) = (/-0.000421d0, -0.00504d0, 0.000004d0/)
x( 8, 1:3) = (/-0.000120d0, 0.00078d0, -0.000006d0/)
x( 9, 1:3) = (/ 0.000017d0, 0.00005d0, -0.000010d0/)
x(10, 1:3) = (/ 0.000015d0, 0.00003d0, 0.000000d0/)
x(11, 1:3) = (/-0.000021d0, 0.00016d0, -0.000012d0/)
x(12, 1:3) = (/-0.000064d0, -0.00001d0, 0.000007d0/)
x(13, 1:3) = (/ 0.000028d0, 0.00011d0, 0.000016d0/)
x(14, 1:3) = (/ 0.000047d0, -0.00006d0, -0.000005d0/)
x(15, 1:3) = (/-0.000014d0, -0.00014d0, -0.000008d0/)
x(16, 1:3) = (/-0.000017d0, 0.00004d0, 0.000002d0/)
x(17, 1:3) = (/ 0.000005d0, 0.00004d0, 0.000002d0/)
x(18, 1:3) = (/ 0.000007d0, -0.00001d0, -0.000001d0/)
r(0) = 0.000000d0
r(1) = 18.601928d0
r(2) = 4.074040d0
r(3) = -0.000007d0
r(4) = -0.000004d0
r(5) = 0.000002d0
r(6) = 0.000000d0
r(7) = -0.000003d0
r(8) = 0.000001d0
else
if ((243 <= tp) .and. (tp <= 366)) then
! 8 / 31 - 12 / 32
a = 243
b = 366
x( 0, 1:3) = (/ 0.000000d0, 0.00000d0, 0.000000d0/)
x( 1, 1:3) = (/14.543576d0, -10.64695d0, 0.994687d0/)
x( 2, 1:3) = (/ 4.052451d0, -16.83097d0, -0.013825d0/)
x( 3, 1:3) = (/ 0.15169d0, 3.51949d0, 0.001824d0/)
x( 4, 1:3) = (/ 0.013092d0, 0.97356d0, 0.000717d0/)
x( 5, 1:3) = (/-0.013453d0, -0.02695d0, -0.000036d0/)
x( 6, 1:3) = (/-0.002111d0, -0.02322d0, -0.000021d0/)
x( 7, 1:3) = (/ 0.000314d0, -0.00527d0, 0.000001d0/)
x( 8, 1:3) = (/ 0.000181d0, 0.00005d0, -0.000014d0/)
x( 9, 1:3) = (/-0.000011d0, 0.00051d0, -0.000002d0/)
x(10, 1:3) = (/-0.000010d0, -0.00011d0, -0.000004d0/)
x(11, 1:3) = (/-0.000056d0, 0.00019d0, -0.000003d0/)
x(12, 1:3) = (/-0.000016d0, 0.00011d0, 0.000018d0/)
x(13, 1:3) = (/ 0.000057d0, -0.00036d0, 0.000003d0/)
x(14, 1:3) = (/ 0.000011d0, -0.00001d0, -0.000012d0/)
x(15, 1:3) = (/-0.000025d0, 0.00022d0, -0.000001d0/)
x(16, 1:3) = (/-0.000003d0, 0.00000d0, 0.000004d0/)
x(17, 1:3) = (/ 0.000011d0, -0.00009d0, 0.000001d0/)
x(18, 1:3) = (/ 0.000003d0, 0.00001d0, -0.000001d0/)
r(0) = 0.000000d0
r(1) = 2.651397d0
r(2) = 4.041161d0
r(3) = 0.000016d0
r(4) = 0.000003d0
r(5) = -0.000001d0
r(6) = 0.000002d0
r(7) = -0.000002d0
r(8) = 0.000002d0
end if
end if
end if
! 視赤経・視赤緯および地心距離の計算
theta = acos((2.0d0 * tp - (a + b)) / (b - a)) * RAD
ra = 0.0d0
dc = 0.0d0
ds = 0.0d0
do i = 0, 17
ra = ra + (x(i + 1, 1) * cos(mod((theta * i), 360.0d0) / RAD))
dc = dc + (x(i + 1, 2) * cos(mod((theta * i), 360.0d0) / RAD))
ds = ds + (x(i + 1, 3) * cos(mod((theta * i), 360.0d0) / RAD))
end do
if (ra >= 24) then
ra = ra - 24
else
if (ra < 0) then
ra = ra + 24
end if
end if
pos(1) = ra
pos(2) = dc
pos(3) = ds
pos(4) = SD / ds
! Rの計算
theta = acos((2.0d0 * (t + f) - (a + b)) / (b - a)) * RAD
rr = 0.0d0
do i = 0, 7
rr = rr + (r(i + 1) * cos(mod((theta * i), 360.0d0) / RAD))
end do
ra = rr - ra
pos(5) = ra
return
end subroutine PlanetaryPosition▶ メインルーチン pcob_main.f90
テストのためのメインルーチンです。Rから求めた赤経は使用していません。
program pcob_main
! 天体の視位置を求める
implicit none
double precision, dimension(5) :: pos
double precision :: dy, dt
double precision :: year, month, day, hour, minute, second
double precision :: rs, ds, ss
integer :: rh, rm, dh, dm, sm
character(1) :: sign
character(255) :: sfmt
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(*, *)
sfmt = '(i5, "年 ", i2, "月 ", i2, "日 ", i2, "時 ", i2, "分 ", i2, "秒")'
write(*, sfmt) int(year), int(month), int(day), int(hour), int(minute), int(second)
call PlanetaryPosition(year, month, day, hour, minute, second, pos)
rh = int(pos(1))
rs = 60.0d0 * (pos(1) - rh)
rm = int(rs)
rs = 60.0d0 * (rs - rm)
sign = '+'
if (pos(2) < 0) then
sign = '-'
end if
dh = abs(int(pos(2)))
ds = 60.0d0 * (abs(pos(2)) - dh)
dm = int(ds)
ds = 60.0d0 * (ds - dm)
sm = int(pos(4))
ss = 60.0d0 * (pos(4) - sm)
write(*, *)
write(*, *) ' 赤経 赤緯 地心距離 視半径'
write(*, *) ' h m s 。 , AU , ,,'
write(*, '(" ", i2, " ", i2, " ", f6.3, " ", a1, i2, " ", i2, " ", f5.2, " ", f9.7, " ", i2, " ", f5.2)') &
rh, rm, rs, sign, dh, dm, ds, pos(3), sm, ss
write(*, *)
stop
end program pcob_mainコンパイル
gfortran pcob_main.f90 pcob_sub.f90実行
./a.out実行すると日付と時刻を入力するよう求められます。日付が2022年５月４日なら20220504、時刻が15時24分37秒なら152437とします。時刻は必ず秒まで入力してください。なお、日付と時刻はカンマ（,）で区切って入力してください（実行結果参照）。
実行結果
DATE AND TIME(JST) ? 20220504,152437
2022年 5月 4日 15時 24分 37秒
赤経 赤緯 地心距離 視半径
h m s 。 , AU , ,,
2 45 19.205 +15 58 34.87 1.0082528 15 53.33国立天文台 天文情報センター暦計算室（暦年表ー太陽の地心座標）の値は以下のようになっています（時刻系の指定は+09時間(日本)）。
h m s 。 , ,, AU , ,,
2 45 19.208 15 58 34.89 1.0082504 15 53.31 