大阪市立科学館学芸員 江越航のホームページ  大阪市立科学館 学芸員

 江越 航のホームページ

Python で軌道計算

PyEphem を用いると、人工衛星の軌道計算を行うことができます。
ここでは、国際宇宙ステーション(ISS)「きぼう」が見える位置を計算してみます。

必要モジュールをインポートします。

import ephem
import datetime

観測者のインスタンス(ここでは osaka)を作成します。
作成した観測者 osaka に、現在日時(世界時)を設定します。

osaka = ephem.Observer()
osaka.lat = '34.6883'
osaka.lon = '135.4933'
osaka.date = datetime.datetime.utcnow()

国際宇宙ステーションのインスタンス(ここでは iss)を作成し、軌道要素を設定します。

人工衛星の軌道要素のフォーマットとして、NORAD(北アメリカ航空宇宙防衛司令部)が提供する「2行軌道要素形式(TLE:Two-line elements)」がよく用いられます。
PyEphem では readtle メソッドを使うことで、TLEのパラメータをそのまま読み込むことができます。

なお、各衛星のTLEは、以下のNORADのホームページに掲載されています。
http://www.celestrak.com/NORAD/elements/stations.txt

iss = ephem.readtle('ISS',
                    '1 25544U 98067A   19251.51298924  .00000370  00000-0  14350-4 0  9995',
                    '2 25544  51.6453 307.9614 0007973  18.9504 129.2301 15.50436909188185')

以上で準備はおしまいです。

あとは、iss.compute(osaka) とすれば、大阪で現在時刻において、国際宇宙ステーションが見える方向が計算されます。
iss.az で方位角、iss.alt で高度の値を知ることができます。
角度はラジアンなので、度の単位にするには、180/π を掛けます。

iss.compute(osaka)
print(iss.az * (180/3.1415))
print(iss.alt * (180/3.1415))

表示結果(例)

271.91030965875177
-31.207935191179825

上記結果は、計算をした時刻において、方位角272度(北を基準に時計回り、つまりほぼ真西)、地平線下31度の方向に国際宇宙ステーションが見える、ということになります。

応用その1

上記方法では、地平線下で見えない場合も計算されます。

実際に大阪上空を通るパスを計算するには、next_pass メソッドを使います。
このメソッドを使えば、次に国際宇宙ステーションが上空を通るパスの

(出現時刻、出現方位角、最大高度となる時刻、最大高度、没時刻、没方向)

を一度に得ることができます。

なお下記では日本標準時で表示するために ephem.localtime メソッドを、ラジアンを度に変換するために np.rad2deg メソッドを使っています。

rise_t, az_rise, max_t, alt_max, set_t, az_set = osaka.next_pass(iss)

print(u"見え始め時刻: %s 方位角: %.1f" % (ephem.localtime(rise_t), np.rad2deg(az_rise)))
print(u"最大迎角時刻: %s 最大迎角: %.1f" % (ephem.localtime(max_t), np.rad2deg(alt_max)))
print(u"見え終わり時刻: %s 方位角: %.1f" % (ephem.localtime(set_t), np.rad2deg(az_set)))

表示結果(例)

見え始め時刻: 2019-09-12 07:13:00.000003 方位角: 285.3
最大迎角時刻: 2019-09-12 07:17:15.000003 最大迎角: 8.4
見え終わり時刻: 2019-09-12 07:21:30.000003 方位角: 184.8

応用その2

上記方法ではまだ、昼間で実際には見えないパスも表示されます。

肉眼で国際宇宙ステーションを見るためには、日の入り後ある程度時間が経った時、もしくは日の出前で、辺りが暗く、かつ上空にある国際宇宙ステーションには太陽光が当たっている条件を探す必要があります。
また、国際宇宙ステーションが地平線近くでも見ることができません。

そこで、

太陽高度が地平線下6度以下(np.rad2deg(sun.alt) < -6)
国際宇宙ステーションには太陽光が当たっている(iss.eclipsed is False)
国際宇宙ステーションの最大迎角が10度以上(np.rad2deg(alt_max) > 10)

という条件を探します。

条件に合わない場合は、90分後を初期値として、再度計算します。
(国際宇宙ステーションは約90分で地球を一周するため)

また、太陽高度の計算の前には、あらかじめ太陽のインスタンス(ここではsun)も準備しておきます。

sun = ephem.Sun()

while True:
    rise_t, az_rise, max_t, alt_max, set_t, az_set = osaka.next_pass(iss)

    osaka.date = max_t
    sun.compute(osaka)
    iss.compute(osaka)
        
    if iss.eclipsed is False and np.rad2deg(sun.alt) < -6 and np.rad2deg(alt_max) > 10:
        break
    else:
        osaka.date = rise_t + (1/(24*60)*90)

print(u"見え始め時刻: %s 方位角: %.1f" % (ephem.localtime(rise_t), np.rad2deg(az_rise)))
print(u"最大迎角時刻: %s 最大迎角: %.1f" % (ephem.localtime(max_t), np.rad2deg(alt_max)))
print(u"見え終わり時刻: %s 方位角: %.1f" % (ephem.localtime(set_t), np.rad2deg(az_set)))

表示結果(例)

見え始め時刻: 2019-09-13 04:46:43.000004 方位角: 318.6
最大迎角時刻: 2019-09-13 04:52:11.000004 最大迎角: 37.7
見え終わり時刻: 2019-09-13 04:57:40.000004 方位角: 117.1

なお、上記条件で検索すると、最大迎角に達する前に地球の影に入った場合、あるいは最大迎角に達した後に地球の影から出てきた場合は、見逃してしまうことに注意が必要です。

全体は次のようになります。

import ephem
import datetime
import numpy as np

# 観測地設定
osaka = ephem.Observer()
osaka.lat = '34.6883'
osaka.lon = '135.4933'
osaka.date = datetime.datetime.utcnow()

# 国際宇宙ステーションの軌道要素
iss = ephem.readtle('ISS',
                    '1 25544U 98067A   19251.51298924  .00000370  00000-0  14350-4 0  9995',
                    '2 25544  51.6453 307.9614 0007973  18.9504 129.2301 15.50436909188185')

# 太陽
sun = ephem.Sun()

# 国際宇宙ステーションが見える条件が見つかるまで繰り返す
while True:
	# 次回、大阪上空を通るパス
    rise_t, az_rise, max_t, alt_max, set_t, az_set = osaka.next_pass(iss)

	# 最大高度になる時刻での、国際宇宙ステーションと太陽の位置を計算
    osaka.date = max_t
    sun.compute(osaka)
    iss.compute(osaka)
        
	# 国際宇宙ステーションが見える条件の判定
    if iss.eclipsed is False and np.rad2deg(sun.alt) < -6 and np.rad2deg(alt_max) > 10:
        break
    else:
        osaka.date = rise_t + (1/(24*60)*90)

print(u"見え始め時刻: %s 方位角: %.1f" % (ephem.localtime(rise_t), np.rad2deg(az_rise)))
print(u"最大迎角時刻: %s 最大迎角: %.1f" % (ephem.localtime(max_t), np.rad2deg(alt_max)))
print(u"見え終わり時刻: %s 方位角: %.1f" % (ephem.localtime(set_t), np.rad2deg(az_set)))