天体の位置計算 – 歳差による赤経、赤緯の変化
「恒星位置のずれ-恒星による赤経、赤緯の変化」49-55頁
地球の地軸は垂直に立っているのではなく、23.4度ほど傾いて自転しています。この地軸は常に一定方向を指しているのではなく、公転面に対して半径23.4度の円を描いて、26000年ほどの周期で回転しています。この動きを歳差運動といいます。歳差運動が起こるのは太陽や月、惑星といった天体が引力によって地軸を引き起こそうとすることが原因です。さらに地軸はごくわずかですが、振動することがあります。これを章動といいます。歳差や章動はコマに例えられることが多く、コマを一定傾けて回転させると同じようなことが起こります。
天体の位置計算では歳差や章動を考慮して計算します。というのは、たとえば惑星の軌道要素の分点がJ2000とします。現在は2023年ですので、今見えている惑星の位置を求めるには23年間分の歳差を考慮しなければなりません。今回は歳差のみを考慮した天体の位置計算を行っていきます。
計算式
最初に歳差定数（\(\zeta,z,\theta\)）を計算し、計算された歳差定数から方向余弦を求めますが、ニューカムが理論と観測によって以下の式を導き出しています
t0およびtはベッセル世紀（36524.22日）を単位として表したものです。t0は1900.0ベッセル年初（2415020.314）から元期のベッセル年初までの経過日数をベッセル世紀単位として表したもの。下記の式は、ｔは元期1950.0となるベッセル年初（2433292.42346）から計算日までの経過日数をベッセル世紀単位で表したものになります。t0とtは以下の式で求めることができます。
\[ \begin{eqnarray} t_{0} &=& (JDE\,-\,2415020.314) / 36524.22\\ t &=& (JD\,-\,2433292.42346) / 36524.22 \end{eqnarray} \]歳差常数(元期 1950.0)
これらによって、式3-12は元期を1950.0として計算するものとしてt0を0.5としています。上記の計算式では厳密には0.5とならず0.4999…となりますが、精度を考えても微小なので0.5として考えても問題ありません。1950.0からの経過日数はサブルーチン BesselianYearで求めています。計算では、以下の式を使用して歳差による恒星位置の変化を求めます。
方向余弦(\(L_{2}, M_{2}, N_{2}\))
\[ \begin{pmatrix} L_{2}\\ M_{2}\\ N_{2} \end{pmatrix} = \begin{pmatrix} -\sin z & -\cos z & 0\\ \cos z & -\sin z & 0\\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0\\ 0 & \cos\theta & \sin\theta\\ 0 & -\sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} \sin\zeta_{0} & \cos\zeta_{0} & 0\\ -\cos\zeta_{0} & \sin\zeta_{0} & 0\\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} L_{1}\\ M_{1}\\ N_{1} \end{pmatrix} \] \[ 引用元:「天体の位置計算 増補版」54頁 式3-13 \]行列を展開した計算式は次のようになります。行列計算をサポートしていない電卓であれば、こちらを使って計算するとよいでしょう。
\[ \begin{eqnarray} L_{2} &=& (\cos\zeta_{0}\cos z\cos\theta\,-\,\sin\zeta_{0}\sin z)L_{1} + (-\sin\zeta_{0}\cos z\cos\theta\,-\,\cos\zeta_{0}\sin z)M_{1} + (-\cos z\sin\theta)N_{1}\\ M_{2} &=& (\cos\zeta_{0}\sin z\cos\theta + \sin\zeta_{0}\cos z)L_{1} + (-\sin\zeta_{0}\sin z\cos\theta + \cos\zeta_{0}\cos z)M_{1} + (-\sin z\sin\theta)N_{1}\\ N_{2} &=& \cos\zeta_{0}\sin\theta L_{1} + (-\sin\zeta_{0}\sin\theta)M_{1} + \cos\theta N_{1} \end{eqnarray} \] \[ 引用元:「天体の位置計算 増補版」54頁 式3-14 \]Fortranによるコード
▶ サブルーチン Precession
サブルーチンに渡す引数を固有運動で移動した赤経・赤緯を渡します。
function Precession(t, ra, dc) result(star_position)
! 歳差の計算
! 天体の位置計算 増補版 49-55頁
double precision, intent(in) :: t
double precision, intent(in) :: ra, dc
double precision :: star_position(3)
double precision :: x1, x2, x3
double precision :: l, m, n, l2, m2, n2
double precision :: ma(3, 3), mb(3, 3), mc(3, 3), md(3, 1), me(3, 1)
! 赤経・赤緯から方向余弦を計算する
l = cos(dc) * cos(ra)
m = cos(dc) * sin(ra)
n = sin(dc)
! 歳差定数の計算(EPOCH=1950)
x1 = 2304.9480d0 * t + 0.302d0 * t**2 + 0.018d0 * t**3
x2 = 2304.9480d0 * t + 1.093d0 * t**2 + 0.019d0 * t**3
x3 = 2004.2555d0 * t - 0.426d0 * t**2 - 0.042d0 * t**3
x1 = deg2rad(x1 / 3600d0)
x2 = deg2rad(x2 / 3600d0)
x3 = deg2rad(x3 / 3600d0)
! 方向余弦の計算
ma(1, 1:3) = (/-sin(x2), -cos(x2), 0d0/)
ma(2, 1:3) = (/ cos(x2), -sin(x2), 0d0/)
ma(3, 1:3) = (/0d0, 0d0, 1d0/)
mb(1, 1:3) = (/1d0, 0d0, 0d0/)
mb(2, 1:3) = (/0d0, cos(x3), sin(x3)/)
mb(3, 1:3) = (/0d0, -sin(x3), cos(x3)/)
mc(1, 1:3) = (/ sin(x1), cos(x1), 0d0/)
mc(2, 1:3) = (/-cos(x1), sin(x1), 0d0/)
mc(3, 1:3) = (/0d0, 0d0, 1d0/)
md(1, 1) = l
md(2, 1) = m
md(3, 1) = n
me = matmul(matmul(ma, mb), matmul(mc, md))
l2 = me(1, 1)
m2 = me(2, 1)
n2 = me(3, 1)
star_position(1) = l2
star_position(2) = m2
star_position(3) = n2
end function Precession▶ メインルーチン pcc_precession.f90
program pcc_precession
use pcclib
implicit none
! 恒星データ
character(16), parameter :: SNM = 'α CMa Sirius ' ! 恒星名
double precision, parameter :: SRA(1:3) = (/6d0, 42d0, 56.714d0/) ! 赤経
double precision, parameter :: SDC(1:3) = (/-16d0, 38d0, 46.36d0/) ! 赤緯
double precision, parameter :: SRA_PROPER_MOTION = -0.03791d0 ! 赤経方向の固有運動量
double precision, parameter :: SDC_PROPER_MOTION = -1.2114d0 ! 赤緯方向の固有運動量
double precision, parameter :: SRV = -7.6d0 ! 視線速度
double precision, parameter :: SAP = 0.377d0 ! 年周視差
!検証用恒星データ
! character(16), parameter :: SNM = '61 Cyg '; ! 恒星名
! double precision, parameter :: SRA(1:3) = (/21d0, 4d0, 39.935d0/) ! 赤経
! double precision, parameter :: SDC(1:3) = (/38d0, 29d0, 59.10d0/) ! 赤緯
! double precision, parameter :: SRA_PROPER_MOTION = 0.35227d0 ! 赤経方向の固有運動量
! double precision, parameter :: SDC_PROPER_MOTION = 3.1847d0 ! 赤緯方向の固有運動量
! double precision, parameter :: SRV = -64.3d0 ! 視線速度
! double precision, parameter :: SAP = 0.296d0 ! 年周視差
double precision :: dy, dt, julian, t
double precision :: year, month, day, hour, minute, second
double precision :: ra_proper_motion, dc_proper_motion
double precision :: ra_deg, dc_deg, ra_rad, dc_rad, ap
double precision :: ra_tmp, dc_tmp
double precision :: l, m, n, ss, cc, tt
double precision :: star_position(3)
character(255) :: sfmt, dfmt
character(1) :: sg, sgr
! 日付と時刻の入力
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.375d0
! 1959年ベッセル年初からの経過年数
t = BesselianYear(julian)
t = t * 100d0
! 固有運動補正
ra_rad = deg2rad(hms2deg(SRA(1), SRA(2), SRA(3)))
dc_rad = deg2rad(dms2deg(SDC(1), SDC(2),SDC(3)))
ra_proper_motion = SRA_PROPER_MOTION / 3600d0 * 15d0;
dc_proper_motion = SDC_PROPER_MOTION / 3600d0;
ap = deg2rad(SAP / 3600d0)
star_position = ProperMotion(t, ra_rad, dc_rad, ra_proper_motion, dc_proper_motion, ap, SRV)
l = star_position(1)
m = star_position(2)
n = star_position(3)
! 固有運動の方向余弦から赤経・赤緯の計算
ss = m
cc = l
tt = Quadrant(ss, cc)
ra_deg = rad2deg(tt)
dc_deg = rad2deg(asin(n))
ra_tmp = ra_deg
dc_tmp = dc_deg
! 歳差補正
t = t / 100d0
ra_rad = deg2rad(ra_deg)
dc_rad = deg2rad(dc_deg)
star_position = Precession(t, ra_rad, dc_rad);
l = star_position(1)
m = star_position(2)
n = star_position(3)
! 歳差の方向余弦から赤経・赤緯の計算
ss = m
cc = l
tt = quadrant(ss, cc);
ra_deg = rad2deg(tt);
dc_deg = rad2deg(asin(n));
sgr = '+'
if (SDC(1) < 0.0d0) then
sgr = '-'
end if
sg = '+'
if (dc_deg < 0.0d0) then
sg = '-'
end if
! 出力
! 角度の譚雨を表す°記号は出力に合わせて調整すること
sfmt = '(i4, "年 ", i2, "月 ", i2, "日 ", i2, "時 ", i2, "分 ", i2, "秒 (JST)")'
dfmt = '(a16, " ", f10.6, " ", a1, f9.6, " ", f10.6, " ", a1, f9.6)'
write(*, *)
write(*, sfmt) int(year), int(month), int(day), int(hour), int(minute), int(second)
write(*, '("Julian Date = ", f14.5, " UTC")') julian
write(*, *)
write(*, *) '星名 赤経α (固有運動補正) 赤緯δ 赤経α (歳差補正) 赤緯δ'
write(*, *) '--------------------------------------------------------------------------------'
write(*, *) ' 。 。 。 。'
write(*, dfmt) SNM, ra_tmp, sgr, abs(dc_tmp), ra_deg, sg, abs(dc_deg)
write(*, *)
stop
end program pcc_precession
コンパイル
gfortran pcc_lib.f90 pcc_precession.f90実行
./a.out実行結果
DATE AND TIME(JST) ? 19781010,203500
1978年 10月 10日 20時 35分 0秒 (JST)
Julian Date = 2443791.98264 UTC
星名 赤経α (固有運動補正) 赤緯δ 赤経α (歳差補正) 赤緯δ
--------------------------------------------------------------------------------
。 。 。 。
α CMa Sirius 100.731763 -16.655894 101.053151 -16.686164