マイコン宇宙講座-惑星の位置計算Ⅱ

歳差を加えた惑星の位置を求めます。全記事の「惑星の位置計算Ⅰ」では位置の分点が1950.0でしたが、歳差を加えることによって指定した日付の惑星の位置を求めることができます。ただ、「マイコン宇宙講座」にも記載されていますが、本書で使われている軌道要素では水星から火星までの精度は５桁、木星から冥王星まで３桁程度です。これは当時の状況を考えれば仕方のないことです。PC-8001という８ビットパソコンであり、現在のように高精度で計算できるほどのスペックではありませんでした。ですから精度がよくないといっても観測に支障はないと思います。

また、本書で用いられている歳差の計算式は1950年分点です。新しい計算式は「天体の位置計算 増補版」(地人書館 1985 長沢 工著)の「９.新しい計算システム」に掲載されています。ただし、与える時刻引数Tは計算する日付の時刻から元期2000年1月1日正午UTを差し引いて、ユリウス世紀（36525日）単位で表したものになります。計算式は以下のようになります。

T = (計算する日付のユリウス日 – 2451545) / 36525
T2 = T * T
T3 = T2 * T
ζ = 2306″.2181 * T+ 0″.30188 * T2 + 0″.017998 * T3
Ζ = 2306″.2181 * T+ 1″.09468 * T2 + 0″.018203 * T3
θ = 2004″.3109 * T – 0″.42665 * T2 – 0″.041833 * T3

その後は、「天体の位置計算 増補版」54頁 式3-13または式3-14用いてL、M、Nを求めて赤経、赤緯を計算します。ただし、プログラムに組み込む場合は、1度電卓を用いて計算してみて問題ないことを確かめてからにしてください。またサブルーチン precessionの式を上記の式に変更したからといって問題ないというわけではありません。

メインルーチン m37.py

# m37.py
# マイコン宇宙講座
# 3-7 惑星の視位置計算プログラム
import math
import lib


M2PI = 2.0 * math.pi

lib.PL[3] = 'Sun'

x = [[0 for i in range(10)] for j in range(4)]
f = [0, 0, 0, 0]
q = [0, 0, 0, 0]

print()
std = input('DATE AND TIME(JST) ? ')
dy, dt = std.split(',')
dy = float(dy)
dt = float(dt)

jd, yy, mm, dd, hh, ms, ss = lib.mjd(dy, dt)

print()
print('  Date= %5d 年 %2d 月 %2d 日  Time= %2d 時 %2d 分 %2.1f 秒\n' % (yy, mm, dd, hh, ms, ss))
print('         MJD=%12.2f ET\n' % (jd))
print(' Planets        R.A.  (Date)   Decl.         Distance')
print(' ------------------------------------------------------')
print('                h   m            。  ,              AU')

t1 = jd - 33281.92334
t1 = t1 * (2.737909288e-5 + 1.260132857e-17 * t1)
t2 = t1 * t1

ty = yy + (mm - 1) / 12.0 + dd / 365.0
if yy < 0:
    ty = yy - (1 - ((mm - 1) / 12.0 + dd / 365.0))

for pn in range(1, 10):
    e, m, p, n, i, a, rd = lib.mean_elements(pn, t1, t2)
    f, q = lib.eph_const(p, n, i, lib.K)

    ec = e
    mo = M2PI * (m / (M2PI) - int(m / (M2PI)))

    ss, cc, ff = lib.kepler(mo, ec)

    b = a * math.sqrt(1.0 - e * e)
    for n in range(1, 4):
        f[n] = a * f[n]
        q[n] = b * q[n]

    # 惑星座標の計算
    for n in range(1, 4):
        x[n][pn] = ff * f[n] + ss * q[n]

for n in range(1, 4):
    x[n][3] = - (x[n][3])

for pn in range(1, 10):
    if pn != 3:
        for n in range(1, 4):
            x[n][pn] = x[n][pn] + x[n][3]

# 惑星の赤経・赤緯の計算
for pn in range(1, 10):
    cc = x[1][pn]
    ss = x[2][pn]

    tt = lib.quadrant(ss, cc)

    ra = tt
    cc = x[1][pn] / math.cos(ra)
    ss = x[3][pn]

    dc = math.atan(ss / cc)

    r1 = ra
    d1 = dc

    ra1, dc1 = lib.precession(r1, d1, ty, 0, 0, lib.K)

    if ra1 > M2PI:
        ra1 = ra1 - M2PI
    if ra1 < 0:
        ra1 = ra1 + M2PI
    ra = ra1
    dc = dc1

    ra *= lib.K[3] / 15.0
    dc *= lib.K[3]
    ds = math.sqrt(x[1][pn]**2 + x[2][pn]**2 + x[3][pn]**2)

    rh = int(ra)
    rm = 60.0 * (ra - rh)
    dh = int(abs(dc))
    dm = 60.0 * (abs(dc) - dh)
    dh = lib.sgn(dc) * dh
    print(' %-7s      %02d  %05.2f       %+03d  %04.1f       %8.5f' % (lib.PL[pn], rh, rm, dh, dm, ds))

print(' ------------------------------------------------------')
print()

例題 2000年1月1日21時30分00秒(JST)における惑星の位置を求めてみよう。

DATE AND TIME(JST) ? 20000101,213000

  Date=  2000 年  1 月  1 日  Time= 21 時 30 分 0.0 秒

         MJD=    51544.52 ET

 Planets        R.A.  (Date)   Decl.         Distance
 ------------------------------------------------------
                h   m            。  ,              AU
 MERCURY      18  08.52       -24  25.4        1.41559
 VENUS        15  59.73       -18  27.5        1.13782
 Sun          18  45.26       -23  01.9        0.98331
 MARS         22  02.17       -13  10.4        1.84970
 JUPITER      01  36.03       +08  39.1        4.62280
 SATURN       02  34.11       +12  31.8        8.63005
 URANUS       21  12.84       -16  48.7       20.69341
 NEPTUNE      20  18.78       -19  21.3       30.96644
 PLUTO        16  48.58       -11  39.2       31.20937
 ------------------------------------------------------

57