"""Verification of every number and identity quoted in the solutions of midsem-2025.qmd (Mid-Semester, 9 October 2025). Q1 (c) Mercury's speed; (l) the two-body errors of Jupiter and Saturn; (m) the shell integral; (p) a from rp and ra Q2 the two possible chords of the parabola Q3 the Hohmann transfer from GSO to LEO with the paper's own data Q4 the velocity in the perifocal frame Q5 Kepler's equation from the area integral Run: python3 courses/sfm/exam/midsem-2025-verify.py """ import numpy as np import sympy as sp from scipy.integrate import quad fails = 0 def ok(name, x, y, tol): global fails good = abs(float(x) - float(y)) <= tol fails += not good print("ok " if good else "FAIL", name, round(float(x), 6), y) # data of the paper G, AU, c = 6.6743e-11, 1.496e11, 3e8 mJ, mS, aJ, aS, aMercury = 1.8982e27, 5.6834e26, 5.2026, 9.5549, 0.387098 mu_E, R_E = 3.986e5, 6371.0 # km^3/s^2, km year = 365.25*86400 # not in the table, general knowledge # Q1 (c): the Sun's mu from the Earth's orbit, which needs only the year mu_sun = 4*np.pi**2*AU**3/year**2 M_sun = mu_sun/G ok("mu_sun, 1e20 m^3/s^2", mu_sun/1e20, 1.327, 5e-4) ok("M_sun, 1e30 kg", M_sun/1e30, 1.989, 5e-4) v_mercury = np.sqrt(mu_sun/(aMercury*AU)) ok("Mercury, km/s", v_mercury/1e3, 47.9, 0.05) ok("Earth's speed, km/s", np.sqrt(mu_sun/AU)/1e3, 29.8, 0.05) ok("Mercury from the Earth's speed", 29.78/np.sqrt(aMercury), 47.9, 0.05) ok("v/c, 1e-4", v_mercury/c/1e-4, 1.6, 0.005) ok("fraction of 0.1 c, percent", v_mercury/(0.1*c)*100, 0.16, 0.005) # Q1 (l): the pull of the other planet at conjunction, as a fraction of the Sun's pull d = aS - aJ ok("separation, AU", d, 4.3523, 5e-5) eJ = mS/M_sun*(aJ/d)**2*100 eS = mJ/M_sun*(aS/d)**2*100 ok("m_S/M_sun, 1e-4", mS/M_sun/1e-4, 2.858, 5e-4) ok("m_J/M_sun, 1e-4", mJ/M_sun/1e-4, 9.546, 5e-4) ok("(a_J/d)^2", (aJ/d)**2, 1.429, 5e-4) ok("(a_S/d)^2", (aS/d)**2, 4.820, 5e-4) ok("error for Jupiter, percent", eJ, 0.041, 5e-4) ok("error for Saturn, percent", eS, 0.46, 5e-3) ok("Saturn's error, as 1 in", 100/eS, 217, 0.5) ok("ratio of the two errors", eS/eJ, 11, 0.5) # Q1 (m): the shell integral, numerically and for a sample density for D, s in ((2.0, 1.0), (5.0, 0.3), (1.2, 1.0)): I, _ = quad(lambda t: 0.5*np.sin(t)/np.sqrt(s*s + D*D - 2*s*D*np.cos(t)), 0, np.pi) ok("mean of 1/l over a shell = 1/D", I, 1/D, 1e-10) s_, D_, ell = sp.symbols("s D ell", positive=True) ok("the integral in l", sp.simplify(sp.integrate(1, (ell, D_ - s_, D_ + s_))/(2*s_*D_) - 1/D_), 0, 0) # a sphere with rho = 1 + s/R: the sum over shells gives -G M/D R = 1.0 rho = lambda s: 1 + s/R M = quad(lambda s: 4*np.pi*s*s*rho(s), 0, R)[0] D = 1.7 pot, _ = quad(lambda s: 4*np.pi*s*s*rho(s)*quad(lambda t: 0.5*np.sin(t)/np.sqrt(s*s + D*D - 2*s*D*np.cos(t)), 0, np.pi)[0], 0, R) ok("whole sphere = M/D", pot, M/D, 1e-8) # Q1 (p) h2mu, e_ = sp.symbols("p e", positive=True) rp, ra = h2mu/(1 + e_), h2mu/(1 - e_) ok("a = p/(1 - e^2)", sp.simplify((rp + ra)/2 - h2mu/(1 - e_**2)), 0, 0) ok("c = a e = (ra - rp)/2", sp.simplify((ra - rp)/2 - e_*h2mu/(1 - e_**2)), 0, 0) # Q2: the parabola q, r1_, r2_ = 7000.0, 8000.0, 16000.0 p = 2*q ok("p", p, 14000, 0) c1, c2 = p/r1_ - 1, p/r2_ - 1 ok("cos f1", c1, 0.75, 1e-12); ok("cos f2", c2, -0.125, 1e-12) f1, f2 = np.arccos(c1), np.arccos(c2) ok("f1, deg", np.degrees(f1), 41.41, 5e-3); ok("f2, deg", np.degrees(f2), 97.18, 5e-3) ok("sin f1", np.sin(f1), 0.6614, 5e-5); ok("sin f2", np.sin(f2), 0.9922, 5e-5) ok("cos(f2 - f1)", np.cos(f2 - f1), 0.5625, 1e-12); ok("cos(f2 + f1)", np.cos(f2 + f1), -0.75, 1e-12) ok("f2 - f1", np.degrees(f2 - f1), 55.77, 5e-3); ok("f2 + f1", np.degrees(f2 + f1), 138.59, 5e-3) same = np.sqrt(r1_**2 + r2_**2 - 2*r1_*r2_*np.cos(f2 - f1)); opp = np.sqrt(r1_**2 + r2_**2 - 2*r1_*r2_*np.cos(f1 + f2)) ok("d, same side", same, 13266, 0.5); ok("d, opposite sides", opp, 22627, 0.5) ok("d^2, same side", same**2, 1.76e8, 1e-3); ok("d^2, opposite sides", opp**2, 5.12e8, 1e-3) ok("cartesian, same side", np.hypot(r2_*np.cos(f2) - r1_*np.cos(f1), r2_*np.sin(f2) - r1_*np.sin(f1)), same, 1e-8) ok("cartesian, opposite", np.hypot(r2_*np.cos(f2) - r1_*np.cos(-f1), r2_*np.sin(f2) - r1_*np.sin(-f1)), opp, 1e-8) ok("the closed forms", 8000*np.sqrt(2.75), same, 1e-8); ok("the closed forms", 16000*np.sqrt(2), opp, 1e-8) # Q3: GSO to LEO with the paper's data r2, r1 = R_E + 35786, R_E + 300 ok("r_GSO", r2, 42157, 0); ok("r_LEO", r1, 6671, 0) at = (r1 + r2)/2 ok("a_t", at, 24414, 0); ok("e_t", (r2 - r1)/(r2 + r1), 0.7268, 5e-5); ok("R", r2/r1, 6.32, 5e-3) vc2, vc1 = np.sqrt(mu_E/r2), np.sqrt(mu_E/r1) vt2, vt1 = np.sqrt(mu_E*(2/r2 - 1/at)), np.sqrt(mu_E*(2/r1 - 1/at)) ok("v_c,GSO", vc2, 3.0749, 5e-5); ok("v_t at GSO", vt2, 1.6073, 5e-5); ok("v_t at LEO", vt1, 10.1575, 5e-5); ok("v_c,LEO", vc1, 7.7299, 5e-5) ok("first burn", vc2 - vt2, 1.4676, 5e-5); ok("second burn", vt1 - vc1, 2.4277, 5e-5) ok("total", vc2 - vt2 + vt1 - vc1, 3.8952, 5e-5) ok("as printed", 1.4676 + 2.4277, 3.8953, 1e-9) ok("time, s", np.pi*np.sqrt(at**3/mu_E), 18982, 0.5); ok("time, h", np.pi*np.sqrt(at**3/mu_E)/3600, 5.27, 5e-3) ok("mass left, I_sp = 315 s, percent", 100*np.exp(-(vc2 - vt2 + vt1 - vc1)/(315*9.81e-3)), 28.4, 0.05) # Q4: v from differentiating r = r (cos f, sin f), with r^2 fdot = h f, h, mu, e = sp.symbols("f h mu e", positive=True) rr = (h**2/mu)/(1 + e*sp.cos(f)) fdot = h/rr**2 vx = sp.simplify(sp.diff(rr*sp.cos(f), f)*fdot) vy = sp.simplify(sp.diff(rr*sp.sin(f), f)*fdot) ok("v_x", sp.simplify(vx + mu/h*sp.sin(f)), 0, 0) ok("v_y", sp.simplify(vy - mu/h*(e + sp.cos(f))), 0, 0) ok("rdot", sp.simplify(sp.diff(rr, f)*fdot - mu/h*e*sp.sin(f)), 0, 0) # Q5: x dy - y dx, then the area E_, a, b = sp.symbols("E a b", positive=True) x, y = a*(sp.cos(E_) - e), b*sp.sin(E_) ok("x y' - y x'", sp.simplify(x*sp.diff(y, E_) - y*sp.diff(x, E_) - a*b*(1 - e*sp.cos(E_))), 0, 0) ok("A(2 pi) = pi a b", sp.simplify(sp.integrate(a*b*(1 - e*sp.cos(E_))/2, (E_, 0, 2*sp.pi)) - sp.pi*a*b), 0, 0) print("\nfailures:", fails) raise SystemExit(1 if fails else 0)