Orbital Mechanics & Astrodynamics

Question 1 (12)

(a) $\dfrac{d}{dt}(\mathbf r\times\dot{\mathbf r}) = \dot{\mathbf r}\times\dot{\mathbf r} + \mathbf r\times\ddot{\mathbf r} = 0 + \mathbf r\times(-\tfrac{\mu}{r^3}\mathbf r) = 0$. (1) for product rule, (1) first term zero, (1) second term zero because $\mathbf r\parallel\ddot{\mathbf r}$. Hence $\mathbf h$ constant → motion planar.

(b) Take $\ddot{\mathbf r}\times\mathbf h = -\tfrac{\mu}{r^3}\mathbf r\times(\mathbf r\times\dot{\mathbf r})$. Using BAC–CAB: $\mathbf r\times(\mathbf r\times\dot{\mathbf r}) = \mathbf r(\mathbf r\cdot\dot{\mathbf r}) - \dot{\mathbf r}(r^2)$, and $\mathbf r\cdot\dot{\mathbf r}=r\dot r$. (2) So RHS $= -\mu\left(\dfrac{\dot r}{r^2}\mathbf r - \dfrac{\dot{\mathbf r}}{r}\right) = \mu\dfrac{d}{dt}\!\left(\dfrac{\mathbf r}{r}\right)$. (2) LHS $=\dfrac{d}{dt}(\dot{\mathbf r}\times\mathbf h)$ since $\mathbf h$ const. Integrate: $\dot{\mathbf r}\times\mathbf h = \mu\left(\dfrac{\mathbf r}{r}+\mathbf e\right)$, where $\mathbf e$ is the constant eccentricity vector. (1) Dot with $\mathbf r$: $\mathbf r\cdot(\dot{\mathbf r}\times\mathbf h)=(\mathbf r\times\dot{\mathbf r})\cdot\mathbf h = h^2$. RHS $=\mu(r + \mathbf r\cdot\mathbf e)=\mu(r + re\cos\theta)$. Solve: (1) $r=\frac{h^2/\mu}{1+e\cos\theta}.$ $\mathbf e$ points toward perigee; its magnitude is the eccentricity. (1)

(c) $e=0$ circle; $0<e<1$ ellipse; $e=1$ parabola; $e>1$ hyperbola. (2)

Question 2 (10)

(a) At any point $\varepsilon=\tfrac12 v^2-\tfrac{\mu}{r}$ (energy conservation). (1) Setting $\varepsilon=-\tfrac{\mu}{2a}$ (1): $\tfrac12 v^2-\tfrac{\mu}{r}=-\tfrac{\mu}{2a}$. (1) Solve for $v^2$: $v^2=\mu\!\left(\tfrac{2}{r}-\tfrac{1}{a}\right)$. (1)

(b) Areal velocity $\dfrac{dA}{dt}=\dfrac{h}{2}$ (Kepler II, from $\mathbf h$ const). (1) Over one period area swept = ellipse area: $\pi a b = \dfrac{h}{2}T$. (1) So $T=\dfrac{2\pi ab}{h}$. With $b=a\sqrt{1-e^2}$ and $h^2=\mu p=\mu a(1-e^2)$ → $h=\sqrt{\mu a(1-e^2)}$. (1) $T=\frac{2\pi a\cdot a\sqrt{1-e^2}}{\sqrt{\mu a(1-e^2)}}=\frac{2\pi a^{3/2}}{\sqrt\mu}. \textbf{(2)}$ Square: $T^2=\dfrac{4\pi^2}{\mu}a^3$. (1)

Question 3 (12)

$v_1=\sqrt{\mu/r_1}=\sqrt{3.986\times10^5/6778}=7.6686\ \text{km/s}$. (1) $v_2=\sqrt{\mu/r_2}=\sqrt{3.986\times10^5/42164}=3.0747\ \text{km/s}$. (1)

(b) $a_t=\tfrac12(r_1+r_2)=\tfrac12(48942)=24471\ \text{km}$. (1) Perigee speed $v_p=\sqrt{\mu(2/r_1-1/a_t)}=\sqrt{3.986\times10^5(2/6778-1/24471)}$ $=\sqrt{3.986\times10^5(2.9508\times10^{-4}-4.0865\times10^{-5})}=\sqrt{3.986\times10^5\cdot2.5421\times10^{-4}}=10.064\ \text{km/s}$. (1.5) Apogee speed $v_a=\sqrt{\mu(2/r_2-1/a_t)}=\sqrt{3.986\times10^5(4.7434\times10^{-5}-4.0865\times10^{-5})}=\sqrt{3.986\times10^5\cdot6.569\times10^{-6}}=1.6183\ \text{km/s}$. (1.5)

(c) $\Delta v_1=v_p-v_1=10.064-7.6686=2.395\ \text{km/s}$. (1.5) $\Delta v_2=v_2-v_a=3.0747-1.6183=1.456\ \text{km/s}$. (1.5) Total $\Delta v=3.852\ \text{km/s}$. (1)

(d) $t=\pi\sqrt{a_t^3/\mu}=\pi\sqrt{24471^3/3.986\times10^5}=\pi\sqrt{3.674\times10^7}=\pi\cdot6061\approx19040\ \text{s}\approx5.29\ \text{h}$. (2)

