前回は地球の重力加速度g一定を仮定、空気抵抗０での「弾道」でした。今回は万有引力の法則っす。ニュートン大先生様であります。太陽と地球（ともに質点）をおいて太陽の周りを地球に回っていただこうという趣向です。本当は太陽も微妙に動くのだけれどその動きはネグってしまって原点固定の単純化。そんなんでも１周は長いのよ。どうする？
※「ブロックを積みながら」投稿順 index はこちら
※動作確認にはWindows 11の パソコン(64bit)上にインストールしたScilabの以下バージョンを使用しています。
Scilab 2024.0.0
万有引力の法則と太陽‐地球系の諸元
さて、いかに掲げますのは万有引力の法則、および今回Xcosで計算しようとする太陽、地球系の「定数」であります。
なお、rの太陽、地球間の平均距離、１天文単位(au)の数字がやけに細かいのは、これが「定義値」だからであります。
Xcos上での実装
簡単のため、太陽は３次元空間の原点に固定、としました。すると後は地球にかかる加速度と初期位置、初期速度が決まれば地球の軌道を描けるハズ。地球にかかる加速度 a とすると
Fr = m * a
なので、上記万有引力の法則から、a = -G*M/r^2 と求まります。GとMは定数なので、後はその時その時で変化するハズの位置毎にaを計算していけば軌道が描けるハズ。しかし問題が。
Xcos上での計算終わりませぬ。
基本Xcosの時間単位は秒なので、上記の定数どもをそのままスケーリングもせずにモデルに代入。どうせ倍精度型なのでテキトーに計算できるっしょ、ということで実行してみたところが、なかなか軌道が描けませぬ。実時間の１年はかからないにせよ、目の前でさらっと描いてほしいのだけれど。だいたいな感じ（いい加減だな）で良いので。
そこで、太陽-地球系の１時間＝3600秒を Xcos上での 1秒にマップすることにしました。これで大幅スピードアップ。ただし計算の目の粗さは3600倍になるのでいい加減。上記 a の値はメートル毎秒毎秒なので、それを3600の3600倍することで、メートル毎時毎時になる筈。ホントか？
実験したXcosフロー
フローが以下に。ｘ、ｙ、ｚの３次元空間の中での軌道を描くものです。左上の方にあるのが肝心な地球の３次元座標を計算する部分、右下の方にあるのが万有引力の法則つかって加速度を求めるところです。求めた加速度は再び３次元空間の各成分に分解してからでないと積分できないので、ちょいとメンドイです。
上記のフローに与える定数を計算するコンテキスト設定が以下に。
数字がカイデ～。メンドクセー。 先ほどの計算単位を１時間とする「スケーリング」が入ってよりメンドクなってます。なお、Vtanというのが初期の地球の軌道接線方向速度です。とりあえず軌道一周を円とみて円周を求め、それを時間で割れば速さが求まるっと。
なお初期位置は x、y、zを求める積分の初期値として仕込んでます。xだけ1au相当。yとzはゼロ。上記ではx, z方向の速度も0。フロー自体は３次元計算可能にはなっているけれども、この設定では x、y 平面で回転することになる筈。なんだ、実質２次元じゃん。
シミュレーション結果
１年＝３６５日 x ２４時間 = ８７６０時間ということでシミュレーションした結果が以下に（Xcos上では8760sという設定。）スケーリングのお陰で超速。
ほぼほぼ円になった感じ。ならんと困るが。