"""Verification of every identity and number in q6-sol.qmd (Midsem 2026, Q6). No normalised units anywhere: the algebra is carried in mu, r1, ra, rp, p, e, vc, exactly as in the write-up. Part 1 symbolic identities (sympy): velocity components, apsidal relations, the two transfer routes in terms of r1, ra, rp, the comparison inequalities, small-e series. Part 2 independent numerical checks, in km and km/s with the real mu: the orbits are integrated as an ODE (no conic formulas) and the burns are read off as vector differences; the transfer speeds come from vis-viva. Part 3 the numbers quoted in the write-up (km/s, minutes, table in e). Run: python3 q6-verify.py """ import numpy as np import sympy as sp from scipy.integrate import solve_ivp # ---------------------------------------------------------------- Part 1 e, p, mu, f = sp.symbols('e p mu f', positive=True) r1, ra, rp = sp.symbols('r1 ra rp', positive=True) h = sp.sqrt(mu*p) # F1, F2: velocity components on a conic r = p/(1 + e*sp.cos(f)) fdot = h/r**2 vr = sp.simplify(sp.diff(r, f)*fdot) # r-dot = (dr/df) f-dot vt = sp.simplify(r*fdot) # r f-dot = h/r assert sp.simplify(vr - sp.sqrt(mu/p)*e*sp.sin(f)) == 0 assert sp.simplify(vt - sp.sqrt(mu/p)*(1 + e*sp.cos(f))) == 0 v2 = sp.simplify(vr**2 + vt**2) assert sp.simplify(v2 - mu/p*(1 + 2*e*sp.cos(f) + e**2)) == 0 a = p/(1 - e**2) assert sp.simplify(v2 - mu*(2/r - 1/a)) == 0 # vis-viva agrees # F3, F3b, F4: apsidal radii and h^2 = 2 mu rp ra / (rp + ra) rp_, ra_ = p/(1 + e), p/(1 - e) assert sp.simplify((rp_ + ra_)/2 - a) == 0 assert sp.simplify(rp_*ra_/a - p) == 0 # p = rp ra / a assert sp.simplify(2*rp_*ra_/(rp_ + ra_) - p) == 0 # p = 2 rp ra/(rp+ra) assert sp.simplify(2*mu*rp_*ra_/(rp_ + ra_) - mu*p) == 0 # h^2 # the circle and the target, with p = r1 vc = sp.sqrt(mu/r1) hT = sp.sqrt(mu*r1) # h of the target = h of the circle assert sp.simplify(hT - r1*vc) == 0 tgt = {r1: r1} # (p = r1 is imposed by substitution below) ra_T, rp_T = r1/(1 - e), r1/(1 + e) assert sp.simplify(hT/rp_T - vc*(1 + e)) == 0 # target speed at periapsis assert sp.simplify(hT/ra_T - vc*(1 - e)) == 0 # target speed at apoapsis # ---- route (i): transfer apsides r1 and ra h1 = sp.sqrt(2*mu*r1*ra/(r1 + ra)) # F4 rho = (r1 + ra)/ra # = 1 + r1/ra assert sp.simplify(sp.simplify(h1/r1)**2 - vc**2*2/rho) == 0 # speed at r1: vc sqrt(2/rho) v_r1 = vc*sp.sqrt(2/rho) assert sp.simplify(h1/r1 - v_r1) == 0 assert sp.simplify(h1/ra - (r1/ra)*v_r1) == 0 # speed at ra: (r1/ra) v_r1 (= (rho - 1) v_r1) burnA1 = v_r1 - vc burnB1 = h1/ra - hT/ra # transfer speed at ra minus target speed at ra assert sp.simplify(burnB1 - (r1/ra)*burnA1) == 0 total1 = burnA1 + burnB1 assert sp.simplify(total1 - vc*(sp.sqrt(2*rho) - rho)) == 0 # general result in r1, ra # ... with ra = r1/(1-e): assert sp.simplify(rho.subs(ra, ra_T) - (2 - e)) == 0 assert sp.simplify((r1/ra).subs(ra, ra_T) - (1 - e)) == 0 # ---- route (ii): transfer apsides rp and r1 h2 = sp.sqrt(2*mu*rp*r1/(rp + r1)) sig = (rp + r1)/rp # = 1 + r1/rp v_r1b = vc*sp.sqrt(2/sig) assert sp.simplify(h2/r1 - v_r1b) == 0 # speed at r1 (apoapsis of the transfer) assert sp.simplify(h2/rp - (r1/rp)*v_r1b) == 0 # speed at rp: (r1/rp) v = (sigma - 1) v burnA2 = vc - v_r1b burnB2 = hT/rp - h2/rp assert sp.simplify(burnB2 - (r1/rp)*burnA2) == 0 total2 = burnA2 + burnB2 assert sp.simplify(total2 - vc*(sig - sp.sqrt(2*sig))) == 0 assert sp.simplify(sig.subs(rp, rp_T) - (2 + e)) == 0 assert sp.simplify((r1/rp).subs(rp, rp_T) - (1 + e)) == 0 # ---- shape of the transfer ellipses T_i, T_ii (e_i = e1, p_i = p1, e_ii = e2, p_ii = p2 here; NOT the target's e, p = r1) assert sp.simplify((ra_T - rp_T)/(ra_T + rp_T) - e) == 0 # target: e from its apsides e1, p1 = (ra - r1)/(ra + r1), h1**2/mu e2, p2 = (r1 - rp)/(r1 + rp), h2**2/mu assert sp.simplify(e1 - (2 - rho)/rho) == 0 and sp.simplify(p1 - 2*r1/rho) == 0 assert sp.simplify(p1/(1 + e1) - r1) == 0 and sp.simplify(p1/(1 - e1) - ra) == 0 # orbit eq. returns the apsides assert sp.simplify(e2 - (sig - 2)/sig) == 0 and sp.simplify(p2 - 2*r1/sig) == 0 assert sp.simplify(p2/(1 + e2) - rp) == 0 and sp.simplify(p2/(1 - e2) - r1) == 0 assert sp.simplify(e1.subs(ra, ra_T) - e/(2 - e)) == 0 # e1 = e/(2-e) assert sp.simplify(e2.subs(rp, rp_T) - e/(2 + e)) == 0 # e2 = e/(2+e) # in units of vc, as functions of e dv1 = sp.sqrt(2*(2 - e)) - (2 - e) dv2 = (2 + e) - sp.sqrt(2*(2 + e)) # the difference used to rank the two routes: (ii) - (i) = 4 - sqrt2 (sqrt(2+e) + sqrt(2-e)) S = sp.sqrt(2 + e) + sp.sqrt(2 - e) assert sp.simplify(dv2 - dv1 - (4 - sp.sqrt(2)*S)) == 0 assert sp.simplify(sp.expand(S**2) - (4 + 2*sp.sqrt(4 - e**2))) == 0 # small-e behaviour: e/2 -/+ e^2/16 + ... print('series (i) :', sp.series(dv1, e, 0, 4)) print('series (ii):', sp.series(dv2, e, 0, 4)) assert sp.limit(dv1, e, 1) == sp.sqrt(2) - 1 assert sp.limit(dv2, e, 1) == 3 - sp.sqrt(6) # single impulse at X (f = 90 deg): radial, energy and cosine-law cross-checks, all dimensional Vr, Vt = vr.subs(f, sp.pi/2).subs(p, r1), vt.subs(f, sp.pi/2).subs(p, r1) assert sp.simplify(Vr - e*vc) == 0 and sp.simplify(Vt - vc) == 0 eps_c = -mu/(2*r1) eps_T = -mu/(2*(r1/(1 - e**2))) # a = p/(1-e^2), p = r1 assert sp.simplify((eps_T - eps_c) - e**2*vc**2/2) == 0 # radial burn: d(eps) = dv^2/2 vT = vc*sp.sqrt(1 + e**2) # speed on the ellipse at X cos_g = vc/vT # cos(gamma) = v_theta / v assert sp.simplify(vc**2 + vT**2 - 2*vc*vT*cos_g - e**2*vc**2) == 0 print('Part 1: all symbolic identities hold') # ---------------------------------------------------------------- Part 2 (km, km/s) MU, R1 = 3.986e5, 1.0e4 VC = np.sqrt(MU/R1) def two_body(t, y): x, yy, vx, vy = y rr = np.hypot(x, yy) return [vx, vy, -MU*x/rr**3, -MU*yy/rr**3] def propagate(y0, t_end, events=None): return solve_ivp(two_body, [0, t_end], y0, rtol=1e-12, atol=1e-10, events=events, dense_output=True) def vis(rr, aa): return np.sqrt(MU*(2/rr - 1/aa)) def target_state_at_x(ev): """Target ellipse, started at periapsis (speed from vis-viva); state where it crosses x = 0 (f = +90 deg).""" rpT, raT = R1/(1 + ev), R1/(1 - ev) aT = (rpT + raT)/2 y0 = [rpT, 0, 0, vis(rpT, aT)] cross = lambda t, y: y[0] cross.direction, cross.terminal = -1, True return propagate(y0, 1e5, [cross]).y_events[0][0] for ev in (0.1, 0.3, 0.6, 0.9): X = target_state_at_x(ev) assert abs(np.hypot(X[0], X[1]) - R1) < 1e-4 # crossing lies on the circle r = r1 = p assert abs(X[2] + VC) < 1e-7 and abs(X[3] - ev*VC) < 1e-7 # v = (-vc, e vc): v_theta = vc, v_r = e vc dvv = X[2:] - np.array([-VC, 0.0]) # target velocity minus circular velocity assert abs(dvv[0]) < 1e-7 and abs(np.hypot(*dvv) - ev*VC) < 1e-7 # purely radial, magnitude e vc print('Part 2a: single impulse = e vc, radial (ODE-integrated target orbit)') def route_costs(ev): """Both routes: integrate the transfer ellipse, read each burn as a vector difference.""" rpT, raT = R1/(1 + ev), R1/(1 - ev) aT = (rpT + raT)/2 # route (i): leaves (r1, 0) on an ellipse with apsides (r1, ra) a1 = (R1 + raT)/2 v1 = vis(R1, a1) s1 = propagate([R1, 0, 0, v1], np.pi*np.sqrt(a1**3/MU), None) xa = s1.y[:, -1] # should be the apoapsis (-ra, 0) assert abs(xa[0] + raT) < 1e-3 and abs(xa[1]) < 1e-3 target_ap = np.array([0, -vis(raT, aT)]) # target velocity at its apoapsis, ccw A1, B1 = abs(v1 - VC), np.linalg.norm(target_ap - xa[2:]) # route (ii): leaves (-r1, 0) on an ellipse with apsides (rp, r1), ends at (rp, 0) a2 = (rpT + R1)/2 v2 = vis(R1, a2) s2 = propagate([-R1, 0, 0, -v2], np.pi*np.sqrt(a2**3/MU), None) xp = s2.y[:, -1] assert abs(xp[0] - rpT) < 1e-3 and abs(xp[1]) < 1e-3 target_pe = np.array([0, vis(rpT, aT)]) A2, B2 = abs(VC - v2), np.linalg.norm(target_pe - xp[2:]) return (A1, B1), (A2, B2), np.pi*np.sqrt(a1**3/MU), np.pi*np.sqrt(a2**3/MU) for ev in (0.1, 0.3, 0.6, 0.9): (A1, B1), (A2, B2), _, _ = route_costs(ev) rho_n, sig_n = 2 - ev, 2 + ev assert abs(A1 - VC*(np.sqrt(2/rho_n) - 1)) < 1e-7 and abs(B1 - (rho_n - 1)*A1) < 1e-7 assert abs(A2 - VC*(1 - np.sqrt(2/sig_n))) < 1e-7 and abs(B2 - (sig_n - 1)*A2) < 1e-7 assert abs(A1 + B1 - VC*float(dv1.subs(e, ev))) < 1e-7 assert abs(A2 + B2 - VC*float(dv2.subs(e, ev))) < 1e-7 print('Part 2b: burns of both routes from the integrated transfer orbits match the closed forms') es = np.linspace(1e-3, 1 - 1e-3, 2000) d1 = np.sqrt(2*(2 - es)) - (2 - es) d2 = (2 + es) - np.sqrt(2*(2 + es)) assert np.all(d1 < es) and np.all(d2 < es) and np.all(d1 < d2) print('Part 2c: (i) < (ii) < e on a 2000-point scan of 0 < e < 1') # ---------------------------------------------------------------- Part 3 ev = 0.3 print('\n--- numbers for e = 0.3, r1 = p = 10 000 km ---') a_t, rp_t, ra_t = R1/(1 - ev**2), R1/(1 + ev), R1/(1 - ev) print(f'a = {a_t:.1f} km, r_p = {rp_t:.1f} km, r_a = {ra_t:.1f} km') print(f'v_c = {VC:.4f} km/s; v_T(X) = {VC*np.hypot(1, ev):.4f} km/s; gamma = {np.degrees(np.arctan(ev)):.2f} deg') print(f'h_T = r1 vc = {R1*VC:.1f} km^2/s; target speeds: periapsis {VC*(1 + ev):.4f}, apoapsis {VC*(1 - ev):.4f} km/s') print(f'single impulse = {ev*VC:.4f} km/s') rho_n, sig_n = (R1 + ra_t)/ra_t, (rp_t + R1)/rp_t print(f'\nroute (i): r1 + ra = {R1 + ra_t:.1f} km, rho = (r1+ra)/ra = {rho_n:.4f}, sqrt(2/rho) = {np.sqrt(2/rho_n):.5f}') h1n = np.sqrt(2*MU*R1*ra_t/(R1 + ra_t)) print(f' h = {h1n:.1f} km^2/s; speed at r1 = {h1n/R1:.4f} km/s; speed at ra = {h1n/ra_t:.4f} km/s; target at ra = {R1*VC/ra_t:.4f} km/s') A1n = VC*(np.sqrt(2/rho_n) - 1) print(f' burn A = {A1n:.4f} km/s = {A1n/VC:.5f} vc; burn B = (rho-1) A = {(rho_n - 1)*A1n:.4f} km/s = {(rho_n - 1)*A1n/VC:.5f} vc;' f' total = {VC*(np.sqrt(2*rho_n) - rho_n):.4f} km/s = {np.sqrt(2*rho_n) - rho_n:.5f} vc') print(f'route (ii): rp + r1 = {rp_t + R1:.1f} km, sigma = (rp+r1)/rp = {sig_n:.4f}, sqrt(2/sigma) = {np.sqrt(2/sig_n):.5f}') h2n = np.sqrt(2*MU*rp_t*R1/(rp_t + R1)) print(f' h = {h2n:.1f} km^2/s; speed at r1 = {h2n/R1:.4f} km/s; speed at rp = {h2n/rp_t:.4f} km/s; target at rp = {R1*VC/rp_t:.4f} km/s') A2n = VC*(1 - np.sqrt(2/sig_n)) print(f' burn A = {A2n:.4f} km/s = {A2n/VC:.5f} vc; burn B = (sigma-1) A = {(sig_n - 1)*A2n:.4f} km/s = {(sig_n - 1)*A2n/VC:.5f} vc;' f' total = {VC*(sig_n - np.sqrt(2*sig_n)):.4f} km/s = {sig_n - np.sqrt(2*sig_n):.5f} vc') print(f'\nshapes (target: e = {ev}, p = r1 = {R1:.0f} km):') print(f' T1: e1 = (ra - r1)/(ra + r1) = {(ra_t - R1)/(ra_t + R1):.4f} = e/(2-e) = {ev/(2 - ev):.4f}; p1 = 2 r1/rho = {2*R1/rho_n:.1f} km') print(f' T2: e2 = (r1 - rp)/(r1 + rp) = {(R1 - rp_t)/(R1 + rp_t):.4f} = e/(2+e) = {ev/(2 + ev):.4f}; p2 = 2 r1/sigma = {2*R1/sig_n:.1f} km') print(f' target: e = (ra - rp)/(ra + rp) = {(ra_t - rp_t)/(ra_t + rp_t):.4f}; (ra - rp, ra + rp) = ({ra_t - rp_t:.1f}, {ra_t + rp_t:.1f}) km') (A1, B1), (A2, B2), t1, t2 = route_costs(ev) c1, c2 = A1 + B1, A2 + B2 print(f'\nsaving vs single: (i) {100*(1 - c1/(ev*VC)):.1f} %, (ii) {100*(1 - c2/(ev*VC)):.1f} %') print(f'coast: (i) {t1/60:.1f} min, (ii) {t2/60:.1f} min') print(f'transfer semi-major axes: a1 = {(R1 + ra_t)/2:.1f} km, a2 = {(rp_t + R1)/2:.1f} km') print('\n e single (i) (ii) (ii)/(i) saving(i) saving(ii) [costs in units of vc]') for ev in (0.1, 0.3, 0.5, 0.7, 0.9, 0.99): d1 = np.sqrt(2*(2 - ev)) - (2 - ev) d2 = (2 + ev) - np.sqrt(2*(2 + ev)) print(f'{ev:5.2f} {ev:7.4f} {d1:8.5f} {d2:8.5f} {d2/d1:8.4f} {100*(1 - d1/ev):6.1f} % {100*(1 - d2/ev):6.1f} %') print(f'e -> 1: (i) {np.sqrt(2) - 1:.4f}, (ii) {3 - np.sqrt(6):.4f}') # illustrative propellant fractions, Isp = 300 s (the Isp is an assumption, not given in the paper) g0, Isp, ev = 9.81e-3, 300.0, 0.3 ve = g0*Isp d1 = np.sqrt(2*(2 - ev)) - (2 - ev) d2 = (2 + ev) - np.sqrt(2*(2 + ev)) print(f'\nIsp = {Isp:.0f} s: g0 Isp = {ve:.3f} km/s') for name, d in (('single', ev), ('(i)', d1), ('(ii)', d2)): print(f' {name:7s} dv = {d*VC:.4f} km/s propellant fraction = {100*(1 - np.exp(-d*VC/ve)):.1f} %')