Python で軌道計算

Python で軌道計算 (PyEphem版)

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

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

import ephem
import datetime

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

osaka = ephem.Observer()
osaka.lat = '34.6914'
osaka.lon = '135.4917'
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度の方向に国際宇宙ステーションが見える、ということになります。

応用その１

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

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

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

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

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

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

print('見え始め時刻:', ephem.localtime(rise_t), '方位角: {0:.1f}'.format(np.rad2deg(az_rise)))
print('最大仰角時刻:', ephem.localtime(max_t), '最大仰角: {0:.1f}'.format(np.rad2deg(alt_max)))
print('見え終わり時刻:', ephem.localtime(set_t), '方位角: {0:.1f}'.format(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

応用その２

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

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

そこで、

太陽高度が地平線下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('見え始め時刻:', ephem.localtime(rise_t).strftime('%Y-%m-%d %H:%M:%S'), '方位角: {0:.1f}'.format(np.rad2deg(az_rise)))
print('最大仰角時刻:', ephem.localtime(max_t).strftime('%Y-%m-%d %H:%M:%S'), '最大仰角: {0:.1f}'.format(np.rad2deg(alt_max)))
print('見え終わり時刻:', ephem.localtime(set_t).strftime('%Y-%m-%d %H:%M:%S'), '方位角: {0:.1f}'.format(np.rad2deg(az_set)))

表示結果（例）

見え始め時刻: 2019-09-13 04:46:43 方位角: 318.6
最大仰角時刻: 2019-09-13 04:52:11 最大仰角: 37.7
見え終わり時刻: 2019-09-13 04:57:40 方位角: 117.1

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

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

import ephem
import datetime
import numpy as np

# 観測地設定
osaka = ephem.Observer()
osaka.lat = '34.6914'
osaka.lon = '135.4917'
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('見え始め時刻:', ephem.localtime(rise_t).strftime('%Y-%m-%d %H:%M:%S'), '方位角: {0:.1f}'.format(np.rad2deg(az_rise)))
print('最大仰角時刻:', ephem.localtime(max_t).strftime('%Y-%m-%d %H:%M:%S'), '最大仰角: {0:.1f}'.format(np.rad2deg(alt_max)))
print('見え終わり時刻:', ephem.localtime(set_t).strftime('%Y-%m-%d %H:%M:%S'), '方位角: {0:.1f}'.format(np.rad2deg(az_set)))

江越 航のホームページ