天体の位置計算 – 光行差による恒星位置のずれ
「恒星位置のずれ-光行差による恒星位置のずれ」74-82頁
光行差は恒星以外にも惑星や彗星、小惑星などの天体の位置計算にも用いられます。光行差は、観測者が運動しているために生ずる見かけ上の恒星位置のずれのことをいいます。光行差には、いくつか種類がありますが、恒星に対しては地球の公転によって起こる年周光行差、惑星や彗星、小惑星などの近距離天体に対しては、天体からの光が届く時間を考慮した惑星光行差が思いられます。ただし、惑星光行差は年周光行差とは別種とも言えるもので、基本原理は違います。
惑星光行差は、たとえば、木星の光が地球に届くまでの時間は、距離によって違いますが、おおよそ30分から40分と言われています。つまり、地球から木星を見た場合、30分から40分前の姿を見ているわけです。木星は公転していますから、30分から40分の間に、木星の位置は変化します。この差を光行差（惑星光行差）といいます。したがって、計算するには、距離を求め、光が届く時間を計算します。その後、計算する日付に光が届く時間を補正して、再度、位置計算をすることで正確な位置計算を行うことができます。
年周光行差は、観測者が運動しているために生じる見かけ上の恒星位置のずれといいましたが、この運動とは地球の公転のことをいいます。公転方向に移動する地球の運動によって、恒星の見かけ上位置がずれる現象です。わかりやすくいうと、雨が降っているとき、歩くと雨が前方から降っている感じというイメジーになるでしょうか。雨が真上から降っているにも斜めになりますね。このずれのことを年周光行差というわけです。
このずれは、厳密に言えば一定ではないものの、光行差定数として決められています。値はおおよそ20″.5となります。
計算式
計算に必要な値は、太陽黄経\(\lambda_{s}\)、平均黄道傾斜角\(\epsilon\)、光行差定数\(\kappa\)、年周視差を考慮した恒星位置(\(L_{4},M_{4},N_{4}\))になります。
最初に光行差定数は、そのままの値では計算に用いることはできませんので、角度に換算したあとラジアン値に換算しておきます。計算式では20″.49552（1976年IAUによる）を使用しています。ちなみに天文年鑑2022では、天文基礎データとして光行差定数は20″.49551となっています。
\[ \kappa = (20.49552 / 3600) / (180 / \pi) \]年周視差を考慮した恒星位置の方向余弦(\(L_{4},M_{4},N_{4}\))を求めます。
\[ \begin{eqnarray} L_{4} &=& \cos\delta\cos\alpha\\ M_{4} &=& \cos\delta\sin\alpha\\ N_{4} &=& \sin\delta \end{eqnarray} \]年周光行差を考慮した恒星位置(\(L_{5},M_{5},N{5}\))は以下の式で求めます。
\[ \begin{pmatrix} L_{5}\\ M_{5}\\ N_{5} \end{pmatrix} = \begin{pmatrix} L_{4}\\ M_{4}\\ N_{4} \end{pmatrix} + \kappa \begin{pmatrix} 1\,-\,M_{4}^{2} & -L_{4}M_{4} & -L_{4}N_{4}\\ -M_{4}L_{4} & 1\,-\,M_{4}^{2} & -M_{4}N_{4}\\ -N_{4}L_{4} & -N_{4}M_{4} & 1\,-\,N_{4}^{2} \end{pmatrix} \begin{pmatrix} \sin\lambda_{s}\\ -\cos\lambda_{s}\cos\epsilon\\ -\cos\lambda_{s}\sin\epsilon \end{pmatrix} \]行列を用いない計算式は以下のとおりになります。
\[ \begin{eqnarray} \Delta L_{1} &=& [(1-L_{4}^{2})\sin\lambda_{s}+L_{4}M_{4}\cos\lambda_{s}\cos\epsilon+L_{4}N_{4}\cos\lambda_{s}\sin\epsilon]\\ \Delta M_{1} &=& [-L_{4}M_{4}\sin\lambda_{s}-(1-M_{4}^{2})\cos\lambda_{s}\cos\epsilon+M_{4}N_{4}\cos\lambda_{s}\sin\epsilon]\\ \Delta N_{1} &=& [-L_{4}N_{4}\sin\lambda_{s}+M_{4}N_{4}\cos\lambda_{s}\cos\epsilon-(1-N_{4}^{2})\cos\lambda_{s}\sin\epsilon]\\ \end{eqnarray} \] \[ \begin{eqnarray} L_{5} &=& L_{4} + \kappa\Delta L_{1}\\ M_{5} &=& M_{4} + \kappa\Delta M_{1}\\ N_{5} &=& N_{4} + \kappa\Delta N_{1} \end{eqnarray} \]Fortranによるコード
▶ サブルーチン AnnualAberration
function AnnualAberration(t0, t1, ra, dc) result(star_position)
! 光行差による恒星位置のずれ
! 「天体の位置計算 増補版」 74-82頁
double precision, parameter :: KP = 20.49552d0
double precision, intent(in) :: t0, t1, ra, dc
double precision :: star_position(3)
double precision :: l, m, n, l5, m5, n5
double precision :: ls, ms, ns
double precision :: ecl, sun_ecl
double precision :: k
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 = sin(sun_ecl)
ms = -cos(sun_ecl) * cos(ecl)
ns = -cos(sun_ecl) * sin(ecl)
! 赤経・赤緯から方向余弦を計算する
l = cos(dc) * cos(ra)
m = cos(dc) * sin(ra)
n = sin(dc)
! 光行差定数をラジアン値に換算
k = deg2rad(KP / 3600d0)
! 年周視差を考慮した恒星位置の計算
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 + k * (matmul(mb, mc))
l5 = me(1, 1)
m5 = me(2, 1)
n5 = me(3, 1)
star_position(1) = l5
star_position(2) = m5
star_position(3) = n5
end function AnnualAberration▶ メインルーチン pcc_aberration.f90
program pcc_aberration
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))
! 視差補正
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))
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)
star_position = AnnualAberration(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_aberrationコンパイル
gfortran pcc_lib.f90 pcc_aberration.f90実行
./a.out実行結果
DATE AND TIME(JST) ? 19781010,203500
1978年 10月 10日 20時 35分 0秒 (JST)
Julian Date = 2443791.98264 UTC
星名 赤経α (視差補正) 赤緯δ 赤経α (光行差補正) 赤緯δ
--------------------------------------------------------------------------------
。 。 。 。
α CMa Sirius 101.052943 -16.688546 101.053630 -16.684970