天体の位置計算-視差による恒星位置の変化Ⅱ
「恒星位置のずれ-視差による恒星位置の変化」66-72頁
前回は、視差を考慮した恒星の位置を求めるために、必要な太陽の黄経を求めました。今回は実際に視差によるずれを補正する計算式を使って位置を求めてみます。計算に必要なデータは、恒星の位置（章動を含めた赤経・赤緯）、年周視差、太陽黄経、平均黄道傾斜角です。計算は恒星の位置に視差の補正値を加えるというものになります。
計算式
今までと同じように方向余弦を使って求めます。最初に恒星の位置（赤経・赤緯）を、「天体の位置計算 増補版」34頁の2−18式を使って方向余弦(\(L_{3},M_{3},N_{3}\))に変換します。
次に補正値\(A(L_{s},M_{s},N_{s}\))を求めます。
\[ A(L,M,N) = \begin{pmatrix} 1\,-\,M_{3}^{2} & -L_{3}M_{3} & -L_{3}N_{3}\\ -M_{3}L_{3} & 1\,-\,M_{3}^{2} & -M_{3}N_{3}\\ -N_{3}L_{3} & -N_{3}M_{3} & 1\,-\,N_{3}^{2} \end{pmatrix} \]赤経・赤緯の方向余弦\(L_{3},M{3},N_{3}\)に年周視差\(\pi\)と補正値\(A(L,M,N)\)、太陽黄経の方向余弦\(L_{s},M_{s},N_{s}\)を除算したものを加算して求めます。
\[ \begin{eqnarray} \begin{pmatrix} L_{4}\\ M_{4}\\ N_{4} \end{pmatrix} = \begin{pmatrix} L_{3}\\ M_{3}\\ N_{3} \end{pmatrix} + \pi A(L,M,N) \begin{pmatrix} L_{s}\\ M_{s}\\ N_{s} \end{pmatrix} \end{eqnarray} \]行列を使わないのであれば次式で同じように求めることができます。
\[ \begin{eqnarray} L_{4} &=& L_{3}\,\, + \pi[(1-L_{3}^{2})L_{s}-L_{3}M_{3}M_{s}-L_{3}N_{3}N_{s}]\\ M_{4} &=& M_{3} + \pi[-L_{3}M_{3}L_{s}+(1-M_{3}^{2})M_{s}-M_{3}N_{3}N_{s}]\\ N_{4} &=& N_{3}\, + \pi[-L_{3}N_{3}L_{s}-M_{3}N_{3}M_{s}+(1-N_{3}^{2})N_{s}] \end{eqnarray} \]
Fortranによるコード
▶ サブルーチン AnnualParallax
function AnnualParallax(t0, t1, ra, dc, ap) result(star_position)
! 視差による恒星位置のずれ
! 「天体の位置計算 増補版」66-72頁
double precision, intent(in) :: t0, t1, ra, dc, ap
double precision :: star_position(3)
double precision :: l, m, n, l4, m4, n4
double precision :: ls, ms, ns
double precision :: ecl, sun_ecl
double precision :: sun_ecliptic(2), ma(3, 1), mb(3, 3), mc(3, 1), me(3, 1)
! 平均黄道傾斜角
ecl = deg2rad(iauMOE50(t1))
! 太陽の黄経の方向余弦
sun_ecliptic = SunEclipticLongitude(t0)
sun_ecl = deg2rad(sun_ecliptic(1))
ls = cos(sun_ecl)
ms = sin(sun_ecl) * cos(ecl)
ns = sin(sun_ecl) * sin(ecl)
! 赤経・赤緯から方向余弦を計算する
l = cos(dc) * cos(ra)
m = cos(dc) * sin(ra)
n = sin(dc)
! 年周視差を考慮した恒星位置の計算
ma(1, 1) = l
ma(2, 1) = m
ma(3, 1) = n
mb(1, 1:3) = (/1d0-l**2, -l*m, -l*n/)
mb(2, 1:3) = (/-m*l, 1d0-m**2, -m*n/)
mb(3, 1:3) = (/-n*l, -n*m, 1d0-n**2/)
mc(1, 1) = ls
mc(2, 1) = ms
mc(3, 1) = ns
me = ma + ap * (matmul(mb, mc))
l4 = me(1, 1)
m4 = me(2, 1)
n4 = me(3, 1)
star_position(1) = l4
star_position(2) = m4
star_position(3) = n4
end function AnnualParallax▶ メインルーチン pcc_parallax.f90
時刻引数がt0とt1になっているのは、サブルーチン内で太陽の黄経（2000年分点）と平均黄道傾斜角（1950年分点）のためです。サブルーチンを呼び出す前に、どちらかを求めておき、値を引数で渡すということにした方がよいかもしれませんが、ここでは、そのままにしました。
program pcc_parallax
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));
! 章動補正
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))
ra_tmp = ra_deg
dc_tmp = dc_deg
! 視差補正
t0 = (julian - 2451545d0) / 36525d0
t1 = (julian - 2415020d0) / 36525d0
ra_rad = deg2rad(ra_deg)
dc_rad = deg2rad(dc_deg)
ap = deg2rad(SAP / 3600d0)
star_position = AnnualParallax(t0, t1, ra_rad, dc_rad, ap)
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_parallaxコンパイル
gfortran pcc_parallax.f90 pcc_library.f90実行
./a.out実行結果
DATE AND TIME(JST) ? 19781010,203500
1978年 10月 10日 20時 35分 0秒 (JST)
Julian Date = 2443791.98264 UTC
星名 赤経α (章動補正) 赤緯δ 赤経α (視差補正) 赤緯δ
--------------------------------------------------------------------------------
。 。 。 。
α CMa Sirius 101.052835 -16.688532 101.052943 -16.688546