Python で軌道計算

Python で軌道計算 (Skyfield版)

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

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

from skyfield.api import Topos, load, EarthSatellite

Timescale オブジェクトを読み込み、 t0 に現在時刻を設定します。日時は、世界時で指定します。

ts = load.timescale()
t0 = ts.now()

観測者の位置を指定する Topos オブジェクトを設定します。
観測者の位置は、軌道計算ライブラリの設定の関係から、測心座標（topocentric coordinate）で設定します。
（(earth + Topos) の形では設定しません。）

osaka = Topos('34.6914 N', '135.4917 E')

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

人工衛星の軌道要素のフォーマットとして、NORAD（北アメリカ航空宇宙防衛司令部）が提供する「2行軌道要素形式（TLE:Two-line elements）」がよく用いられます。

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

軌道要素は次のように、EarthSatellite メソッドから設定できます。

line1 = '1 25544U 98067A   20226.06311231  .00000634  00000-0  19556-4 0  9992'
line2 = '2 25544  51.6462  66.9823 0001637  29.7739 108.2756 15.49160058240839'
iss = EarthSatellite(line1, line2, 'ISS (ZARYA)', ts)

または、tle() メソッドのパラメータに NORAD の url を直接指定すれば、ネットワークからTLEのパラメータをそのまま読み込むことができます。
読み込んだパラメータの中から、ISS の要素をインスタンス（iss）に設定します。

satellites = load.tle('http://celestrak.com/NORAD/elements/stations.txt')
iss = satellites['ISS (ZARYA)']

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

あとは、(iss - osaka).at(t0).altaz() とすれば、大阪で現在時刻において、国際宇宙ステーションが見える方向が計算されます。
alt, az, distance に高度、方位角、距離が代入されます。

alt, az, distance = (iss - osaka).at(t0).altaz()

print('高度:{0:.1f} 度'.format(alt.degrees))
print('方位角:{0:.1f} 度'.format(az.degrees))
print('距離:{0:.0f} km'.format(distance.km))

表示結果（例）

高度:-7.3 度
方位角:43.8 度
距離:3313 km

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

応用その１

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

実際に大阪上空を通るパスを計算するには、find_events メソッドを使います。
このメソッドを使えば、指定期間内に国際宇宙ステーションが大阪上空を通るパスのうち、指定仰角以上となるすべてのパスの

（出現時刻、最大高度となる時刻、没時刻）

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

下記は現在時刻から1日後までの間に、最大仰角が10度以上になるパスの時刻を検索した例です。

find_events メソッドにより、events には見え始めの場合 0、最大仰角の場合 1、見え終わりの場合 2 が、t にはそれぞれの時刻が、順番に配列形式で代入されます。

なお、最初に timedelta および timezone メソッドをインポートしておく必要があります。

t0 = ts.now()
t1 = ts.utc(t0.utc_datetime() + timedelta(days=1))
tz = timezone('Asia/Tokyo')

t, events = iss.find_events(osaka, t0, t1, altitude_degrees=10.0)

for ti, event in zip(t, events):
    if event == 0:
        print('見え始め', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
    if event == 1:
        print('最大仰角', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
    if event == 2:
        print('見え終わり', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')

表示結果（例）

見え始め: 2020-08-14 07:43:10 JST
最大仰角: 2020-08-14 07:45:23 JST
見え終わり: 2020-08-14 07:47:36 JST
見え始め: 2020-08-14 09:18:47 JST
・・・

応用その２

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

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

そこで、

太陽高度が地平線下6度以下（observe メソッド）
国際宇宙ステーションには太陽光が当たっている（is_sunlit メソッド）
国際宇宙ステーションの最大仰角が10度以上（find_events メソッド）

という条件を探します。

太陽位置の計算の前には、あらかじめ天体暦を読み込み、太陽、地球のオブジェクト（ここではsun, earth）を作成しておきます。

太陽光が当たっているかどうかは、is_sunlit メソッドで調べることができます。

時刻 t が配列で与えられた場合、それに対応して結果も配列で出力されます。

eph = load('de421.bsp')
sun, earth = eph['sun'], eph['earth']

t, events = iss.find_events(osaka, t0, t1, altitude_degrees=10.0)

sun_alt = (earth + osaka).at(t).observe(sun).apparent().altaz()[0].degrees
sun_lit = iss.at(t).is_sunlit(eph)
 
for ti, event, s_alt, s_lit in zip(t, events, sun_alt, sun_lit):
    if s_alt < -6 and s_lit==True:
        if event == 0:
            print('見え始め:', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
        if event == 1:
            print('最大仰角:', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
        if event == 2:
            print('見え終わり:', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')

表示結果（例）

見え始め: 2020-08-07 19:47:11 JST
最大仰角: 2020-08-07 19:49:55 JST
見え終わり: 2020-08-07 19:52:38 JST

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

from skyfield.api import Topos, load
from datetime import timedelta
from pytz import timezone

# 時刻設定
ts = load.timescale()
t0 = ts.now()
t1 = ts.utc(t0.utc_datetime() + timedelta(days=1))
tz = timezone('Asia/Tokyo')

# 天体暦設定
eph = load('de421.bsp')
sun, earth = eph['sun'], eph['earth']

# 観測地設定
osaka = Topos('34.6914 N', '135.4917 E')

# 国際宇宙ステーションの軌道要素設定
satellites = load.tle('http://celestrak.com/NORAD/elements/stations.txt')
iss = satellites['ISS (ZARYA)']

# 指定時刻間に大阪上空を通るパスの検索
t, events = iss.find_events(osaka, t0, t1, altitude_degrees=10.0)

# 太陽高度、および太陽光が当たっているかの計算
sun_alt = (earth + osaka).at(t).observe(sun).apparent().altaz()[0].degrees
sun_lit = iss.at(t).is_sunlit(eph)

# 国際宇宙ステーションが見える条件の判定 
for ti, event, s_alt, s_lit in zip(t, events, sun_alt, sun_lit):
    if s_alt < -6 and s_lit==True:
        if event == 0:
            print('見え始め:', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
        if event == 1:
            print('最大仰角:', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
        if event == 2:
            print('見え終わり:', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')

江越 航のホームページ