天体の位置計算-視差による恒星位置の変化Ⅰ
「恒星位置のずれ-視差による恒星位置の変化」66-72頁
今回は２回に分けて視差による恒星位置の変化を求めていきます。
地球は太陽から１億５千万Km離れており、その公転軌道を約１年かけて回っています。たとえば、ある日の公転軌道上にある地球は、半年後には太陽を挟んで反対側の位置にあります。これは、ある恒星を地球から見た場合、公転軌道を円とすると、その直径分である３億Kmほど地球の位置にずれがでます。そのため、半年前に見た恒星の見かけ上の位置にわずかながら影響がでてきます。これを視差といいます。
恒星には観測によって求められた視差が存在します。たとえば、うしかい座α星であるArcturusは80.8ミリ秒という視差があります。視差は恒星の精密な位置を求めるのに必要なものであり、観測によって検出された視差をもとに恒星の位置が計算されます、もちろん、恒星の位置を決定するのに必要な要素のひとつであり、これだけでは決定できるものではありません。
計算式
この計算では、太陽の赤経・赤緯もしくは黄経および平均黄道傾斜角が必要になります。
太陽黄経を求める
平均黄道傾斜角はいいとして太陽の黄経を求めるには大変です。「天体の位置計算 増補版」では天体位置表から補間して求めていました。
現在は天体位置表はありません。そこで、海上保安庁海洋情報部の略算式を使用して求めることにします。ただし、太陽の黄経は、係数の桁数からせいぜい小数点4桁までとします。厳密に考えなければ、同書の計算精度から大きな差は出ません。
\[ \begin{array}{r|rrr} \hline i & P_{i} & Q_{i} & R_{i}\\ \hline 1 & 1.1947 & 35999.05 & 267.52\\ 2 & 0.0200 & 71998.10 & 265.10\\ 3 & 0.0200 & 32964.00 & 158.00\\ 4 & 0.0018 & 19.00 & 159.00\\ 5 & 0.0018 & 445267.00 & 208.00\\ 6 & 0.0015 & 45038.00 & 254.00\\ 7 & 0.0013 & 22519.00 & 352.00\\ 8 & 0.0007 & 65929.00 & 45.00\\ 9 & 0.0007 & 3035.00 & 110.00\\ 10 & 0.0007 & 9038.00 & 64.00\\ 11 & 0.0006 & 33718.00 & 316.00\\ 12 & 0.0005 & 155.00 & 118.00\\ 13 & 0.0005 & 2281.00 & 221.00\\ 14 & 0.0004 & 29930.00 & 48.00\\ 15 & 0.0004 & 31557.00 & 161.00\\ 16 & -0.0048 & 35999.00 & 267.52\\ 17 & 0.0048 & 1934.00 & 145.00\\ 18 & -0.0004 & 72002.00 & 111.00\\ \hline \end{array} \]計算は以下の式で求めます。複雑に見えますが、簡単にいうと各項目ごとに求め、最後に総和をとります。表計算ソフトのほうがいいかもしれませんね。
\[ \theta = \displaystyle \sum_{i=1}^{15}P_{i}\cos(Q_{i}T + R_{i}) + \displaystyle \sum_{i=16}^{18}P_{i}T\cos(Q_{i}T + R_{i}) + 36000.7695T + 280.4602 \]
Fortranによるコード
▶ サブルーチン SunEclipticLongitude
function SunEclipticLongitude(t) result(sun_ecliptic)
! 太陽の黄経計算
! 海上保安庁水路部(現 海洋情報部)の略算式による計算
double precision, intent(in) :: t
double precision :: sun_ecliptic(2)
double precision :: xpl(18, 3), rd(8, 3)
double precision :: s, c, r, sun_ecl
integer :: i
! 黄経
xpl( 1, 1:3) = (/ 1.9147d0, 35999.05d0, 267.52d0/)
xpl( 2, 1:3) = (/ 0.0200d0, 71998.1d0, 265.10d0/)
xpl( 3, 1:3) = (/ 0.0020d0, 32964.0d0, 158.00d0/)
xpl( 4, 1:3) = (/ 0.0018d0, 19.0d0, 159.00d0/)
xpl( 5, 1:3) = (/ 0.0018d0,445267.0d0, 208.00d0/)
xpl( 6, 1:3) = (/ 0.0015d0, 45038.0d0, 254.00d0/)
xpl( 7, 1:3) = (/ 0.0013d0, 22519.0d0, 352.00d0/)
xpl( 8, 1:3) = (/ 0.0007d0, 65929.0d0, 45.00d0/)
xpl( 9, 1:3) = (/ 0.0007d0, 3035.0d0, 110.00d0/)
xpl(10, 1:3) = (/ 0.0007d0, 9038.0d0, 64.00d0/)
xpl(11, 1:3) = (/ 0.0006d0, 33718.0d0, 316.00d0/)
xpl(12, 1:3) = (/ 0.0005d0, 155.0d0, 118.00d0/)
xpl(13, 1:3) = (/ 0.0005d0, 2281.0d0, 221.00d0/)
xpl(14, 1:3) = (/ 0.0004d0, 29930.0d0, 48.00d0/)
xpl(15, 1:3) = (/ 0.0004d0, 31557.0d0, 161.00d0/)
xpl(16, 1:3) = (/-0.0048d0, 35999.00d0, 267.52d0/)
xpl(17, 1:3) = (/ 0.0048d0, 1934.0d0, 145.00d0/)
xpl(18, 1:3) = (/-0.0004d0, 72002.0d0, 111.00d0/)
! 地心距離
rd(1, 1:3) = (/ 0.016706d0, 35999.05d0, 177.53d0/)
rd(2, 1:3) = (/ 0.000139d0, 71998.00d0, 175.00d0/)
rd(3, 1:3) = (/ 0.000031d0, 445267.00d0, 298.00d0/)
rd(4, 1:3) = (/ 0.000016d0, 32964.00d0, 68.00d0/)
rd(5, 1:3) = (/ 0.000016d0, 45038.00d0, 164.00d0/)
rd(6, 1:3) = (/ 0.000005d0, 32519.00d0, 233.00d0/)
rd(7, 1:3) = (/ 0.000005d0, 33718.00d0, 226.00d0/)
rd(8, 1:3) = (/-0.000042d0, 35999.00d0, 178.00d0/)
! 黄経の計算
s = 36000.7695d0 * t + 280.4602d0
do i = 1, 15
s = s + xpl(i, 1) * cos(deg2rad(xpl(i, 2) * t + xpl(i, 3)))
end do
c = 0.0d0
do i = 16, 18
c = c + xpl(i, 1) * t * cos(deg2rad(xpl(i, 2) * t + xpl(i, 3)))
end do
sun_ecl = s + c + 0.00569d0
! 動径の計算
r = 1.00014d0
do i = 1, 7
r = r + rd(i, 1) * cos(deg2rad(rd(i, 2) * t + rd(i, 3)))
end do
r = r + rd(8, 1) * t * cos(deg2rad(xpl(8, 2) * t + xpl(8, 3)))
sun_ecliptic(1) = Mla2(sun_ecl)
sun_ecliptic(2) = r
end function SunEclipticLongitude▶ メインルーチン pcc_sunecliptic.f90
program pcc_sunecliptic
use pcclib
implicit none
double precision :: dy, dt, julian
double precision :: year, month, day, hour, minute, second
double precision :: t, ecl, rd
double precision :: sun_ecliptic(2), sun_pos(3), dtime(3)
character(255) :: dfmt
! 日付と時刻の入力
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.375
! ユリウス日から年月日の計算（日の小数）
dtime = Julian2date(julian)
year = dtime(1)
month = dtime(2)
day = dtime(3)
! 太陽の黄経と動径の計算
t = (julian - 2451545d0) / 36525d0
sun_ecliptic = SunEclipticLongitude(t)
ecl = sun_ecliptic(1)
rd = sun_ecliptic(2)
sun_pos = deg2dms(ecl)
! 出力
dfmt = '(" ", i4, " ", i2, " ", f8.5, " ", f14.5, " ", f9.5,"(",i3, "° ", i2, "′ ", f3.1, "″)", f14.5)'
write(*, *)
write(*, *) 'Date / UTC Julian Date 黄経 動径'
write(*, *) '-----------------------------------------------------------------------------------'
write(*, dfmt) int(year), int(month), day, julian, ecl, int(sun_pos(1)), int(sun_pos(2)), sun_pos(3), rd
write(*, *)
stop
end program pcc_sunecliptic
コンパイル
gfortran pcc_lib.f90 pcc_sunecliptic.f90実行
./a.out実行結果
DATE AND TIME(JST) ? 19780610,212000 Date / UTC Julian Date 黄経 動径 ----------------------------------------------------------------------------------- 1978 6 10.51389 2443670.01389 79.26929( 79° 16′ 9.4″) 1.01528
