ロケットの飛行解析ソフト（OpenTsiolkovsky）をPythonに移植してみた

## 1. はじめに
こんにちは。
この度、インターステラテクノロジズ株式会社様がオープンソースで公開されているロケットの飛行解析ツール「OpenTsiolkovsky」をPythonに移植してみましたので、詳細について解説していきたいと思います。
ツールは以下GitHubの方に公開しておりますのでご覧ください。

https://github.com/shamoji-os/OpenTsiolkovsky-Python-

また、出典は以下となりますのであわせてご確認くださいませ。

https://github.com/istellartech/OpenTsiolkovsky

https://qiita.com/ina111/items/a56ec9221a2ecaa11ae0

本記事では、Qiitaらしく(？)主にコーディングやプログラムの設計、今回の移植にあたって工夫したところといった、技術的な内容に焦点を当ててまとめていきたいと思います。
そもそも「飛行解析って何？」「どんなことが出来るの？」といった疑問をお持ちの方には、手前味噌で大変恐縮ですが、以下ブログの方に今回の概要をまとめた記事を別途作成しておりますので、ご一読されてから本記事をお読みいただけると理解が深まるかと思います。

https://shamoji-space.hatenablog.com/entry/2023/11/25/102224


## 2. 移植の目的
今回、Python移植を行った主な目的は勉強目的です。
ロケットの飛行解析は、ロケットの飛行に関わる様々なデータ・物理法則・座標系等を往ったり来たりしながら行う複雑な解析になります。それらについて、ソースコードレベルの詳細な設計方法等を学びたく、今回オープンソースで開発されている飛行解析用ソフト「OpenTsiolkovsky」を移植するという手法を採りました。
「OpenTsiolkovsky」はC++でコーディングされたメインのロジック部と、Pythonスクリプト等から構成されています。今回はC++でコーディングされた部分をPythonでコーディングし直すことが主な目標となります。

なお、折角勉強目的で移植するということもあり、今回の移植にあたっては「ソースコードの可読性」を最も重要視しました。
特にプログラム全体のデザインを整理しておくことは後の内容の理解にも大きく関わってくると思いましたので、まずはコードをベタ移植しながらプログラムの挙動等を確認しつつ、そこから全体の流れを見てクラス設計を考え、ある程度まとまってきたタイミングでリファクタリングする、といったやり方でなるべく見通しの良いプログラムにしようと試みました。

## 3. プログラムの設計

最終的な本プログラムのクラス設計図は以下のようになりました。

![OpenTsiolkovsky(Python)_クラス図.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490279/75ab1da5-5498-09ac-9a27-75770b254478.png)

全体的にオブジェクト指向っぽく考えました。
ロケット、周囲の環境(大気、風、重力)、プログラム全体の進行管理の役割を担うクラス等を作成し、それらを相互に作用させることで解析を進める形としています。
ロケットの部分については、ロケットの各段(第1段、第2段、第3段など)を表現する「RocketStage」クラスと、それらを束ねる「Rocket」クラスで表現されるものとして設計しました。各段の状況(パラメータ)は「RocketState」クラスから生成したオブジェクトにそれぞれ保存することとしています。ここでRocketStateクラスはメンバ変数のみからなる、C++でいうところの謂わば構造体のような役割を担っており、解析結果の出力、途中経過の保存、クラス間でのデータ授受等の際には、本クラスから生成したオブジェクトが活用できるように意識して設計しました。

プログラム全体の処理の流れとしては大まかに以下のような挙動をしています。その処理の役割を担っているpythonスクリプトの名前もあわせて入れておきましたので、適宜上のクラス図と照らし合わせながらご覧頂ければと思います。

1. プログラム実行： OpenTsiolkovskyPy.py
2. 入力ファイル(json5)の読み取り： define.py
3. シミュレーション開始準備の指示： Simulator.py → Rocket.py
4. 2の入力情報を基にロケットを構成、各段を初期化： Rocket.py, RocketStage.py
5. 飛行解析(以下のループ処理)を開始： Simulator.py
    5-1. 現在の解析対象のステージを確認： Rocket.py
    5-2. 位置情報を計算して更新： RocketStage.py
    5-3. 高度が0mを下回った場合は解析中止： Rocket.py
    5-4. 更新した位置情報から周辺の環境情報(重力、大気、風など)を計算： Rocket.py ↔ Environment.py
    5-5. その他の状況(ロケットの飛行状況、推進系の状態など)を更新する： RocketStage.py
    5-6. 解析状況の確認、特にイベント(ステージの分離等、非連続な事象)の処理等を行う： Rocket.py
    5-7. 解析途中の結果を保存： Rocket.py - observer() → Logger.py
6. 終了処理。解析完了。

