@MarkAdler's comment led me to ask Why would a slow spiral from a C3 of zero take about 2.4 times as much ΔV as an impulsive maneuver? which resulted in this tidy and efficient @MarkAdler answer which points to another thoughtful answer about slowly spiraling out from a circular orbit to escape in the limit of very weak prograde propulsion, which (at first counterintuitively) slows you down while raising your orbit.
Below that answer is yet another easter-egg-like comment gem.
Always aligned with the velocity vector. That is the most efficient use of the thrust in order to increase the specific energy. The final γ is 31°.
In this answer @Julio provides a diagram showing definitions for both $\beta$ and $\gamma$ angles which measure the angle between the instantaneous velocity vector and the radial and the tangential directions, respectively.
In this answer @TomSpilker elaborates on these angles, and in this answer I give a little more information on how to calculate them.
Now I've gone back and calculated an outwardly spiraling orbit under low thrust using various conditions. Invariably I end up with a final angle $\gamma$ (gamma) of about 39 degrees when checking the moment where C3 = 0, not 31 degrees.
I'm doing a unitless calculation where GM = 1.0 and the period of an r=1.0 orbit is $2 \pi$. In this case 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()
