ホーマン軌道の計算

質問

ホーマン軌道の計算について質問しようとしたのですが、
「当サイトの利用規約・ルール」に反する、という内容の mail が送られてきて、
質問自体が削除されてしまいました。
たぶん、盛んに URL 参照先を書いたのがいけなかったのでしょう。

それはともかく、MAXIMA (Web 検索で「MAXIMA 数式処理」を探して下さい) を使いました。
mu : 1 $ r_1 : 1 $
a_1 : 1 $ a_h : k*a_1 $
V_1h : sqrt(mu*(2/r_1 - 1/a_h)) $
V_1 : sqrt(mu*(2/r_1 - 1/a_1)) $

で、ΔＶ ＝ V_1h - V_1 を求めた場合と、r = a*(1 - e*cos(E)) 　(E は離心近点離角)
で求めた場合の数値が異なります。どうして違うのか、その理由がわかりません。
たぶん勘違いしているのだと思いますが、ヒントなど、お教えいただければ幸いです。

詳細は、URL 参照先を書けないので、おいおいなんとか抜け道を見つけます。あしからず。

回答No.4

ojisan7 · ベストアンサー

たびたび、すみません。

＞ホーマン遷移軌道 (楕円軌道) の長軸の長さの半分は、最終目標円軌道の半径と等しくはありません。

と、ありますが、最終目標円軌道の半径と等しいと思います。
大切なのは、
a_h＝(r_max+r_min)/2　・・・◎
です。従って、hogehoneさんが、MAXIMAを使って得た、
★１と★２の計算は、計算として正しいと思います。
★１と★２の結果が違って見えたのは、単に、◎の式を考慮していなかったためだと思われます。

お礼

ありがとうございます。

Web 検索の、Web 文書検索でも Webイメージ (image) 検索でも、ホーマン軌道を検索すると、

「 a_h/a_1 = r_max/r_min = 1.17770　・・ 中心天体からの距離最大・最小の比 ★２ 」
かと思ってしまう。それは誤解で、そうではなくて、

その解釈は間違いなのですね。

こっちのことを言っているのではないかと誤解してしまいます。(そのような説明図がたくさんある)

回答No.3

ojisan7

すみません。タイプミスをしてしまいました。

【・・・ということですが、（２）は、
v_max/v_min ＝1.177700348432056 　
であり、（２）は、計算すると、確かに、
a_h/a_1 ＝1.08885・・
となりますが、・・・】の部分を、

【・・・ということですが、★２は、
r_max/r_min ＝1.177700348432056 　
であり、　★１は、計算すると、確かに、
a_h/a_1 ＝1.08885・・
となりますが、
a_h＝(r_max+r_min)/2ですから、
a_h/a_1 ＝(r_max+r_min)/2 *1/r_min
＝(1+r_max/r_min )/2=(1+1.177700348432056)/2 　
＝1.08885・・
となり、つじつまが合います。
】と、訂正させてください。

さて、hogehoneさんの【a_h】の解釈はちょっと、勘違いされているのかも知れません。【a_h】はホーマン遷移軌道の長軸の半分の長さ（最終目標円軌道の半径）だと思います。

お礼

ありがとうございます。

＞ 「hogehoneさんの【a_h】の解釈はちょっと、勘違いされているのかも知れません。
＞ 【a_h】はホーマン遷移軌道の長軸の半分の長さ（最終目標円軌道の半径）だと思います。」

　しかし、

ホーマン遷移軌道 (楕円軌道) の長軸の長さの半分は、最終目標円軌道の半径と等しくはありません。

　また、

Web 検索、文書検索でもイメージ (image) 検索でも、ホーマン軌道を検索すると、そうではなくて、

「 a_h/a_1 = r_max/r_min = 1.17770　・・ 中心天体からの距離最大・最小の比 ★２ 」

こっちのことを言っているのではないかと誤解してしまいます。(そのような説明図がたくさんある)
もう一回考え直して計算し、解釈し直してみます。

(まず第一段階の、円軌道からホーマン軌道に移る場面の話をしています。このへんもよろしく。)

回答No.2

ojisan7

> (1) ΔＶ = V_1h - V_1 = sqrt(mu*(2/r_1 - 1/a_h)) - sqrt(mu*(2/r_1 - 1/a_1))
で計算した値と、Web 上でよく見つかる、
(2) 行列計算による方法　の結果が異なるのです。