運動方程式はルンゲクッタ法を用いて解いています。
Simulator.pyのメインループ内でcalculator.py中のrunge_kutta関数が呼ばれていますが、ここに引数として渡しているRocketオブジェクトがrunge_kutta関数から呼び出されると、Rocketオブジェクト内に定義されている__call__関数が呼び出されて都度の計算結果を返すようになっています。
この辺りは、scipyのode等便利なライブラリがあったものの、解析途中でログの保存だったりと言った別の処理を割り込ませたかったため、ライブラリの使用は断念しています。この辺は移植元のC++コードではライブラリの方で柔軟に対応できていたようで、Python移植にあたり少し詰まったところでした。

また、実際に使うときは頻繁にアクセスすることになるであろう入力ファイルの書式は少し拘ったところです。
元のコードではjsonが使用されていましたが、パラメータの設定などを詳しく説明しようと思うとどうしてもコメント機能は欲しくなるところ。調べてみた範囲では「YAML」や「json5」等を使用するのが解決策として挙げられていましたが、今回はシンプルにjsonにコメント機能が拡張されただけというのが分かりやすくまたライブラリも既にあったためjson5を採用しました。結果として入力ファイルの書式の自由度が上がり、かなり利便性が高まったと感じました！
更に、入力ファイルを作成していると良く出くわす問題が「色々なところに同じ数字を入れる問題」でした。例えばイベント時刻（「第1段の燃焼が〇〇秒に終了」の〇〇秒に当たる情報）等は良く参照されやすい値ですが、コピペで運用しているとどこかを変更し忘れてしまう等して計算を誤ってしまう可能性がありました。これではあまり良くないので、別途json5のラッパークラスを作成し、"&<参照したい変数名>" の記法で任意の他の設定項目の値を参照する機能を追加してみました。
これにより、例えば第1段燃焼終了(B1BO)のイベント時刻を記した場所があった場合、"&B1BO"と書くことで当該の情報を読み取りに行くような動作をしてくれますので、繰り返し同じ情報を書かなくても良くなりました。

その他も細々とした工夫をしていますが、大きなところとしては上に挙げた内容でほとんどだと思います。もし実装箇所で不明な点等がありましたらコメント等で教えて頂ければ幸いです。

## 4. どんな成果があったか？
今回移植してみて以下の効果があったと思いました。

* 今まで漠然と「運動方程式を解いているんだよねー」と思っていた解析の中身が良く分かるようになった！
* Pythonが遅い、というのはやはりあると思った。一方でcProfile等のライブラリも充実しているので高速化すべき対象の処理を絞り込むのは意外と簡単にできた。
    * 特にfor文を用いた処理は本当に遅い。numpyを用いたり、行列計算等に置き換えたりすると早くなることが多かった。
* 全体の設計次第でプログラムの効率や可読性は大きく変わると思った。一度で完璧な設計を行うのは難しいので、何度かリファクタリングを繰り返したり、将来性を見越して拡張性を持たせた設計をすることが重要だと思った。
* 全体を通して一つのプログラムを書ききったことで、少し自信が付いた！

飛行解析は実際やってみると面白い分野なのですが、自由に扱えるツールやHow to的な情報等が不足しているとも感じています。筆者としては、本記事ないし本プログラムから少しでも興味を持っていただけたようであれば嬉しい限りです。


## 5. おわりに
この記事を書いている今ですが、実際にコードを書き上げたときからは大分時間が経ってしまっておりまして、記憶が抜けてしまっているところも多々あり色々と苦労しながら何とか書き上げたところでした。そのためもしかすると不正確な内容等があるかもしれませんが、GitHub上のソースコードの方が正ですので、もし差分や誤り等があった場合はこっそり教えて頂けると嬉しいです。

当時は、自分が忘れたときのことも考えつつなるべく初見の人が見ても分かるようなプログラムを書こう！という思いでやっていたはずなのですが、今見てみると説明が不足していたり、これどうなってるんだっけ？みたいな部分が意外と多かったりで、自分の甘さを再認識することが出来ました。
単体の関数の説明、特に引数・返り値の詳細等が抜けているところは勿論そうですが、ソースコード内のコメントとしては書かないような内容、例えばプログラム全体の構成や大まかな処理の流れといったところは別途ドキュメントとして残しておかないと結構混乱してしまうなという教訓を得たところです。
また良かれと思ってやったInterfaceを用いたクラス設計ですが、コードが煩雑になってイマイチ何をやり取りしているか分かりづらくなってしまったかなという印象を今見ると持ちました。この辺の上手い設計についてももう少し勉強したいと思ったところです。

なお、参照させて頂いた元記事の方にご紹介がある通り、飛行解析ツールを用いることでロケットの飛行経路等を求めることが出来るようになります。もしご興味がありましたら是非実際にローカル環境等で実行してPC上で自由にロケットを飛ばしてみてください！

最後になりますが、本記事を執筆することが出来たのも、素晴らしいソフトウェアをオープンソースという形で公開してくださった先人の方がいらっしゃたからこそと感じております。
製作者のインターステラテクノロジズ株式会社様にはこの場を借りてお礼申し上げたいと思います。

それではこの辺にて、ここまでの長文にお付き合い頂き、誠にありがとうございました。