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

惑星の位置を表す赤経・赤緯（1950.0分点）を求めます。惑星の位置計算は、小惑星や彗星の位置計算にも役に立てることができますので計算方法を覚えておくとよいでしょう。ただし、彗星は基本的に楕円軌道となりますので、計算方法が少し異なります。なお、このプログラムで求められる位置は、1950年分点ですので現在に照らし合わせると合いません。あくまでも書籍の出版時のデータに基づいた計算になります。ただ、計算方法は今も昔も変わりません。もし厳密に計算するのであれば、惑星光行差、歳差、章動、最新の軌道要素などを含めて計算する必要があります。

サブルーチン eph_const

※サブルーチンはlib.pyに記述してください。

# EPH CONST
def eph_const(pe, nd, ic, K):
    f = [0, 0, 0, 0]
    q = [0, 0, 0, 0]

    r1 = math.sin(pe)
    r2 = math.sin(nd)
    r3 = math.sin(ic)
    r4 = math.cos(pe)
    r5 = math.cos(nd)
    r6 = math.cos(ic)

    f[1] = r5 * r4 - r1 * r6 * r2
    q[1] = -r1 * r5 - r4 * r6 * r2
    r7 = r6 * r5 * K[4]
    r8 = r2 * K[4]
    r9 = r3 * K[5]
    f[2] = r1 * r7 + r4 * r8 - r1 * r9
    q[2] = r4 * r7 - r1 * r8 - r4 * r9
    r7 = r3 * K[4]
    r8 = r6 * r5 * K[5]
    r9 = r2 * K[5]
    f[3] = r1 * r7 + r1 * r8 + r4 * r9
    q[3] = r4 * r7 + r4 * r8 - r1 * r9

    return f, q

メインルーチン m35.py

# m35.py
# マイコン宇宙講座
# 3-5 惑星の位置計算プログラム
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.  (1950.0)  Decl.         Distance')
print(' ------------------------------------------------------')
print('                h   m            。  ,              AU')

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

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)
    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年１月１日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.  (1950.0)  Decl.         Distance
 ------------------------------------------------------
                h   m            。  ,              AU
 MERCURY      18  05.45       -24  25.9        1.41559
 VENUS        15  56.85       -18  19.1        1.13782
 Sun          18  42.23       -23  05.1        0.98331
 MARS         21  59.48       -13  24.9        1.84970
 JUPITER      01  33.40       +08  23.8        4.62280
 SATURN       02  31.39       +12  18.6        8.63005
 URANUS       21  10.05       -17  01.0       20.69341
 NEPTUNE      20  15.89       -19  30.7       30.96644
 PLUTO        16  45.80       -11  34.0       31.20937
 ------------------------------------------------------

「マイコン宇宙講座」の出力結果と較べると、天王星から冥王星までのDistance値に最大±0.00004の差が出ています。N-BASICとPythonの計算精度の差によるものかなと思っています。

57