ロケットの瞬間落下点の計算（OpenVerne） #Python

ロケットが飛行中に変な方向に飛んでいかないようにすることは当然ですが非常に重要です。

ロケットの飛翔中は上空を高速で飛翔しています。観測ロケットで最大でマッハ４とか、軌道投入すると最大でマッハ20とか。そのような高速飛翔物体が万一の際に、突然推力が無くなって落下する位置というのは安全確保上必要なデータです。しかし、飛翔中のロケットの現在位置・速度だけ見ていても落下予測地点を瞬時に理解するのは難しいです。
緊急停止などの重要な判断にはそのような瞬間落下点（IIP:Instantaneous Impact Point)が必要になることから、地上で監視している人はロケット飛翔中にIIPをモニタしています。

文科省の文章によるとIIPは下記の通りです。

IIP：IIPはロケットの推力をある時点で瞬時に停止した時のロケットの落下予想点を表すもので、飛翔経路の正常・異常を瞬時に判断するのに適している情報である。

重要な割に、どういう計算しているのか国の文章には書いていなかったり実装例が無いので、pythonでの実装例を作って公開しました。
https://github.com/istellartech/OpenVerne

実行例1

example_00.py

# -*- coding: utf-8 -*-
""" Simplest example """

from OpenVerne import IIP
import numpy as np

posLLH_ = np.array([35.0, 140.0, 100])
velNED_ = np.array([10, 0, 0])

_IIP = IIP(posLLH_, velNED_)
print(_IIP)

とOpenVerne.pyと同じフォルダ内で実行してみると、北緯３５度、東経１４０度、高度１００mの地点から速度北方向に１０m/sに物体が進んでいるときのIIPは落下までの時間5.18秒で52.2m移動することが出力結果からわかります。試しにGoogleEarthで表示すると下記の通り。

実行例2

入力値に対してIIPが計算できるということは、逆に2点間をケプラー運動的に動いたときに必要な初期速度が計算できます。
つまり、任意の地点から地球表面の任意の点に”物体”をお届けするのに必要な速度が計算できます。例えば東京から大阪にお届けするための速度をOpenVerneを用いて計算すると・・・

example_03.py

# -*- coding: utf-8 -*-
""" This example is for setting the falling point to an arbitrary position on the Earth.
"""

from OpenVerne import IIP
import numpy as np
from scipy.optimize import minimize

if __name__ == '__main__':

    def f(velNED_, pos_init, pos_goal):
        _IIP = IIP(pos_init, velNED_)
        lat1 = _IIP.posLLH_IIP_deg[0]
        lon1 = _IIP.posLLH_IIP_deg[1]
        lat2 = pos_goal[0]
        lon2 = pos_goal[1]
        return (lat1 - lat2) ** 2 + (lon1 - lon2) ** 2

    pos_Tokyo = np.array([35.67288, 139.74336, 10000])
    pos_Osaka = np.array([34.68963, 135.53010, 0])

    velNED_init = np.array([100, -1000, -1000])

    ans = minimize(f, velNED_init, args=(pos_Tokyo, pos_Osaka))
    print(ans.message)

    _IIP = IIP(pos_Tokyo, ans.x)
    a = 340.29  # 地上での音速 [m/s]
    print("必要な速度\n北方向 = %.1f [m/s], 東方向 = %.1f [m/s], 上方向 = %.1f [m/s]" % (ans.x[0], ans.x[1], -ans.x[2]))
    print("北方向 = %.1f [km/h], 東方向 = %.1f [km/h], 上方向 = %.1f [km/h]" % (ans.x[0]*3.6, ans.x[1]*3.6, -ans.x[2]*3.6))
    print("マッハ数 : %.2f " % (np.sqrt((ans.x[0]/a)**2 + (ans.x[1]/a)**2 + (ans.x[2]/a)**2)))
    # print(_IIP)

必要な速度
北方向 = -378.9 [m/s], 東方向 = -1365.1 [m/s], 上方向 = 1314.9 [m/s]
北方向 = -1363.9 [km/h], 東方向 = -4914.2 [km/h], 上方向 = 4733.8 [km/h]
マッハ数 : 5.68

と出力されます。ただ単に届けばいい軌道を解いているだけなので、最小エネルギーの軌道を作っているわけではないですが、かなりそれでもマッハ5程度で投げ上げてあげないと届かないことがわかります。怪力の人がハンマー投げしても難しいですね。

このように瞬間的な速度があっても遠くまで物体を届けるのが難しいので、大砲やレールガンなどの大きな初期速度を得るような加速器を用いた宇宙輸送システムは実用化されません。

ちなみに、もっと遠くの地点を同じ方法で解こうとすると（例えば東京からWashingtonDCなど地球の裏側に違いところ）到達する軌道が無数に作れてしまうので、かなり最適化問題を工夫してあげないと解が安定しません。

このような重力のみの力でケプラー運動する物体の2地点間境界問題はLambert問題と呼ばれていて、これはこれで既に効率的な解法がたくさん公開されているものですので、普通はLambert問題として解きます。上記例はOpeｎVerneの使い方としての参考例です。

由来

OpenVerneの名前はジュール・ヴェルヌの月世界旅行から来ています。
月世界旅行の中では大砲を使って月に向かいます。初速だけでどういう軌道や落下地点になるかの計算をするスクリプトということでジュール・ヴェルヌらしさを感じました。

参考

IIPの計算内容は、高校物理でやる物体投げ上げの放物線運動の延長です。IIPの計算自体はいくつかの計算方法があり、OpenVerneでの実装の特徴は下記の通りです。韓国のKARI（韓国のNASA的組織）ではIIPに関する新しい研究がされてたりするので、KARIの人の論文を参考にしています。

物体の初期位置（緯度経度高度）と速度（NED座標系）を入力するとIIPの緯度経度出力
落下物体はケプラーの法則に従う地球中心を焦点とする楕円運動を行う
地球は楕円体近似したモデル
空気抵抗等の重力以外の力は無視
落下地点の高度は楕円体の高度ゼロ地点
ケプラー運動ではあるが、地球の自転は考慮している
慣性座標系での運動であるので、コリオリの力や遠心力は考慮しなくて問題無し