海上保安庁海洋情報部-続・天体の位置計算式をコード化してみた
前記事「海上保安庁海洋情報部-天体の位置計算式をコード化してみた」では、位置を求める係数をプログラム中に記述していました。今回は、前記事「海上保安庁海洋情報部-天体の位置計算データファイルの作成」で作成したデータファイルを利用して、日付と惑星番号を選択するだけで位置を求めるようにコードを変更しました。なお、恒星については対応していません。近日中に公開する予定です。
CSVファイルについての注意
データファイルはCSV形式で保存したものを読み込んでいます。ですので、あらかじめ表計算のワークシートごとにCSV形式で出力して保存しておく必要があります。このとき、1行目には、ダミーデータとして０を追加してください（見やすいように加工をしています）。
0,0,0 22.714111, -5.77971, 0.993148 3.895588, 20.02854, 0.012659 -0.102703, 1.72530, 0.002294 0.034943, -0.98556, -0.000642 0.006578, 0.01260, -0.000037 -0.002269, 0.00796, 0.000013 0.00009, -0.00351, 0.000009 0.00009, 0.00081, 0.000005 -0.000024, 0.00022, -0.000007 0.000024, 0.00010, 0.000002 0.000022, 0.00009, -0.000014 -0.000066, -0.00033, -0.000006 -0.000019, -0.00024, 0.000014 0.00004, 0.00022, 0.000004 0.000006, 0.00014, -0.000006 -0.000014, -0.00008, -0.000001 0.000004, -0.00003, 0.000002 0.000002, 0.00003, 0.000000
CSV形式のファイル名は以下の形式になるようにします。
惑星名-data-番号.csv
たとえば、太陽でしたら1月0日から5月1日の係数データのファイル名は、sun-data-01.csvになります。次はsun-data-02.csv、sun-data-03となります。この形式はプログラム中でも利用していますので、変更する際はコードを書き換えてください。月は一月ごとになりますので、moon-data-01.csvからmoon-data-12.csvまでとなります。
CSVファイルを保管するフォルダ名とディレクトリ
CSVファイルを格納するディレクトリ構造は以下の図のとおりです。よけいなファイルがありますが気にしないでください。必要なのは、pcob_main.f90とpcob_sub.f90だけです。
親ディレクトリ名は何でも構いませんが、data、2022と各惑星名のフォルダ名は変更しないようにしてください。大文字小文字混在もNGです。ディレクトリ構造も変更しないでください。もし変更する場合は、コードも対応したものに書き換えてください。
Fortranによるコード
▶ サブルーチン pcob_sub.f90
前記事「海上保安庁海洋情報部-天体の位置計算式をコード化してみた」で作成したファイルの名前を変更したものです。ですが、コードは削除と追加を行っています。言うまでもないことですが、係数データは、その年で有効なものです。過去や未来の日付を入力しても正しく求めることができません。計算する日付（年）の係数データを用意する必要があります。
subroutine getPlanetaryPosition(planet_number, year, month, day, hour, minute, second, pos)
! 視赤経・視赤緯・地心距離および視半径の計算
! 海上保安庁 海洋情報部による略算式
! https:!www1.kaiho.mlit.go.jp/KOHO/index.html
!
! 引 数: 観測年月日＋時分秒(日本標準時)
! year : 年
! month : 月
! day : 日
! hour : 時
! minute : 分
! second : 秒
! 戻り値: 視赤経・視赤緯・地心距離・視半径を返す
! pos(1): 視赤経
! pos(2): 視赤緯
! pos(3): 地心距離
! pos(4): 視半径
! pos(5): R
! pos(6): 平均黄道傾斜角
!
! 計算精度について:
! 略算式とはいえ、精度は国立天文台天文情報センター暦計算室の値と較べても
! 非常に高い一致をみる。ただ、使用にあたっては、係数の桁数の関係で視赤経
! および視赤緯、地心距離は小数点以下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, dimension(6), parameter :: SO = (/0.0d0, 16.02d0, 0.1383d0, 0.078d0, 1.535d0, 1.23d0/) ! ベースとなる太陽の視半径
character(7), parameter :: CELESTIA(6) = (/'moon ', 'sun ', 'venus ', 'mars ', 'jupiter', 'saturn '/)
double precision, intent(in) :: year, month, day
double precision, intent(in) :: hour, minute, second
double precision, intent(out) :: pos(6)
integer, intent(in) :: planet_number
double precision, dimension(12, 2) :: check = 0.0d0
double precision, dimension(0:30, 1:3) :: x = 0.0d0
double precision, dimension(0:10, 1:2) :: r = 0.0d0
double precision :: f, tp, theta, rr
double precision :: ra ! 赤経
double precision :: dc ! 赤緯
double precision :: ds ! 地心距離
double precision :: ec ! 黄道傾斜角
integer :: i, n, a, b, p, q, y, s, t
character(20) :: path
character(19) :: filename, rfilename
character(2) :: filenumber
character(4) :: syear
! 各種引数の計算
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
! 惑星番号によって読み込むファイル名を構成する
write(syear, '(i4)') int(year)
path = './data/'//syear//'/'//trim(CELESTIA(planet_number))//'/'
select case(planet_number)
case(1) ! 月
n = 30
check( 1, 1:2) = (/ 0.0d0, 32.0d0/)
check( 2, 1:2) = (/ 31.0d0, 60.0d0/)
check( 3, 1:2) = (/ 59.0d0, 91.0d0/)
check( 4, 1:2) = (/ 90.0d0, 121.0d0/)
check( 5, 1:2) = (/120.0d0, 152.0d0/)
check( 6, 1:2) = (/151.0d0, 182.0d0/)
check( 7, 1:2) = (/181.0d0, 213.0d0/)
check( 8, 1:2) = (/212.0d0, 244.0d0/)
check( 9, 1:2) = (/243.0d0, 274.0d0/)
check(10, 1:2) = (/273.0d0, 305.0d0/)
check(11, 1:2) = (/304.0d0, 335.0d0/)
check(12, 1:2) = (/334.0d0, 366.0d0/)
do i = 1, 12
if ((check(i, 1) <= tp) .and. (tp <= check(i, 2))) then
a = int(check(i,1))
b = int(check(i, 2))
write(filenumber, '(i2.2)') i
filename = trim(CELESTIA(planet_number))//'-data-'//filenumber//'.csv'
exit
end if
end do
case(2:6) ! 太陽・金星・火星・木星・土星
n = 18
check(1, 1:2) = (/ 0.0d0, 121.0d0/)
check(2, 1:2) = (/120.0d0, 244.0d0/)
check(3, 1:2) = (/243.0d0, 366.0d0/)
do i = 1, 3
if ((check(i, 1) <= tp) .and. (tp <= check(i, 2))) then
a = int(check(i,1))
b = int(check(i, 2))
write(filenumber, '(i2.2)') i
filename = trim(CELESTIA(planet_number))//'-data-'//filenumber//'.csv'
rfilename = 'r-data-'//filenumber//'.csv'
exit
end if
end do
end select
! 構成されたファイル名から係数データを読み込む
write(*, *)
open(10, file=trim(path)//trim(filename), status="old")
do i = 0, n
read(10, *) x(i, 1:3)
end do
close(10)
write(*, *)
! Rおよび黄道傾斜角の係数を読み込む(太陽・金星・火星・木星・土星)
if ((2 <= planet_number) .and. (planet_number <= 6)) then
open(20, file='./data/'//syear//'/R/'//trim(rfilename), status="old")
do i = 0, 8
read(20, *) r(i, 1:2)
end do
close(20)
end if
write(*, *)
! 視赤経・視赤緯および地心距離の計算
theta = acos((2.0d0 * tp - (a + b)) / (b - a)) * RAD
ra = 0.0d0
dc = 0.0d0
ds = 0.0d0
do i = 0, n
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)) ! 月の場合はH.P.
end do
if (int(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
! Rと黄道傾斜角の計算(太陽・金星・火星・木星・土星)
if ((2 <= planet_number) .and. (planet_number <= 6)) then
theta = acos((2.0d0 * (t + f) - (a + b)) / (b - a)) * RAD
rr = 0.0d0
ec = 0.0d0
do i = 0, 7
rr = rr + (r(i + 1, 1) * cos(mod((theta * i), 360.0d0) / RAD))
ec = ec + (r(i + 1, 2) * cos(mod((theta * i), 360.0d0) / RAD))
end do
ra = rr - ra
pos(4) = SO(planet_number) / ds
pos(5) = ra ! Rを考慮した赤経
pos(6) = ec ! 黄道傾斜角
else
! 月のときは０が入る
pos(3) = ds
ds = sin(ds / RAD)
ds = 0.2725d0 * ds
ds = asin(ds) * RAD
pos(4) = ds * 60.0d0
pos(5) = 0.0d0
pos(6) = 0.0d0
end if
return
end subroutine getPlanetaryPosition
▶ メインルーチン pcob_main.f90
program pcob_main
! 天体の視位置を求める
implicit none
character(7), parameter :: CELESTIA(6) = (/'Moon ', 'Sun ', 'Venus ', 'Mars ', 'Jupiter', 'Saturn '/)
double precision, dimension(6) :: pos
double precision :: year, month, day, hour, minute, second
double precision :: rs, ds, ss
integer :: rh, rm, dh, dm, sm
character(1) :: sign
character(132) :: sfmt
double precision :: dy = 0.0d0
double precision :: dt = 0.0d0
integer :: planet_number = 0
write(*, *)
write(*, fmt='(a)', advance='no') 'DATE AND TIME(JST) ? '
read(*, *) dy, dt
write(*, fmt='(a)', advance='no') '(1: Moon, 2: Sun, 3:Venus, 4: Mars, 5: Jupiter, 6: Saturn) ? '
read(*, *) planet_number
if ((planet_number < 1) .or. (planet_number > 6)) then
stop 'Bad Planet Number!'
end if
! 年月日および時刻を抽出する
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)
call getPlanetaryPosition(planet_number, 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)
write(*, '(i5, "年 ", i2, "月 ", i2, "日 ", i2, "時 ", i2, "分 ", i2, "秒")') &
int(year), int(month), int(day), int(hour), int(minute), int(second)
write(*, *)
write(*, '(a, "の視位置")') trim(CELESTIA(planet_number))
write(*, *)
if (planet_number == 1) then
! 月の場合は黄道傾斜角を表示しない
sm = int(pos(4))
ss = 60.0d0 * (pos(4) - sm)
sfmt = '(i2, " ", i2, " ", f6.3, " ", a1, i2, " ", i2, " ", f5.2, " ", f9.7, " ", i3, " ", f5.2)'
write(*, *) ' 赤経 赤緯 H.P. 視半径'
write(*, *) ' h m s 。 , , ,,'
write(*, sfmt) rh, rm, rs, sign, dh, dm, ds, pos(3), sm, ss
else
sm = int(pos(4))
ss = 60.0d0 * (pos(4) - sm)
sfmt = '(i2, " ", i2, " ", f6.3, " ", a1, i2, " ", i2, " ", f5.2, " ", f9.7, " ", i2, " ", f5.2, " ", f10.7)'
write(*, '("")')
write(*, *) ' 赤経 赤緯 地心距離 視半径 黄道傾斜角'
write(*, *) ' h m s 。 , AU , ,, 。'
write(*, sfmt) rh, rm, rs, sign, dh, dm, ds, pos(3), sm, ss, pos(6)
end if
write(*, *)
stop
end program pcob_mainコンパイル
gfortran pcob_main.f90 pcob_sub.f90実行
./a.out実行結果
2022年5月24日9時00分00秒における視位置を求めてみました。
太陽
DATE AND TIME(JST) ? 20220524,090000
(1: Moon, 2: Sun, 3:Venus, 4: Mars, 5: Jupiter, 6: Saturn) ? 2
2022年 5月 24日 9時 0分 0秒
Sunの視位置
赤経 赤緯 地心距離 視半径 黄道傾斜角
h m s 。 , AU , ,, 。
4 3 6.284 +20 43 25.32 1.0125655 15 49.27 23.4378669月
Moonの視位置
赤経 赤緯 H.P. 視半径
h m s 。 , , ,,
23 23 53.038 - 9 4 32.02 0.9561266 15 37.92金星
Venusの視位置
赤経 赤緯 地心距離 視半径 黄道傾斜角
h m s 。 , AU , ,, 。
1 34 4.675 + 7 41 57.77 1.1649751 0 7.12 23.4378669火星
Marsの視位置
赤経 赤緯 地心距離 視半径 黄道傾斜角
h m s 。 , AU , ,, 。
0 0 3.800 - 1 51 22.71 1.5005738 0 3.12 23.4378669木星
Jupiterの視位置
赤経 赤緯 地心距離 視半径 黄道傾斜角
h m s 。 , AU , ,, 。
0 10 38.521 - 0 5 17.56 5.3888037 0 17.09 23.4378669土星
Saturnの視位置
赤経 赤緯 地心距離 視半径 黄道傾斜角
h m s 。 , AU , ,, 。
21 51 9.696 -14 10 5.64 9.7020342 0 7.61 23.4378669