ということですが、（２）は、
v_max/v_min ＝1.177700348432056 　
であり、（２）は、計算すると、確かに、
a_h/a_1 ＝1.08885・・
となりますが、
a_h＝(v_max+v_min)/2ですから、
a_h/a_1 ＝(v_max+v_min)/2 *1/v_min
＝(1+v_max/v_min )/2=(1+1.177700348432056)/2 　
＝1.08885・・
となり、つじつまが合います。

お礼

ありがとうございます。ご回答の、
＞・・・であり、（２）は、計算すると、確かに、a_h/a_1 ＝1.08885・・となりますが、

これは、(2) での計算では a_h/a_1 = r_max/r_min = 1.17770・・
中心天体からの距離最大・最小の比 ★２　の誤りではありませんでしょうか。

ホーマン軌道を使った軌道制御の第一段階の話をしています。よく紹介されている、
ΔＶ = V_1h - V_1 = sqrt(mu*(2/r_1 - 1/a_h)) - sqrt(mu*(2/r_1 - 1/a_1))　― (1)

これは、半径 a_1 の円軌道から a_h (a_h は増速後の楕円軌道の遠地点の中心天体からの距離)
軌道に移るのに必要な増速量 ΔＶ_1 を表している式ではない、とおっしゃっている意味に
取れてしまうのですが、これは解釈間違いでしょうか。

注： (2) 式と数値積分による方法で得られた値は等しいのですが、しかしそれらは (1)
式で求めた値と異なる、というのが、質問の発端です。過去レス等も併せてご覧下さい。

補足

間違い訂正させてください。
＞軌道に移るのに必要な増速量 ΔＶ_1 を表している式ではない

＞軌道に移るのに必要な増速量 ΔＶ を表している式ではない

(ついつい、ホーマン軌道を使った軌道修正の第２段階の話が出てきそうで　"_1"　などと書いてしまいました。ごめんなさい。)

回答No.1

sqw-99

初歩的なミスで、度数をラジアン変換せずに打ち込んでないかい？
１度＝０．０１７RADだと思った。
（Ｅ*（π／３６０））
でラジアン変換出来る。

違ったらごめんね。＾＾

お礼

ありがとうございます。
print(sin(%pi/2), sin(%pi), cos(%pi/2), cos(%pi));　を実行すると、
1 0 0 -1　と答えてくれます。

● どなたか答えていたければ、URL 含みのレスもできる～ かな。
(1) ΔＶ = V_1h - V_1 = sqrt(mu*(2/r_1 - 1/a_h)) - sqrt(mu*(2/r_1 - 1/a_1))
で計算した値と、Web 上でよく見つかる、
(2) 行列計算による方法　の結果が異なるのです。

★１ と ★２ の結果の数値が違う、その理由がよくわかりません。よろしくお願い致します。

なお、別に数値積分による方法も試してみたのですが、その結果は (2) の場合と等しく
なりました。

補足

[回答に対するお礼] を書いたのですが、その中の参照 URL は削除されてしまいました。
しょうがないので、ここに展開し書きます。(MAXIMA を使いました。
Web 検索で「数式処理 MAXIMA」とかやると、けっこう Hit します)

● (1) の方法
/*
5.txt
円軌道からホーマン軌道に移る際に必要な増速量ΔＶと
a_h/a_1 (軌道半径の比) の関係を求める表を作る。
*/
mu : 1 $
r_1 : 1 $
a_1 : 1 $
a_h : k*a_1 $
V_1h : sqrt(mu*(2/r_1 - 1/a_h)) $
V_1 : sqrt(mu*(2/r_1 - 1/a_1)) $
fortran(ev( V_1h - V_1 )) $
fortran(solve(dltV=SQRT(2 - 1/k) - 1, k));
positive;
fa : -1/(dltV^2 + 2*dltV - 1) $
/* plot2d(fa, [dltV, 0, 0.05]) $ */
print("dltV a_h/a_1") $
for dltV : 0 thru 0.05 step 0.01 do print(dltV, ev(fa)) $