Question 4 (10)

(a) Mean anomaly $M=n(t-t_p)$ increases uniformly, $n=2\pi/T$. (1) Define eccentric anomaly $E$ via the auxiliary circle of radius $a$. Using Kepler II, area swept from perigee is proportional to time: $\dfrac{M}{2\pi}=\dfrac{\text{area}}{\pi a^2}$. (1) The area (sector on ellipse) computed from the circle geometry = $\tfrac12 a^2(E-e\sin E)$. (1) Equating: $M=E-e\sin E$. (1)

(b) Pseudocode: (4, 1 each for guess, formula, loop, stop)

def kepler(M, e, tol=1e-10):
    E = M if e < 0.8 else 3.14159   # initial guess
    while True:
        dE = (E - e*sin(E) - M) / (1 - e*cos(E))
        E -= dE
        if abs(dE) < tol:
            return E

Iteration: $E_{n+1}=E_n-\dfrac{E_n-e\sin E_n-M}{1-e\cos E_n}$.

(c) $E_0=0.8$. $f=0.8-0.3\sin0.8-0.8=-0.3\sin0.8=-0.3(0.71736)=-0.21521$. $f'=1-0.3\cos0.8=1-0.3(0.69671)=0.79099$. (1) $E_1=0.8-(-0.21521/0.79099)=0.8+0.27208=1.0721\ \text{rad}$. (1)

Question 5 (10)

(a) Formula stated. (1) The oblate bulge adds mass near the equator; the resulting non-central torque-free but non-Keplerian potential makes the orbital plane precess about Earth's spin axis, so $\Omega$ drifts. (2) At $i=90°$ (polar), $\cos i=0$ ⇒ symmetry makes net nodal torque zero, $\dot\Omega=0$. (1)

(b) SSO requires $\dot\Omega = +1.991\times10^{-7}\ \text{rad/s}$ (matching Earth's mean motion about Sun, $2\pi/\text{year}$), so the orbital plane keeps a fixed angle to the Sun. (2) Since $\dot\Omega>0$ requires $\cos i<0$, the orbit must be retrograde, $i>90°$ (typically ~98°). (1)

(c) $n=\sqrt{\mu/a^3}=\sqrt{3.986\times10^5/7078^3}=\sqrt{3.986\times10^5/3.546\times10^{11}}=1.0603\times10^{-3}\ \text{rad/s}$. (1) $\cos i = \dfrac{\dot\Omega}{-\tfrac32 J_2 n (R_\oplus/a)^2}$ (with $e=0$). $(R_\oplus/a)^2=(6378/7078)^2=0.8119$. Denominator $=-\tfrac32(1.0826\times10^{-3})(1.0603\times10^{-3})(0.8119)=-1.3980\times10^{-6}$. (1) $\cos i=1.991\times10^{-7}/(-1.3980\times10^{-6})=-0.14242$ ⇒ $i=\arccos(-0.14242)=98.19°$. (1)

Question 6 (6)

(a) Bi-elliptic uses an intermediate apogee $r_b\gg r_2$; the plane-change/velocity-change burns happen where orbital speed is small, cutting the required $\Delta v$ despite three burns. It beats Hohmann when the ratio $r_2/r_1\gtrsim 11.94$ (and always for $>15.58$), at the cost of much longer transfer time. (3)

(b) In the planet's frame the encounter is a hyperbola; energy is conserved, so incoming and outgoing $v_\infty$ have equal magnitude — only the direction rotates. In the Sun's frame you add the planet's heliocentric velocity to $\mathbf v_\infty$; because the direction changed, the vector sum has a different magnitude, so the spacecraft gains (or loses) heliocentric speed while the planet loses (gains) a negligible amount. (3)

[
  {"claim":"Hohmann v1 circular LEO", "code":"mu=3.986e5; r1=6778.0; v1=sqrt(mu/r1); result = abs(float(v1)-7.6686)<0.01"},
  {"claim":"Hohmann total delta-v ~3.852 km/s", "code":"mu=3.986e5; r1=6778.0; r2=42164.0; at=(r1+r2)/2; v1=sqrt(mu/r1); v2=sqrt(mu/r2); vp=sqrt(mu*(2/r1-1/at)); va=sqrt(mu*(2/r2-1/at)); dv=(vp-v1)+(v2-va); result = abs(float(dv)-3.852)<0.01"},
  {"claim":"Transfer time ~19040 s", "code":"import sympy as sp; mu=3.986e5; at=24471.0; t=sp.pi*sqrt(at**3/mu); result = abs(float(t)-19040)<50"},
  {"claim":"Kepler one Newton step E1~1.0721", "code":"M=0.8; e=0.3; E0=0.8; f=E0-e*sin(E0)-M; fp=1-e*cos(E0); E1=E0-f/fp; result = abs(float(E1)-1.0721)<0.001"},
  {"claim":"SSO inclination ~98.19 deg", "code":"mu=3.986e5; a=7078.0; Re=6378.0; J2=1.0826e-3; Om=1.991e-7; n=sqrt(mu/a**3); c=Om/(-1.5*J2*n*(Re/a)**2); i