エスケープする低推力スパイラルは、C3 = 0の飛行経路角（ガンマ）は常に39度ですか？

@ MarkAdlerのコメントは、
0のC3からの遅いスパイラルは、衝動的な操作の約2.4倍のΔVを必要とするのはなぜですか？この綿密で効率的な@MarkAdler
答えが、別の思いやりのある答え。

その答えの下に、もう一つの

イースターエッグ（コメント宝石。

常に速度ベクトルに揃えられます。これは、比エネルギーを増加させるために推力を最も効率的に使用することです。最後のγは31°です。

この回答では、@
Julioは、測定する$ beta $と$ gamma
$の両方の角度の定義を示す図を提供しています瞬時速度ベクトルと半径方向および接線方向との間の角度をそれぞれ示す。

この回答で@TomSpilkerはこれらの角度を詳しく説明し、この回答私はそれらを計算する方法についてもう少し詳しく説明します。

今私は戻って様々な条件を使って低推力下で外向きに螺旋軌道を計算しました。いつもC3 =
0（31度ではない）の瞬間をチェックすると、最終的な角度$ gamma $（ガンマ）が約39度になります。

私は、GM = 1.0でr = 1.0軌道の周期が$ 2 pi $である無単位計算をしています。この場合、C3 = v ^
2 – 2/rである。

note: For this calculation, thrust is always in
the same direction as velocity $mathbf{v}$, rather than in the
tangential direction (perpendicular to $mathbf{r}$) and I’m
beginning to wonder if herein lies the difference between 31 and 39
degrees.

Question: Is this ~39 degrees at C3=0 correct,
and is it expected to be invariant like this?

      starting conditions                              at C3 = 0
-------------------------------     ------------------------------------------
rstart  vstart    C3    thrust      time   delta-v  gamma(deg)    r       v        C3
 1.0     1.0    -1.0    0.01        74.5    0.745     38.9       8.78    0.477   0.000
 1.0     1.0    -1.0    0.001       856.3   0.856     39.2      27.80    0.268   0.000
 1.0     1.0    -1.0    0.0001      9192.1  0.919     39.2      87.91    0.151   0.000
 4.0     0.5    -0.25   0.0001      4192.1  0.419     39.1      87.90    0.151   0.000

def deriv(X, t):
    x, v  = X.reshape(2, -1)
    vnorm = v/np.sqrt((v**2).sum())
    acc_g = -x * ((x**2).sum())**-1.5
    acc_t = thrust * vnorm
    return np.hstack((v, acc_g + acc_t))

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
degs, rads        = 180/pi, pi/180

T    = 16 * twopi        # or 160, 1600

ntot = 20001
time = np.linspace(0, T, ntot)

rstart = 1.0             # or 4.0
vstart = np.sqrt(1./rstart)

X0     = np.array([rstart, 0, 0, vstart])

thrust = 0.01            # or 0.001, 0.0001

answer, info = ODEint(deriv, X0, time, full_output= True)

xx, vv = answer.T.reshape(2, 2, -1)

r   = np.sqrt((xx**2).sum(axis=0))
vsq =         (vv**2).sum(axis=0)
C3 = vsq - 2./r

nstop = np.argmax(C3>0) + 1

dotted     = (xx*vv).sum(axis=0)
rabs, vabs = [np.sqrt((thing**2).sum(axis=0)) for thing in (xx, vv)]
gamma      = np.arcsin(dotted/(rabs*vabs))   # Per Tom Spilker's answer Eq. 3

print 'C3 min, max: ', C3.min(), C3.max()
print 'nstop, ntot: ', nstop, ntot
if True:
    plt.figure()

    plt.subplot(1, 2, 1)
    plt.plot(xx[0, :nstop], xx[1, :nstop])

    plt.subplot(3, 2, 2)
    plt.plot(time[:nstop], r[:nstop])
    plt.ylabel('r')

    plt.subplot(3, 2, 4)
    plt.plot(time[:nstop], C3[:nstop])
    plt.plot(time[:nstop], np.zeros_like(C3)[:nstop], '-k')
    plt.ylabel('C3')

    plt.subplot(3, 2, 6)
    plt.plot(time[:nstop], degs*gamma[:nstop])
    plt.ylabel('gamma (deg)')

    plt.suptitle('thrust = 0.0001, start at r=4, time=4192.1, gamma=39.12 deg, r=87.90', fontsize=16)

    plt.show()