天体の位置計算 – 章動による赤経、赤緯の変化
「恒星位置のずれ-章動による赤経、赤緯の変化」56-66頁
章動は歳差と同じように自転軸が僅かに周期的に変化する運動のことをいいます。変動量は10″にも満たない量であり、あまり赤経や赤緯に大きななずれを与えるものではありません。ですが、グリニッジ視恒星時や視位置といった計算を行う場合は、必ず必要となってきます。計算量が多いのでコンピュータを使った計算向きといえます。「天体の位置計算 増補版」では、63頁から65頁にかけて章動表が掲載されています。長周期章動23項目、短周期章動が46項目となっていて、これらを計算することになります。
ただし、Sofa（Standars of Fundamental Astronomy）で公開されている章動表では「天体の位置計算 増補版」に掲載されているものよりも遥かに多い計算量となります。
計算式
計算は引数をあらかじめ求めておき、章動表を使って計算していきます。引数の計算桁数が多いため、10桁の関数電卓では不足します。関数付きではなくとも16桁の電卓が欲しいところです。電卓がなくても表計算ソフトがあれば下記のように計算することができます。表はダウンロードできるようにしていますので、どういった計算しているのか見てみるのもいいかもしれません。
＃記号が表示されている箇所がありますが、単にセル幅が桁数に合ってないだけです。計算には差し支えありません。
時刻引数
1900年1月0日正午からの経過日数\(d\)および経過日数を36525単位としたものを求めます。
引数
求めた値は360の整数倍を差し引いておきます。
引数を求めたら、章動表の\(\Delta\phi, \Delta\epsilon\)を下記の式で項目ごとに求め、最後に総和をとることで章動値を求めることができます。
\[ \begin{array}{c|c|c} \hline 引\hspace{20pt}数& \Delta\phi& \Delta\epsilon\\ l\hspace{10pt}l^{‘}\hspace{10pt}F\hspace{10pt}D\hspace{10pt}\Omega & \sinの係数 & \cosの係数\\ \hline a\hspace{10pt}b\hspace{10pt}c\hspace{10pt}e\hspace{10pt}f & A & B\\ \hline \end{array} \] \[ \begin{eqnarray} \Delta\phiに対してA\sin(al + bl^{‘} + cF + eD + f\Omega)\\ \Delta\epsilonに対してA\cos(al + bl^{‘} + cF + eD + f\Omega) \end{eqnarray} \]平均黄道傾斜角
\[ \epsilon = 84428.2584\,-\,46.8450t\,-\,0.005904t^{2} + 0.0018108t^{3} \]方向余弦(\(L_{3},M_{3},N_{3}\))
\[ \begin{pmatrix} L_{3}\\ M_{3}\\ N_{3} \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0\\ 0 & \cos(\epsilon + \Delta\epsilon) & -\sin(\epsilon + \Delta\epsilon)\\ 0 & \sin(\epsilon + \Delta\epsilon) & \cos(\epsilon + \Delta\epsilon) \end{pmatrix} \] \[ \hspace{35pt}\times \begin{pmatrix} \cos\Delta\phi & -sin\Delta\phi & 0\\ \sin\Delta\phi & \cos\Delta\phi & 0\\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0\\ 0 & cos\,\epsilon & sin\,\epsilon\\ 0 & -\sin\epsilon & \cos\,\epsilon \end{pmatrix} \begin{pmatrix} L_{2}\\ M_{2}\\ N_{2} \end{pmatrix} \] \[ 引用元:「天体の位置計算 増補版」54頁 式3-17 \]行列を展開した計算式は次のようになります。行列計算をサポートしていない電卓であれば、こちらを使って計算するとよいでしょう。
\[ \begin{eqnarray} L_{3} &=& \cos\Delta\phi L_{2}\,-\,\cos\,\epsilon\sin\Delta\phi M_{2}\,-\,\sin\,\epsilon\sin\Delta\phi N_{2}\\ M_{3} &=& \cos(\epsilon + \Delta\epsilon)\sin\Delta\phi L_{2} + \{\cos(\epsilon + \Delta\epsilon)\cos\,\epsilon\cos\Delta\phi + \sin(\epsilon + \Delta\epsilon)\sin\,\epsilon\}M_{2}\\ & & + \{\cos(\epsilon + \Delta\epsilon)\sin\,\epsilon\cos\Delta\phi\,-\,\sin(\epsilon + \Delta\epsilon)\cos\,\epsilon\}N_{2}\\ N_{3} &=& \sin(\epsilon + \Delta\epsilon)\sin\phi L_{2} + \{\sin(\epsilon + \Delta\epsilon)\cos\,\epsilon\cos\Delta\phi\,-\,\cos(\epsilon + \Delta\epsilon)\sin\,\epsilon \}M_{2}\\ & & + \{\sin(\epsilon + \Delta\epsilon)\sin\,\epsilon\cos\Delta\phi + \cos(\epsilon + \Delta\epsilon)\cos\,\epsilon\}N_{2} \end{eqnarray} \] \[ 引用元:「天体の位置計算 増補版」54頁 式3-18 \]Fortranによるコード
サブルーチンは従来から２つに分けました。黄経の章動⊿φと黄道傾斜角の章動⊿εを求めるiauNut50と、章動計算するNutationです。こうしておくと、iauNut50からグリニッジ視恒星時の計算に必要な分点差Eqを求めることができます。
▶ サブルーチン iauNut50
function iauNut50(dt, t) result(eps)
! 章動の計算
! 「天体の位置計算 増補版」63-65頁
double precision, intent(in) :: dt, t
double precision :: xpl(69, 9)
double precision :: arg, s, c, dt2, dt3
double precision :: dp ! 黄経の章動
double precision :: de ! 黄道傾斜角の章動
double precision :: el ! 月の平均近点角
double precision :: elp ! 太陽の平均近点角
double precision :: f ! 昇交点から測った月の平均黄経
double precision :: d ! 太陽と月の平均離角
double precision :: om ! 月の平均昇交点黄経
double precision :: eps(2)
integer :: i
xpl( 1, 1:9) = (/ 0d0, 0d0, 0d0, 0d0, 1d0, -17.2327d0, -0.01737d0, 9.2100d0, 0.00091d0/)
xpl( 2, 1:9) = (/ 0d0, 0d0, 2d0, -2d0, 2d0, -1.2729d0, -0.00013d0, 0.5522d0, -0.00029d0/)
xpl( 3, 1:9) = (/ 0d0, 0d0, 0d0, 0d0, 2d0, 0.2088d0, 0.00002d0, -0.0904d0, 0.00004d0/)
xpl( 4, 1:9) = (/ 0d0, 1d0, 0d0, 0d0, 0d0, 0.1261d0, -0.00031d0, 0.0000d0, 0.00000d0/)
xpl( 5, 1:9) = (/ 0d0, 1d0, 2d0, -2d0, 2d0, -0.0497d0, 0.00012d0, 0.0216d0, -0.00006d0/)
xpl( 6, 1:9) = (/ 0d0, -1d0, 2d0, -2d0, 2d0, 0.0214d0, -0.00005d0, -0.0093d0, 0.00003d0/)
xpl( 7, 1:9) = (/ 0d0, 0d0, 2d0, -2d0, 1d0, 0.0124d0, 0.00001d0, -0.0066d0, 0.00000d0/)
xpl( 8, 1:9) = (/-2d0, 0d0, 2d0, 0d0, 1d0, 0.0045d0, 0.00000d0, -0.0024d0, 0.00000d0/)
xpl( 9, 1:9) = (/ 2d0, 0d0, 0d0, -2d0, 0d0, 0.0045d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(10, 1:9) = (/ 0d0, 0d0, 2d0, -2d0, 0d0, -0.0021d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(11, 1:9) = (/ 0d0, 2d0, 0d0, 0d0, 0d0, 0.0016d0, -0.00001d0, 0.0000d0, 0.00000d0/)
xpl(12, 1:9) = (/ 0d0, 2d0, 2d0, -2d0, 2d0, -0.0015d0, 0.00001d0, 0.0007d0, 0.00000d0/)
xpl(13, 1:9) = (/ 0d0, 1d0, 0d0, 0d0, 1d0, -0.0015d0, 0.00000d0, 0.0008d0, 0.00000d0/)
xpl(14, 1:9) = (/ 0d0, -1d0, 0d0, 0d0, 1d0, -0.0010d0, 0.00000d0, 0.0005d0, 0.00000d0/)
xpl(15, 1:9) = (/ 2d0, 0d0, -2d0, 0d0, 0d0, 0.0010d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(16, 1:9) = (/-2d0, 0d0, 0d0, 2d0, 1d0, -0.0005d0, 0.00000d0, 0.0003d0, 0.00000d0/)
xpl(17, 1:9) = (/ 0d0, -1d0, 2d0, -2d0, 1d0, -0.0005d0, 0.00000d0, 0.0003d0, 0.00000d0/)
xpl(18, 1:9) = (/ 0d0, -2d0, 2d0, -2d0, 1d0, -0.0004d0, 0.00000d0, 0.0002d0, 0.00000d0/)
xpl(19, 1:9) = (/ 2d0, 0d0, 0d0, -2d0, 1d0, 0.0004d0, 0.00000d0, -0.0002d0, 0.00000d0/)
xpl(20, 1:9) = (/ 0d0, 1d0, 2d0, -2d0, 1d0, 0.0003d0, 0.00000d0, -0.0002d0, 0.00000d0/)
xpl(21, 1:9) = (/-2d0, 0d0, 2d0, 0d0, 2d0, -0.0003d0, 0.00000d0, 0.0002d0, 0.00000d0/)
xpl(22, 1:9) = (/ 1d0, 0d0, 0d0, -1d0, 0d0, -0.0003d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(23, 1:9) = (/ 1d0, -1d0, 0d0, -1d0, 0d0, -0.0002d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(24, 1:9) = (/ 0d0, 0d0, 2d0, 0d0, 2d0, -0.2037d0, -0.00002d0, 0.0884d0, -0.00005d0/)
xpl(25, 1:9) = (/ 1d0, 0d0, 0d0, 0d0, 0d0, 0.0675d0, 0.00001d0, 0.0000d0, 0.00000d0/)
xpl(26, 1:9) = (/ 0d0, 0d0, 2d0, 0d0, 1d0, -0.0342d0, -0.00004d0, 0.0183d0, 0.00000d0/)
xpl(27, 1:9) = (/ 1d0, 0d0, 2d0, 0d0, 2d0, -0.0261d0, 0.00000d0, 0.0113d0, -0.00001d0/)
xpl(28, 1:9) = (/ 1d0, 0d0, 0d0, -2d0, 0d0, -0.0149d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(29, 1:9) = (/-1d0, 0d0, 2d0, 0d0, 2d0, 0.0114d0, 0.00000d0, -0.0050d0, 0.00000d0/)
xpl(30, 1:9) = (/ 0d0, 0d0, 0d0, 2d0, 0d0, 0.0060d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(31, 1:9) = (/ 1d0, 0d0, 0d0, 0d0, 1d0, 0.0058d0, 0.00000d0, -0.0031d0, 0.00000d0/)
xpl(32, 1:9) = (/-1d0, 0d0, 0d0, 0d0, 1d0, -0.0057d0, 0.00000d0, 0.0030d0, 0.00000d0/)
xpl(33, 1:9) = (/-1d0, 0d0, 2d0, 2d0, 2d0, -0.0052d0, 0.00000d0, 0.0022d0, 0.00000d0/)
xpl(34, 1:9) = (/ 1d0, 0d0, 2d0, 0d0, 1d0, -0.0044d0, 0.00000d0, 0.0023d0, 0.00000d0/)
xpl(35, 1:9) = (/ 0d0, 0d0, 2d0, 2d0, 2d0, -0.0032d0, 0.00000d0, 0.0014d0, 0.00000d0/)
xpl(36, 1:9) = (/ 2d0, 0d0, 0d0, 0d0, 0d0, 0.0028d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(37, 1:9) = (/ 1d0, 0d0, 2d0, -2d0, 2d0, 0.0026d0, 0.00000d0, -0.0011d0, 0.00000d0/)
xpl(38, 1:9) = (/ 2d0, 0d0, 2d0, 0d0, 2d0, -0.0026d0, 0.00000d0, 0.0011d0, 0.00000d0/)
xpl(39, 1:9) = (/ 0d0, 0d0, 2d0, 0d0, 0d0, 0.0025d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(40, 1:9) = (/-1d0, 0d0, 2d0, 0d0, 1d0, 0.0019d0, 0.00000d0, -0.0010d0, 0.00000d0/)
xpl(41, 1:9) = (/-1d0, 0d0, 0d0, 2d0, 1d0, 0.0014d0, 0.00000d0, -0.0007d0, 0.00000d0/)
xpl(42, 1:9) = (/ 1d0, 0d0, 0d0, -2d0, 1d0, -0.0013d0, 0.00000d0, 0.0007d0, 0.00000d0/)
xpl(43, 1:9) = (/-1d0, 0d0, 2d0, 2d0, 1d0, -0.0009d0, 0.00000d0, 0.0005d0, 0.00000d0/)
xpl(44, 1:9) = (/ 1d0, 1d0, 0d0, -2d0, 0d0, -0.0007d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(45, 1:9) = (/ 0d0, 1d0, 2d0, 0d0, 2d0, 0.0007d0, 0.00000d0, -0.0003d0, 0.00000d0/)
xpl(46, 1:9) = (/ 1d0, 0d0, 0d0, 2d0, 0d0, 0.0006d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(47, 1:9) = (/ 0d0, 0d0, 0d0, 2d0, 1d0, -0.0006d0, 0.00000d0, 0.0003d0, 0.00000d0/)
xpl(48, 1:9) = (/ 0d0, -1d0, 2d0, 0d0, 2d0, -0.0006d0, 0.00000d0, 0.0003d0, 0.00000d0/)
xpl(49, 1:9) = (/ 1d0, 0d0, 2d0, 2d0, 2d0, -0.0006d0, 0.00000d0, 0.0003d0, 0.00000d0/)
xpl(50, 1:9) = (/ 2d0, 0d0, 2d0, -2d0, 2d0, 0.0006d0, 0.00000d0, -0.0002d0, 0.00000d0/)
xpl(51, 1:9) = (/ 0d0, 0d0, 0d0, -2d0, 1d0, -0.0005d0, 0.00000d0, 0.0003d0, 0.00000d0/)
xpl(52, 1:9) = (/ 0d0, 0d0, 2d0, 2d0, 1d0, -0.0005d0, 0.00000d0, 0.0003d0, 0.00000d0/)
xpl(53, 1:9) = (/ 1d0, 0d0, 2d0, -2d0, 1d0, 0.0005d0, 0.00000d0, -0.0003d0, 0.00000d0/)
xpl(54, 1:9) = (/ 0d0, 0d0, 0d0, 1d0, 0d0, -0.0004d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(55, 1:9) = (/ 0d0, 1d0, 0d0, -2d0, 0d0, -0.0004d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(56, 1:9) = (/ 1d0, -1d0, 0d0, 0d0, 0d0, 0.0004d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(57, 1:9) = (/ 1d0, 0d0, -2d0, 0d0, 0d0, 0.0004d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(58, 1:9) = (/ 2d0, 0d0, 2d0, 0d0, 1d0, -0.0004d0, 0.00000d0, 0.0002d0, 0.00000d0/)
xpl(59, 1:9) = (/ 1d0, 0d0, 2d0, 0d0, 0d0, 0.0003d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(60, 1:9) = (/ 1d0, 1d0, 0d0, 0d0, 0d0, -0.0003d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(61, 1:9) = (/ 1d0, -1d0, 2d0, 0d0, 2d0, -0.0003d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(62, 1:9) = (/-2d0, 0d0, 0d0, 0d0, 1d0, -0.0002d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(63, 1:9) = (/-1d0, 0d0, 2d0, -2d0, 1d0, -0.0002d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(64, 1:9) = (/ 2d0, 0d0, 0d0, 0d0, 1d0, 0.0002d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(65, 1:9) = (/-1d0, -1d0, 2d0, 2d0, 2d0, -0.0002d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(66, 1:9) = (/ 0d0, -1d0, 2d0, 2d0, 2d0, -0.0002d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(67, 1:9) = (/ 1d0, 0d0, 0d0, 0d0, 2d0, -0.0002d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(68, 1:9) = (/ 1d0, 1d0, 2d0, 0d0, 2d0, 0.0002d0, 0.00000d0, 0.0000d0, 0.00000d0/)
xpl(69, 1:9) = (/ 3d0, 0d0, 2d0, 0d0, 2d0, -0.0002d0, 0.00000d0, 0.0000d0, 0.00000d0/)
dt2 = dt * dt
dt3 = dt2 * dt
el = 296.104608d0 + 13.0649924465d0 * dt + 6.890d-12 * dt2 + 2.95d-20 * dt3
elp = 358.475833d0 + 0.9856002669d0 * dt - 0.112d-12 * dt2 - 6.80d-20 * dt3
f = 11.250889d0 + 13.2293504490d0 * dt - 2.407d-12 * dt2 - 0.70d-20 * dt3
d = 350.737486d0 + 12.1907491914d0 * dt - 1.076d-12 * dt2 + 3.90d-20 * dt3
om = 259.183275d0 - 0.0529539222d0 * dt + 1.557d-12 * dt2 + 4.60d-20 * dt3
! 360°の整数倍を差し引く
el = Mla2(el)
elp = Mla2(elp)
f = Mla2(f)
d = Mla2(d)
om = Mla2(om)
! Standards of Fundamental Astronomy
! http:!www.iausofa.org/2019_0722_C/CompleteList.html
! from nut80.c
dp = 0d0
de = 0d0
do i = 1, 69
arg = el * xpl(i, 1) + elp * xpl(i, 2) + f * xpl(i, 3) + d * xpl(i, 4) + om * xpl(i, 5)
s = (xpl(i, 6) + (xpl(i, 7) * t))
c = (xpl(i, 8) + (xpl(i, 9) * t))
if (s /= 0d0) then
dp = dp + s * sin(deg2rad(arg))
end if
if (c /= 0d0) then
de = de + c * cos(deg2rad(arg))
end if
end do
! end of nut80.c
eps(1) = dp
eps(2) = de
end function iauNut50 ▶ サブルーチン Nutation
function Nutation(dt, t, ra, dc) result(star_position)
! 章動による赤経・赤緯の変化
! 「天体の位置計算 増補版」 58-66頁
double precision, intent(in) :: dt
double precision, intent(in) :: t
double precision, intent(in) :: ra, dc
double precision :: ecl, dp, de
double precision:: star_position(3), eps(2)
double precision :: l, m, n, l3, m3, n3
double precision :: ma(3, 3), mb(3, 3), mc(3, 3), md(3, 1), me(3, 1)
! 章動の計算
eps = iauNut50(dt, t)
dp = deg2rad(eps(1) / 3600d0)
de = deg2rad(eps(2) / 3600d0)
! 平均黄道傾斜角の計算(1950)
ecl = deg2rad(iauMOE50(t))
! 赤経・赤緯から方向余弦を計算する
l = cos(dc) * cos(ra)
m = cos(dc) * sin(ra)
n = sin(dc)
! 方向余弦の計算
ma(1, 1:3) = (/1d0, 0d0, 0d0/)
ma(2, 1:3) = (/0d0, cos(ecl + de), -sin(ecl + de)/)
ma(3, 1:3) = (/0d0, sin(ecl + de), cos(ecl + de)/)
mb(1, 1:3) = (/cos(dp), -sin(dp), 0d0/)
mb(2, 1:3) = (/sin(dp), cos(dp), 0d0/)
mb(3, 1:3) = (/0d0, 0d0, 1d0/)
mc(1, 1:3) = (/1d0, 0d0, 0d0/)
mc(2, 1:3) = (/0d0, cos(ecl), sin(ecl)/)
mc(3, 1:3) = (/0d0, -sin(ecl), cos(ecl)/)
md(1, 1) = l
md(2, 1) = m
md(3, 1) = n
me = matmul(matmul(ma, mb), matmul(mc, md))
l3 = me(1, 1)
m3 = me(2, 1)
n3 = me(3, 1)
star_position(1) = l3
star_position(2) = m3
star_position(3) = n3
end function Nutation▶ メインルーチン pcc_nutation.f90
program pcc_nutation
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, t0, t1
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))
! 歳差による赤経・赤緯の変化
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));
ra_tmp = ra_deg
dc_tmp = dc_deg
! 章動補正
t0 = julian - 2415020d0
t1 = t0 / 36525d0
ra_rad = deg2rad(ra_deg)
dc_rad = deg2rad(dc_deg)
star_position = Nutation(t0, t1, 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_nutationコンパイル
gfortran pcc_lib.f90 pcc_nutation.f90実行
./a.out実行結果
DATE AND TIME(JST) ? 19781010,203500
1978年 10月 10日 20時 35分 0秒 (JST)
Julian Date = 2443791.98264 UTC
星名 赤経α (歳差補正) 赤緯δ 赤経α (章動補正) 赤緯δ
--------------------------------------------------------------------------------
。 。 。 。
α CMa Sirius 101.053151 -16.686164 101.052835 -16.688532