マイコン宇宙講座-離心近点離角計算

ケプラーの運動方程式を解くことで、天体の位置計算に必要な離心近点離角を求めます。ケプラーの運動方程式を解くには平均近点離角と離心率が必要となります。これら２つは軌道要素として与えられていますので、サブルーチンkeplerの引数と与えれば、離心近点離角を求めることができます。ここからさまざまな計算過程を経て天体の位置を計算することができるわけです。

観測結果から軌道要素を決定する方法については「マイコン天文学Ⅰ」に書かれていますので、興味のある方は参考にしてみるとよいでしょう。昔はこれらの手順を対数表や三角関数表などを用いて手で計算していたんですね。中には最小二乗法をそろばんで計算していたという方もおられたようです。その後、電卓から関数電卓、コンピュータと移り変わり、もはやコンピュータなしでは計算が成り立たないというほどになりました。また天文ソフトウェアの普及によりアマチュアでも簡単にこれらのソフトウェアを利用して、さまざまな天文現象をディスプレイ上で観たり、天体望遠鏡と接続して天体写真の撮影や観測に役立てられるようになり、プロと変わらない環境を手に入れられるようになりました。

サブルーチン kepler

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

6行目のprint文は、これ以降は必要ないのでコメントにしておいてください。また、return文で返している変数a1も不要になりますので、同じくこれ以降は、return ss, cc, ff #, a1に変更しておいてください。これらは収れんの様子を確認するためのものです。

# KEPLER
def kepler(mo, ec):
    a1 = 0.0
    while True:
        a2 = ec * math.sin(a1 + mo)
        print('    %10.7f     %10.7f' % (a1, a2))
        if abs(a2 - a1) > 1.0e-5:
            a1 = a2
        else:
            a1 = a2 + mo
            break

    ss = math.sin(a1)
    cc = math.cos(a1)
    ff = cc - ec

    return ss, cc, ff, a1

メインルーチン m28.py

# m28.py
# マイコン宇宙講座
# 2-8 離心近点離角計算プログラム
import lib


print()
mo = float(input('Mo ? ')) / lib.K[3]
ec = float(input('E ? '))

print()
print('INPUT Mo = %9.5f' % (mo * lib.K[3]))
print(('INPUT e  =   %9.7f' % (ec)))
print()
# print('          A1          A2')

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

print()
print('E = %9.5f' % (a1 * lib.K[3]))
print()

例題 平均近点離角(Mo)と離心率(e)を与えて離心近点離角を求めてみよう。

Mo ? 10
E ? 0.5

INPUT Mo =  10.00000
INPUT e  =   0.5000000

     0.0000000      0.0868241
     0.0868241      0.1291959
     0.1291959      0.1495402
     0.1495402      0.1592151
     0.1592151      0.1637933
     0.1637933      0.1659543
     0.1659543      0.1669732
     0.1669732      0.1674533
     0.1674533      0.1676795
     0.1676795      0.1677860
     0.1677860      0.1678362
     0.1678362      0.1678598
     0.1678598      0.1678710
     0.1678710      0.1678762

E =  19.61860

57