● (2) の方法
/*
6.txt
楕円軌道の計算
μ (mu), a, e は固定し、離心近点離角 E をパラメーターとして位置・速度を求める。
*/
/*
関数
*/
f1(a,e,E) := a*(cos(E) - e) $
f2(a,e,E) := a*sqrt(1 - e^2)*sin(E) $
f3(a,e,E) := sqrt(f1(a,e,E)^2 + f2(a,e,E)^2) $
f4(a,e,E) := a*(1 - e*cos(E)) $
f5(mu,a,r,E) := -(sqrt(mu*a)/r)*sin(E) $
f6(mu,a,r,E) := (sqrt(mu*a)/r)*sqrt(1 - e^2)*cos(E) $
f7(mu,a,r,E) := sqrt(f5(mu,a,r,E)^2 + f6(mu,a,r,E)^2) $
f8(mu,a,r) := sqrt(mu*(2/r - 1/a)) $
/*
定数
*/
e:0.0816 $
a:1/(1-e) $
mu:1 $
/*
角運動量保存の確認 (E = 0, π/2, 2π/2, 3π/2、以下同様)
*/
print("Anglular momentum (L)") $
for E0:0 thru 3*%pi/2 step %pi/2 do
print(E0/(%pi/2), "|", f1(a,e,E0)*f6(mu,a,f4(a,e,E0),E0) - f2(a,e,E0)*f5(mu,a,f4(a,e,E0),E0)) $
/*
速度ベクトル Vx, Vy と絶対値
*/
print("Velocity") $
for E0:0 thru 3*%pi/2 step %pi/2 do
print(E0/(%pi/2), "|", f5(mu,a,f4(a,e,E0),E0),",",f6(mu,a,f4(a,e,E0),E0),
" :",f7(mu,a,f4(a,e,E0),E0)) $
print("v_max/v_min", f8(mu,a,f4(a,e,0))/f8(mu,a,f4(a,e,%pi))) $
/*
位置ベクトル Px, Py と絶対値
*/
print("Coordinates") $
for E0:0 thru 3*%pi/2 step %pi/2 do
print(E0/(%pi/2), "|", f1(a,e,E0), f2(a,e,E0),":", f3(a,e,E0)) $
print("r_max/r_min", f4(a,e,%pi)/f4(a,e,0)) $

● 実行結果とコメント

計算にはフリーの数式処理プログラム MAXIMA を使用しました。Web 検索「数式処理 MAXIMA」で
結構ヒットします。ダウンロード先はこちらです。http://maxima.sourceforge.net/download.shtml

"5.txt" を MAXIMA のコンソールにコピペして実行すると、次のような結果が得られました。
(一部手動 edit あり)

dltV a_h/a_1
-----------------------
0 1
0.01 1.020512297173181
0.02 1.042100875364735
0.03 1.06484932382068
0.04 1.088850174216028
0.05 1.114206128133705

これによるとΔＶが 0.04 のとき、a_h/a_1 は 1.08885・・ とわかります。★１

続いて "6.txt" を同様に実行します。r = 1 で増速量が 0.04 になるように a, e を設定しています。

Anglular momentum (L) 　角運動量保存の確認。E = 0, π/2, 2π/2, 3π/2 としています。以下同様。
0 | 1.04
1 | 1.04
2 | 1.04
3 | 1.04 　確かに皆同じになっています。

Velocity 　速度ベクトル Vx, Vy と絶対値
0 | 0, 1.04 : 1.04 　　　　　　速度は 1+0.04、0.04 が増速分
1 | -0.95833188405688, 0 : 0.95833188405688
2 | 0, -0.88307692307692 : 0.88307692307692
3 | 0.95833188405688, 0 : 0.95833188405688

v_max/v_min 1.177700348432056 　速度最大・最小の比

Coordinates 　位置ベクトル Px, Py と絶対値
0 | 1.0 0 : 1.0
1 | -0.08885017421603 1.085219032468587 : 1.088850174216028 　注*1 参照
2 | -1.177700348432056 0 : 1.177700348432056 　なんと、1.177・・です。★２
3 | -0.08885017421603 -1.085219032468587 : 1.088850174216028

r_max/r_min 1.177700348432056 　地心距離最大・最小の比 ★２

同じ増速量ΔＶを設定したつもりなのですが、★１と★２の値が異なっています。
この理由がよくわかりません。何かヒントをお教えいただけると幸いです。

注*1　この値はなぜか★１の値と同じです。しかしこの位置は近地点から π/2 しか進んでいません。