天体の位置計算 – 固有運動による恒星位置のずれ
「恒星位置ずれ-固有運動による恒星位置のずれ」41-49頁
恒星は永劫未来、同じ位置にあるのではなく、ほんのわずかづつ位置が変化していきます。位置の変化にはさまざまな要因があります。そのため、恒星の位置を正しく求めるには、これらの要素を計算して求めなければいけません。また、恒星の位置を表す赤経・赤緯はゆっくりと変化するために、星表に掲載されている位置データをそのまま使うことはできません。
恒星の位置のずれを引き起こす要因としては、次のようなものがあります。
- 固有運動
- 歳差や章動
- 極運動
- 視差
- 光行差
- 大気差
これらの要因は、ほんのわずかな量なので一般的にはあまり大きな問題を引き起こすものではありませんが、だからといって全く無視をしていると、実視界の狭い望遠鏡だと視野から外れてしまうということが起こります。今回は、固有運動による恒星位置のずれを求めてみましょう。
固有運動は、Halley彗星の発見で有名なEdmond・Halleyが、1718年に古代ヒッパルコスなどの星表と比較してアークトゥルスやシリウスの位置が変わっていることに気づいて発見しました。この発見は固有運動発見の糸口となりました。固有運動は恒星自身が運動してその位置を変えていくことをいい、理論的に導き出せるものでなく観測によって決められます。そのために、長い観測時間を経て測定して位置を求め、その運動量が決定されます。
星表データには赤経方向の固有運動、赤緯方向の固有運動というものが記載されています。星表によっては、赤経方向をミリ秒/年、赤緯方向をミリ秒/年で表したり、s/100年、″/100年で表したりします。ちなみに、61 Cygの現在の固有運動は、SIMBAD Astronomical Databaseによると、赤経方向の固有運動 4164.209ミリ秒/年、赤緯方向の固有運動 3249.614ミリ秒/年となっています。
近似式
星表に掲載されている赤経および赤緯を\(\alpha_{0}\)、赤緯を\(\delta_{0}\)とすると、\(t\)年後の赤経\(\alpha_{1}\)および赤緯\(\delta_{1}\)は以下の近似式を使って求めることができます。固有運動の大きい恒星を除けば近似式でも十分です。
計算式
近似式を使用しない場合は、次の手順で求めます。
１）1950年のベッセル年初から計算日までの経過年数\(t\)を次式で計算しておきます
\[ \begin{eqnarray} MJD &=& JD – 2400000.5\\ t &=& MJD\,-\,33281.923\\ &=& t(0.00002737909288 + 1.25013^{-17}t)\\ \end{eqnarray} \] \[ 引用元:「マイコン宇宙講座」31頁 式2.3.8 \]２）全固有運動速度\(\mu_{0}\)および固有運動方向\(\psi\)を求めます
\[ \begin{eqnarray} \mu_{0} = \sqrt{\mu_{\alpha}^{2}\cos^{2}\delta_{0}+\mu_{\delta}^{2}}\\ \tan\psi_{0} = \frac{\mu_{\alpha}\cos\delta_{0}}{\mu_{\delta}}\\ \end{eqnarray} \] \[ 引用元:「天体の位置計算 増補版」44頁 式3-2 式3-3 \]３）\(t\)年間における恒星移動の角距離\(\sigma\)を求めます
\[ \sigma = \mu_{0}t\left(1\,-\,\frac{\pi}{1.496e^{8}}R_{a}t\right) \] \[ 引用元:「天体の位置計算 増補版」45頁 式3-6 \]４）固有運動で移動した位置(\(\alpha_{1}, \delta_{1}\))の方向余弦(\(L_{1}, M_{1}, N_{1}\))を求めます
\[ \begin{pmatrix} L_{1}\\ M_{1}\\ N_{1} \end{pmatrix} = \begin{pmatrix} \cos\alpha_{0} & -\sin\alpha_{0} & 0\\ \sin\alpha_{0} & \cos\alpha_{0} & 0\\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \cos\delta_{0} & 0 & -\sin\delta_{0}\\ 0 & 1 & 0\\ \sin\delta_{0} & 0 & \cos\delta_{0} \end{pmatrix} \begin{pmatrix} 1 & 0 & 0\\ 0 & \cos\psi_{0} & \sin\psi_{0}\\ 0 & -\sin\psi_{0} & \cos\psi_{0} \end{pmatrix} \begin{pmatrix} \cos\sigma\\ 0\\ \sin\sigma \end{pmatrix} \] \[ 引用元:「天体の位置計算 増補版」46頁 式3-8 \]行列を展開した計算式は次のようになります。行列計算をサポートしていない電卓であれば、こちらを使って計算するとよいでしょう。
\[ \begin{eqnarray} L_{1} &=& \cos\alpha_{0}\cos\delta_{0}\cos\sigma\,-\,(\sin\alpha_{0}\sin\psi_{0} + \cos\alpha_{0}\sin\delta_{0}\cos\psi_{0})\sin\sigma\\ M_{1} &=& \sin\alpha_{0}\cos\delta_{0}\cos\sigma + (\cos\alpha_{0}\sin\psi_{0}\,-\,\sin\alpha_{0}\sin\delta_{0}\cos\psi_{0})\sin\sigma\\ N_{1} &=& \sin\delta_{0}\cos\sigma + \cos\delta_{0}\cos\psi_{0}\sin\sigma \end{eqnarray} \] \[ 引用元:「天体の位置計算 増補版」46頁 式3-9 \]Fortranによるコード
▶サブルーチン ProperMotion
function ProperMotion(t, ra, dc, ra_proper, dc_proper, ap, rv) result(star_position)
! 固有運動補正による恒星位置のずれ
! 「天体の位置計算 増補版」41-49頁
double precision, intent(in) :: t
double precision, intent(in) :: ra, dc, ra_proper, dc_proper, ap, rv
double precision :: star_position(3)
double precision :: l1, m1, n1
double precision :: mu, phi, sigma, ss, cc
double precision :: ma(3, 3), mb(3, 3), mc(3, 3), md(3, 1), me(3, 1)
! 全固有運動速度の計算
mu = sqrt((ra_proper * cos(dc))**2 + dc_proper**2)
! 固有運動の向きを計算
ss = ra_proper * cos(dc)
cc = dc_proper
phi = Quadrant(ss, cc)
! 1年間における恒星の移動の角距離を計算
sigma = mu * t * (1d0 - (ap / 1.496d8) * rv * (t * (365.2422d0 * 86400.0d0)))
sigma = deg2rad(sigma)
! 方向余弦の計算
ma(1, 1:3) = (/cos(ra), -sin(ra), 0d0/)
ma(2, 1:3) = (/sin(ra), cos(ra), 0d0/)
ma(3, 1:3) = (/0d0, 0d0, 1d0/)
mb(1, 1:3) = (/cos(dc), 0d0, -sin(dc)/)
mb(2, 1:3) = (/0d0, 1d0, 0d0/)
mb(3, 1:3) = (/sin(dc), 0d0, cos(dc)/)
mc(1, 1:3) = (/1d0, 0d0, 0d0/)
mc(2, 1:3) = (/0d0, cos(phi), sin(phi)/)
mc(3, 1:3) = (/0d0, -sin(phi), cos(phi)/)
md(1, 1) = cos(sigma)
md(2, 1) = 0d0
md(3, 1) = sin(sigma)
me = matmul(matmul(ma, mb), matmul(mc, md))
l1 = me(1, 1)
m1 = me(2, 1)
n1 = me(3, 1)
star_position(1) = l1
star_position(2) = m1
star_position(3) = n1
end function ProperMotion▶メインルーチン pcc_proper_motion.f90
program pcc_proper_motion
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 :: 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))
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(*, *) '星名 赤経α (1950) 赤緯δ 赤経α (固有運動補正) 赤緯δ'
write(*, *) '--------------------------------------------------------------------------------'
write(*, *) ' 。 。 。 。'
write(*, dfmt) SNM, hms2deg(SRA(1), SRA(2), SRA(3)), sgr, abs(dms2deg(SDC(1), SDC(2),SDC(3))), ra_deg, sg, abs(dc_deg)
write(*, *)
stop
end program pcc_proper_motion
コンパイル
gfortran pcc_lib.f90 pcc_proper_motion.f90実行
./a.out実行結果
DATE AND TIME(JST) ? 19781010,203500
1978年 10月 10日 20時 35分 0秒 (JST)
Julian Date = 2443791.98264 UTC
星名 赤経α (1950) 赤緯δ 赤経α (固有運動補正) 赤緯δ
----------------------------------------------------------------------------------
。 。 。 。
α CMa Sirius 100.736308 -16.646211 100.731763 -16.655894