MATH307Su26 - Orbital Mechanics and Numerical Integration

# Orbital Mechanics and Numerical Integration > [!abstract] The project in one line > Turn Kepler's second law into a definite integral for orbital transit time, evaluate it with the trapezoid and Simpson rules, verify both the answers and the convergence rates against an exact ground truth, then run the same instrument on orbits of your own choosing. ## Introduction You are tracking the International Space Station as it orbits Earth, watching its position update every few seconds on a real-time map. The station travels at over 7.6 km/s, yet its orbital period is precisely predictable (about 93 minutes per lap). How do mission controllers calculate exactly when the ISS will pass over a specific location, or how long a communication satellite lingers over the far side of its orbit? The mathematics is one of the cleanest applications of integration in this course. When [Johannes Kepler](https://en.wikipedia.org/wiki/Johannes_Kepler) discovered that planets sweep out equal areas in equal times, he described (without knowing it) an integration problem: the time to travel an arc of orbit is an integral of the swept area. For elliptical orbits the integrand varies wildly between the near and far points, and evaluating the integral is a genuine computation, not a table lookup. ![Same angular arc, very different transit times](Media/orbital_geometry.png) ## Project Description The project moves in four passes. First, derive the transit-time integral from Kepler's second law by hand and run the warm-ups to checkable numbers (the ISS period, the integrand's swing, the comet ratio). Second, mirror the hand work in code: implement the trapezoid rule from scratch and confirm the circular case to machine precision through the verification helper. Third, run the real test (Simpson's rule and a convergence study against a Kepler-equation exact value) and generalize to the Molniya orbit, where Kepler's second law bites hardest. Finally, run the whole study on three real orbits you choose, picked for three different eccentricity personalities. **Foundation task:** Derive the transit-time integral from Kepler's second law, implement both quadrature rules from scratch, confirm the circular-orbit case to machine precision, and run a convergence study on an elliptical arc whose exact transit time is known. (This is Steps 1 to 5 below.) **Application task:** Use the verified machinery on the Molniya orbit to quantify how a satellite crawls near apogee and sprints near perigee, connect the numbers to orbital speeds, and investigate how eccentricity drives the cost of accurate quadrature. (This is the Generalization section below.) **Key deliverable:** The same study on three real orbits you choose with three eccentricity personalities (near-circular, working eccentric, and extreme): the difficulty dial read and interpreted before any computing, each transit time verified against its Kepler-equation exact value, and a cross-orbit cost ranking defended from the convergence numbers. Two companions support the milestone forms: the [[MATH307Su26 - Orbital Mechanics and Numerical Integration (Milestone Map)|Milestone Map]] says which artifact from this handout answers each form field, and the [[MATH307Su26 - Orbital Mechanics and Numerical Integration (Verification Table)|Verification Table]] is filled in as you verify and submitted with Milestone 2. ## Mathematical Background: from swept area to a definite integral **Idea.** [Kepler's second law](https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion) plus the polar area element $dA=\tfrac12 r^2\,d\theta$ from Calc 3 turns "how long does this arc take" into a definite integral. We evaluate it with the [trapezoid rule](https://en.wikipedia.org/wiki/Trapezoidal_rule) and [Simpson's rule](https://en.wikipedia.org/wiki/Simpson%27s_rule) from class, and this problem comes with a rare gift: an exact answer (via [Kepler's equation](https://en.wikipedia.org/wiki/Kepler%27s_equation)) to test our numerics against. **Recall.** In polar coordinates the area of a thin sector of opening $d\theta$ at radius $r$ is $dA=\tfrac12 r^2\,d\theta$ (a triangle of base $r\,d\theta$ and height $r$, in the limit). This is the Calc 3 area element, and it is the whole geometric input. **Kepler's second law.** An orbiting body sweeps out equal areas in equal times: the areal velocity is constant, $\frac{dA}{dt}=\frac{h}{2},$ where $h$ is the [specific angular momentum](https://en.wikipedia.org/wiki/Specific_angular_momentum) of the orbit, a conserved quantity with $h=\sqrt{GM\,a(1-e^2)}$ for an ellipse with semi-major axis $a$ and [eccentricity](https://en.wikipedia.org/wiki/Orbital_eccentricity) $e$ around a body of gravitational parameter $GM$. Setting the two expressions for $dA$ equal, $\tfrac12 r^2\,d\theta=\tfrac{h}{2}\,dt \quad\Longrightarrow\quad dt=\frac{r^2}{h}\,d\theta.$ **The orbit shape.** For an ellipse with the attracting body at a focus, the radius at [true anomaly](https://en.wikipedia.org/wiki/True_anomaly) $\theta$ (the angle measured from closest approach) is $r(\theta)=\frac{a(1-e^2)}{1+e\cos\theta}.$ **The transit-time integral.** Integrating $dt$ from $\theta_1$ to $\theta_2$ gives the time to traverse that arc, $\boxed{\ t(\theta_1,\theta_2)=\frac{1}{h}\int_{\theta_1}^{\theta_2} r(\theta)^2\,d\theta =\frac{a^2(1-e^2)^2}{h}\int_{\theta_1}^{\theta_2}\frac{d\theta}{\bigl(1+e\cos\theta\bigr)^2}.\ }$ For $e>0$ the integrand swings between $1/(1+e)^2$ at perigee and $1/(1-e)^2$ at apogee. At $e=0.5$ that is $\tfrac49\approx0.444$ versus $4$, a factor of $9$; at the Molniya value $e=0.74$ it is a factor of $\bigl(\tfrac{1.74}{0.26}\bigr)^2\approx44.8$. Rapidly varying integrands are exactly where quadrature rules earn their keep. > [!note] Where this comes from > Nothing here is new machinery: the area element is Calc 3, the conservation law is supplied physics, and the quadrature rules and their $O(h^2)$/$O(h^4)$ error orders are the Taylor-series analysis from class. The project is the assembly. ### Worked example: the circular orbit (your first ground truth) For $e=0$ the radius is constant, $r=a$, and $h=\sqrt{GM\,a}$, so the integral collapses: $t(\theta_1,\theta_2)=\frac{a^2}{h}(\theta_2-\theta_1)=\sqrt{\frac{a^3}{GM}}\,(\theta_2-\theta_1), \qquad T=2\pi\sqrt{\frac{a^3}{GM}},$ which is [Kepler's third law](https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion#Third_law). For the ISS (altitude $408$ km above Earth's mean radius $6371$ km, so $a=6779$ km, with $GM=398{,}600\ \text{km}^3/\text{s}^2$): $T=2\pi\sqrt{\frac{6779^3}{398600}}=5554.7\ \text{s}\approx 92.6\ \text{min},\qquad h=\sqrt{GM\,a}=51{,}981.8\ \text{km}^2/\text{s},$ and a quarter orbit takes $T/4=1388.67$ s. Confirm these by hand before running anything. **Hand warm-ups (do these before coding; the warm-up code step reproduces every one).** 1. Earth around the Sun ($a=149.6$ million km, $e=0.017$): perihelion $a(1-e)=147.06$ and aphelion $a(1+e)=152.14$ million km, a $3.46\%$ swing. 2. A comet with $e=0.9$: $r_{\max}/r_{\min}=(1+e)/(1-e)=19$. 3. The integrand $1/(1+e\cos\theta)^2$ at $\theta=0,\ \pi/2,\ \pi$ for $e=0.5$: $\tfrac49,\ 1,\ 4$. > [!note] An exact answer exists for the ellipse too (and we will use it) > The substitution to the [eccentric anomaly](https://en.wikipedia.org/wiki/Eccentric_anomaly) $E$ evaluates the transit integral exactly: with $\tan\frac{E}{2}=\sqrt{\frac{1-e}{1+e}}\,\tan\frac{\theta}{2}$, the time from perigee is > $t(\theta)=\sqrt{\frac{a^3}{GM}}\,\bigl(E-e\sin E\bigr),$ > which is [Kepler's equation](https://en.wikipedia.org/wiki/Kepler%27s_equation). So why compute numerically at all? First, the *inverse* problem (where is the satellite at time $t$) is transcendental and needs numerics anyway. Second, the substitution is special to this integrand; add drag, thrust, or a nonspherical Earth and it dies, while quadrature does not care. Here the formula is our gift: an exact ground truth to hold the numerics against. ## Implementation Guidelines The workflow is the same in every language, so we state it once and then implement it step by step. > [!abstract] The workflow (state it once, run it in any language) > Given an orbit $(a,e)$ around a body of gravitational parameter $GM$, an arc $[\theta_1,\theta_2]$, and an interval count $n$: > 1. **Warm up**: define the verification helper, then reproduce the hand numbers (the integrand values $\tfrac49, 1, 4$ at $e=0.5$; Earth's perihelion and aphelion; the comet ratio $19$; the ISS period $T=5554.7$ s and quarter $T/4=1388.672$ s). > 2. **Define and view**: code $r(\theta)=a(1-e^2)/(1+e\cos\theta)$ and $h=\sqrt{GM\,a(1-e^2)}$, and draw the orbit with the body at the focus. > 3. **Trapezoid**: $t\approx\Delta\theta\,\bigl(\tfrac{f_0}{2}+f_1+\cdots+f_{n-1}+\tfrac{f_n}{2}\bigr)$ with $f=r^2/h$; verify the circular case to machine precision. > 4. **Simpson and converge**: $t\approx\tfrac{\Delta\theta}{3}\bigl(f_0+4f_{\text{odd}}+2f_{\text{even}}+f_n\bigr)$ ($n$ even); against the Kepler-equation exact value, doubling $n$ must cut the errors by $4$ and $16$. > 5. **Render**: the log-log convergence plot with $O(h^2)$ and $O(h^4)$ reference slopes. > 6. **Generalize**: perigee versus apogee dwell on the Molniya orbit, then eccentricity as the difficulty dial (its own section below). > 7. **Your own orbits**: three eccentricity personalities, prediction first (the finale). Pick your language below; each callout contains Steps 1 to 5 in order, run top to bottom. Code is four-language and portable (Python is the run-verified reference; Mathematica is numerical only). This project ships no data files: every orbit is a pair of numbers, so there is nothing to download. > [!example]- Python (reference) > > **Step 1: define the helper, then warm up on the hand numbers.** The helper compares two values (or arrays) and reports the largest difference; every check in this project runs through it. > ```python > import numpy as np > import matplotlib.pyplot as plt > def same(name, X, Y, tol=1e-9): # verification helper > d = np.max(np.abs(np.asarray(X) - np.asarray(Y))) > print(f"{name}: max difference {d:.1e}", "OK" if d < tol else "CHECK THIS") > GM = 398600.0 # km^3/s^2 (Earth) > r_of = lambda th, a, e: a*(1 - e**2)/(1 + e*np.cos(th)) > h_of = lambda a, e: np.sqrt(GM*a*(1 - e**2)) > vals = [1/(1 + 0.5*np.cos(t))**2 for t in (0.0, np.pi/2, np.pi)] > same("integrand values equal 4/9, 1, 4 (e = 0.5)", vals, [4/9, 1.0, 4.0]) > print("Earth: perihelion", 149.6*(1 - 0.017), " aphelion", 149.6*(1 + 0.017)) # 147.06 152.14 > same("comet r_max/r_min at e = 0.9 equals 19", (1 + 0.9)/(1 - 0.9), 19.0) > a_iss, e_iss = 6371.0 + 408.0, 0.0 # ISS: a = 6779 km, circular > T_iss = 2*np.pi*np.sqrt(a_iss**3/GM) # Kepler III > print("ISS: T =", round(T_iss, 1), " quarter =", round(T_iss/4, 3), " h =", round(h_of(a_iss, e_iss), 1)) > ``` > *Look for:* the hand numbers exactly: integrand check at $0$, perihelion/aphelion $147.06/152.14$, and ISS $T=5554.7$ s, $T/4=1388.672$ s, $h=51981.8$ km²/s. > > **Step 2: define the orbit and see it.** Draw the Molniya orbit ($a=26{,}560$ km, $e=0.74$) with Earth at the focus; the ellipse is visibly off-center. > ```python > a_m, e_m = 26560.0, 0.74 # Molniya orbit > th = np.linspace(0, 2*np.pi, 800) > r = r_of(th, a_m, e_m) > plt.plot(r*np.cos(th), r*np.sin(th)); plt.plot(0, 0, "o") # Earth at the focus > plt.gca().set_aspect("equal"); plt.title("Molniya orbit"); plt.show() > print(r_of(0, a_m, e_m), r_of(np.pi, a_m, e_m)) # 6905.6 46214.4 > ``` > *Look for:* closest approach $6905.6$ km, farthest $46{,}214.4$ km; Earth sits at the focus, not the ellipse's center, so the orbit is visibly lopsided. > > **Step 3: trapezoid rule, verified on the circular orbit.** For $e=0$ the integrand $r^2/h$ is constant, so the trapezoid rule commits no error at all: this run checks correctness (grid, endpoints, weights), not convergence. > ```python > def transit_trap(t1, t2, a, e, n): > th = np.linspace(t1, t2, n+1); dth = (t2 - t1)/n > y = r_of(th, a, e)**2 / h_of(a, e) > return dth*(0.5*y[0] + y[1:-1].sum() + 0.5*y[-1]) > t_num = transit_trap(0, np.pi/2, a_iss, e_iss, 100) > print(t_num, T_iss/4) # 1388.672 1388.672 > same("circular quarter orbit at machine precision", t_num, T_iss/4) > ``` > *Look for:* the error must be at machine level ($\sim10^{-13}$), not merely small. If it is $10^{-4}$, a weight or endpoint is wrong. > > **Step 4: Simpson's rule and the convergence study.** On the Molniya quarter arc $[0,\pi/2]$ the integrand varies strongly, and Kepler's equation gives the exact transit time $t^\*=1645.4287$ s. Doubling $n$ must cut the trapezoid error by $4$ and the Simpson error by $16$. > ```python > def transit_simp(t1, t2, a, e, n): # n must be even > th = np.linspace(t1, t2, n+1); dth = (t2 - t1)/n > y = r_of(th, a, e)**2 / h_of(a, e) > return dth/3*(y[0] + 4*y[1::2].sum() + 2*y[2:-1:2].sum() + y[-1]) > ecc = 2*np.arctan(np.sqrt((1 - e_m)/(1 + e_m))*np.tan(np.pi/4)) # eccentric anomaly at theta = pi/2 > t_star = np.sqrt(a_m**3/GM)*(ecc - e_m*np.sin(ecc)) # exact: 1645.4287 s > for n in (10, 20, 40, 80, 160, 320): > et = abs(transit_trap(0, np.pi/2, a_m, e_m, n) - t_star) > es = abs(transit_simp(0, np.pi/2, a_m, e_m, n) - t_star) > print(n, f"{et:.3e}", f"{es:.3e}") # errors drop 4x and 16x per doubling > T_m = 2*np.pi*np.sqrt(a_m**3/GM) # full Molniya period, 43077.8 s = 11.966 h > for n in (16, 32, 64): # the periodic surprise (see the note below) > print(n, f"{abs(transit_trap(0, 2*np.pi, a_m, e_m, n) - T_m):.1e}") > ``` > *Look for:* at $n=320$ the trapezoid error is $6.2\times10^{-3}$ s and the Simpson error $5.5\times10^{-8}$ s; the full-orbit trapezoid errors collapse ($2.2$ s, $9\times10^{-6}$ s, $7\times10^{-12}$ s at $n=16, 32, 64$), far faster than $O(h^2)$. > > **Step 5: plot the convergence.** Straight lines parallel to the $O(h^2)$ and $O(h^4)$ references are the visual proof that the implementation meets the theory. > ```python > ns = np.array([10, 20, 40, 80, 160, 320]) > et = [abs(transit_trap(0, np.pi/2, a_m, e_m, n) - t_star) for n in ns] > es = [abs(transit_simp(0, np.pi/2, a_m, e_m, n) - t_star) for n in ns] > plt.loglog(ns, et, "o-", label="trapezoid"); plt.loglog(ns, es, "s-", label="Simpson") > plt.loglog(ns, et[0]*(ns[0]/ns)**2.0, "--", color="0.6", label="$O(h^2)quot;) > plt.loglog(ns, es[0]*(ns[0]/ns)**4.0, ":", color="0.6", label="$O(h^4)quot;) > plt.xlabel("n"); plt.ylabel("absolute error (s)"); plt.legend(); plt.show() > ``` > *Look for:* both curves straight on log-log axes and parallel to their reference slopes. > [!example]- MATLAB > > **Step 1: define the helper, then warm up on the hand numbers.** The helper reports the largest difference; every check runs through it. > ```matlab > same = @(name, X, Y) fprintf('%s: max difference %.1e\n', name, max(abs(X - Y), [], 'all')); > GM = 398600; % km^3/s^2 (Earth) > r_of = @(th, a, e) a*(1 - e^2) ./ (1 + e*cos(th)); > h_of = @(a, e) sqrt(GM*a*(1 - e^2)); > vals = 1 ./ (1 + 0.5*cos([0 pi/2 pi])).^2; > same('integrand values equal 4/9, 1, 4 (e = 0.5)', vals, [4/9 1 4]) > fprintf('Earth: perihelion %.2f aphelion %.2f\n', 149.6*(1 - 0.017), 149.6*(1 + 0.017)) > same('comet r_max/r_min at e = 0.9 equals 19', (1 + 0.9)/(1 - 0.9), 19) > a_iss = 6371 + 408; e_iss = 0; % ISS: a = 6779 km, circular > T_iss = 2*pi*sqrt(a_iss^3/GM); % Kepler III > fprintf('ISS: T = %.1f quarter = %.3f h = %.1f\n', T_iss, T_iss/4, h_of(a_iss, e_iss)) > ``` > *Look for:* $147.06/152.14$; ISS $T=5554.7$ s, $T/4=1388.672$ s, $h=51981.8$ km²/s. > > **Step 2: define the orbit and see it.** > ```matlab > a_m = 26560; e_m = 0.74; % Molniya orbit > th = linspace(0, 2*pi, 800); r = r_of(th, a_m, e_m); > plot(r.*cos(th), r.*sin(th)); hold on; plot(0, 0, 'o'); axis equal; title('Molniya orbit') > fprintf('%.1f %.1f\n', r_of(0, a_m, e_m), r_of(pi, a_m, e_m)) % 6905.6 46214.4 > ``` > *Look for:* the division must be elementwise (`./`); Earth at the focus, visibly lopsided. > > **Step 3: trapezoid rule, verified on the circular orbit.** > ```matlab > n = 100; t1 = 0; t2 = pi/2; > th = linspace(t1, t2, n+1); dth = (t2 - t1)/n; > y = r_of(th, a_iss, e_iss).^2 / h_of(a_iss, e_iss); > t_num = dth*(0.5*y(1) + sum(y(2:end-1)) + 0.5*y(end)); > fprintf('%.3f %.3f\n', t_num, T_iss/4) > same('circular quarter orbit at machine precision', t_num, T_iss/4) > ``` > *Look for:* machine-level error ($\sim10^{-13}$), not merely small. > > **Step 4: Simpson's rule and the convergence study.** > ```matlab > ecc = 2*atan(sqrt((1 - e_m)/(1 + e_m))*tan(pi/4)); % eccentric anomaly at theta = pi/2 > t_star = sqrt(a_m^3/GM)*(ecc - e_m*sin(ecc)); % exact: 1645.4287 s > for n = [10 20 40 80 160 320] > th = linspace(0, pi/2, n+1); dth = (pi/2)/n; > y = r_of(th, a_m, e_m).^2 / h_of(a_m, e_m); > t_tr = dth*(0.5*y(1) + sum(y(2:end-1)) + 0.5*y(end)); > t_si = dth/3*(y(1) + 4*sum(y(2:2:end-1)) + 2*sum(y(3:2:end-2)) + y(end)); > fprintf('%4d %.3e %.3e\n', n, abs(t_tr - t_star), abs(t_si - t_star)); > end > T_m = 2*pi*sqrt(a_m^3/GM); % full period, 43077.8 s = 11.966 h > for n = [16 32 64] % the periodic surprise (see the note below) > th = linspace(0, 2*pi, n+1); dth = 2*pi/n; > y = r_of(th, a_m, e_m).^2 / h_of(a_m, e_m); > fprintf('%3d %.1e\n', n, abs(dth*(0.5*y(1) + sum(y(2:end-1)) + 0.5*y(end)) - T_m)); > end > ``` > *Look for:* the odd-index (weight $4$) and even-index (weight $2$) slices are 1-based here, unlike Python; errors at $n=320$: $6.2\times10^{-3}$ and $5.5\times10^{-8}$ s; full-orbit errors collapse. > > **Step 5: plot the convergence.** > ```matlab > ns = [10 20 40 80 160 320]; et = zeros(size(ns)); es = et; > for i = 1:numel(ns) % reuse the Step 4 loop body > n = ns(i); th = linspace(0, pi/2, n+1); dth = (pi/2)/n; > y = r_of(th, a_m, e_m).^2 / h_of(a_m, e_m); > et(i) = abs(dth*(0.5*y(1) + sum(y(2:end-1)) + 0.5*y(end)) - t_star); > es(i) = abs(dth/3*(y(1) + 4*sum(y(2:2:end-1)) + 2*sum(y(3:2:end-2)) + y(end)) - t_star); > end > figure; loglog(ns, et, 'o-'); hold on; loglog(ns, es, 's-') > loglog(ns, et(1)*(ns(1)./ns).^2, '--', 'Color', [.6 .6 .6]) > loglog(ns, es(1)*(ns(1)./ns).^4, ':', 'Color', [.6 .6 .6]) > xlabel('n'); ylabel('absolute error (s)') > legend('trapezoid', 'Simpson', 'O(h^2)', 'O(h^4)') > ``` > *Look for:* straight parallel lines on log-log axes. > [!example]- R > > **Step 1: define the helper, then warm up on the hand numbers.** The helper reports the largest difference; every check runs through it. > ```r > same <- function(name, X, Y) cat(sprintf("%s: max difference %.1e\n", name, max(abs(X - Y)))) > GM <- 398600 # km^3/s^2 (Earth) > r_of <- function(th, a, e) a*(1 - e^2) / (1 + e*cos(th)) > h_of <- function(a, e) sqrt(GM*a*(1 - e^2)) > vals <- 1 / (1 + 0.5*cos(c(0, pi/2, pi)))^2 > same("integrand values equal 4/9, 1, 4 (e = 0.5)", vals, c(4/9, 1, 4)) > cat("Earth: perihelion", 149.6*(1 - 0.017), " aphelion", 149.6*(1 + 0.017), "\n") > same("comet r_max/r_min at e = 0.9 equals 19", (1 + 0.9)/(1 - 0.9), 19) > a_iss <- 6371 + 408; e_iss <- 0 # ISS: a = 6779 km, circular > T_iss <- 2*pi*sqrt(a_iss^3/GM) # Kepler III > cat("ISS: T =", T_iss, " quarter =", T_iss/4, " h =", h_of(a_iss, e_iss), "\n") > ``` > *Look for:* $147.06/152.14$; ISS $T=5554.7$ s, $T/4=1388.672$ s, $h=51981.8$ km²/s. > > **Step 2: define the orbit and see it.** > ```r > a_m <- 26560; e_m <- 0.74 # Molniya orbit > th <- seq(0, 2*pi, length.out = 800); r <- r_of(th, a_m, e_m) > plot(r*cos(th), r*sin(th), type = "l", asp = 1, main = "Molniya orbit") > points(0, 0, pch = 19) > cat(r_of(0, a_m, e_m), r_of(pi, a_m, e_m), "\n") # 6905.6 46214.4 > ``` > *Look for:* Earth at the focus, visibly lopsided ellipse. > > **Step 3: trapezoid rule, verified on the circular orbit.** > ```r > transit_trap <- function(t1, t2, a, e, n) { > th <- seq(t1, t2, length.out = n+1); dth <- (t2 - t1)/n > y <- r_of(th, a, e)^2 / h_of(a, e) > dth*(0.5*y[1] + sum(y[2:n]) + 0.5*y[n+1]) > } > t_num <- transit_trap(0, pi/2, a_iss, e_iss, 100) > cat(t_num, T_iss/4, "\n") # 1388.672 1388.672 > same("circular quarter orbit at machine precision", t_num, T_iss/4) > ``` > *Look for:* machine-level error ($\sim10^{-13}$), not merely small. > > **Step 4: Simpson's rule and the convergence study.** > ```r > transit_simp <- function(t1, t2, a, e, n) { # n must be even > th <- seq(t1, t2, length.out = n+1); dth <- (t2 - t1)/n > y <- r_of(th, a, e)^2 / h_of(a, e) > dth/3*(y[1] + 4*sum(y[seq(2, n, 2)]) + 2*sum(y[seq(3, n-1, 2)]) + y[n+1]) > } > ecc <- 2*atan(sqrt((1 - e_m)/(1 + e_m))*tan(pi/4)) # eccentric anomaly at theta = pi/2 > t_star <- sqrt(a_m^3/GM)*(ecc - e_m*sin(ecc)) # exact: 1645.4287 s > for (n in c(10, 20, 40, 80, 160, 320)) > cat(n, abs(transit_trap(0, pi/2, a_m, e_m, n) - t_star), > abs(transit_simp(0, pi/2, a_m, e_m, n) - t_star), "\n") > T_m <- 2*pi*sqrt(a_m^3/GM) # full period, 43077.8 s = 11.966 h > for (n in c(16, 32, 64)) # the periodic surprise (see the note below) > cat(n, abs(transit_trap(0, 2*pi, a_m, e_m, n) - T_m), "\n") > ``` > *Look for:* the weight-4 and weight-2 index slices are 1-based here, unlike Python; errors at $n=320$: $6.2\times10^{-3}$ and $5.5\times10^{-8}$ s; full-orbit errors collapse. > > **Step 5: plot the convergence.** > ```r > ns <- c(10, 20, 40, 80, 160, 320) > et <- sapply(ns, function(n) abs(transit_trap(0, pi/2, a_m, e_m, n) - t_star)) > es <- sapply(ns, function(n) abs(transit_simp(0, pi/2, a_m, e_m, n) - t_star)) > plot(ns, et, log = "xy", type = "b", pch = 19, xlab = "n", ylab = "absolute error (s)", > ylim = range(c(et, es))) > lines(ns, es, type = "b", pch = 15, col = "navy") > lines(ns, et[1]*(ns[1]/ns)^2, lty = 2, col = "gray") > lines(ns, es[1]*(ns[1]/ns)^4, lty = 3, col = "gray") > legend("bottomleft", c("trapezoid", "Simpson", "O(h^2)", "O(h^4)"), > lty = c(1, 1, 2, 3), pch = c(19, 15, NA, NA), col = c("black", "navy", "gray", "gray")) > ``` > *Look for:* straight parallel lines on log-log axes. > [!example]- Mathematica (numerical only) > > **Step 1: define the helper, then warm up on the hand numbers.** The helper reports the largest difference; every check runs through it. (`E`, `N`, `D` are protected names, so use `ec`, `gm`, and friends.) > ```wolfram > same[name_, x_, y_] := Print[name, ": max difference ", Max@Abs[x - y]]; > gm = 398600.; (* km^3/s^2; numerical only *) > rOf[th_, a_, ec_] := a (1 - ec^2)/(1 + ec Cos[th]); > hOf[a_, ec_] := Sqrt[gm a (1 - ec^2)]; > vals = 1/(1 + 0.5 Cos[#])^2 & /@ {0., N[Pi/2], N[Pi]}; > same["integrand values equal 4/9, 1, 4 (e = 0.5)", vals, {4./9, 1., 4.}] > Print["Earth: perihelion ", 149.6 (1 - 0.017), " aphelion ", 149.6 (1 + 0.017)] > same["comet r_max/r_min at e = 0.9 equals 19", (1 + 0.9)/(1 - 0.9), 19.] > aIss = 6371. + 408.; eIss = 0.; (* ISS: a = 6779 km, circular *) > tIss = 2 Pi Sqrt[aIss^3/gm]; (* Kepler III *) > Print["ISS: T = ", tIss, " quarter = ", tIss/4, " h = ", hOf[aIss, eIss]] > ``` > *Look for:* $147.06/152.14$; ISS $T=5554.7$ s, $T/4=1388.67$ s, $h=51981.8$ km²/s. > > **Step 2: define the orbit and see it.** > ```wolfram > am = 26560.; em = 0.74; (* Molniya orbit *) > ParametricPlot[rOf[th, am, em] {Cos[th], Sin[th]}, {th, 0, 2 Pi}, > Epilog -> {PointSize[.02], Point[{0, 0}]}, PlotLabel -> "Molniya orbit"] > Print[{rOf[0., am, em], rOf[N[Pi], am, em]}] (* {6905.6, 46214.4} *) > ``` > *Look for:* Earth at the focus, visibly lopsided ellipse. > > **Step 3: trapezoid rule, verified on the circular orbit.** > ```wolfram > transitTrap[t1_, t2_, a_, ec_, n_] := Module[{th, dth, y}, > th = Subdivide[N[t1], N[t2], n]; dth = (t2 - t1)/n; > y = rOf[#, a, ec]^2/hOf[a, ec] & /@ th; > dth (0.5 First[y] + Total[y[[2 ;; -2]]] + 0.5 Last[y])] > tNum = transitTrap[0, Pi/2, aIss, eIss, 100]; > Print[{tNum, tIss/4}] (* 1388.67, 1388.67 *) > same["circular quarter orbit at machine precision", tNum, tIss/4] > ``` > *Look for:* machine-level error ($\sim10^{-13}$), not merely small. > > **Step 4: Simpson's rule and the convergence study.** > ```wolfram > transitSimp[t1_, t2_, a_, ec_, n_] := Module[{th, dth, y}, (* n must be even *) > th = Subdivide[N[t1], N[t2], n]; dth = (t2 - t1)/n; > y = rOf[#, a, ec]^2/hOf[a, ec] & /@ th; > dth/3 (First[y] + 4 Total[y[[2 ;; -2 ;; 2]]] + 2 Total[y[[3 ;; -3 ;; 2]]] + Last[y])] > ecc = N[2 ArcTan[Sqrt[(1 - em)/(1 + em)] Tan[Pi/4]]]; (* eccentric anomaly *) > tStar = Sqrt[am^3/gm] (ecc - em Sin[ecc]); (* exact: 1645.4287 s *) > Do[Print[{n, Abs[transitTrap[0, Pi/2, am, em, n] - tStar], > Abs[transitSimp[0, Pi/2, am, em, n] - tStar]}], {n, {10, 20, 40, 80, 160, 320}}] > tm = 2 Pi Sqrt[am^3/gm]; (* full period, 11.966 h *) > Do[Print[{n, Abs[transitTrap[0, 2 Pi, am, em, n] - tm]}], {n, {16, 32, 64}}] > ``` > *Look for:* errors at $n=320$: $6.2\times10^{-3}$ and $5.5\times10^{-8}$ s; full-orbit errors collapse ($2.2$, $9\times10^{-6}$, $7\times10^{-12}$ s). > > **Step 5: plot the convergence.** > ```wolfram > ns = {10, 20, 40, 80, 160, 320}; > errT = Table[{n, Abs[transitTrap[0, Pi/2, am, em, n] - tStar]}, {n, ns}]; > errS = Table[{n, Abs[transitSimp[0, Pi/2, am, em, n] - tStar]}, {n, ns}]; > ListLogLogPlot[{errT, errS}, Joined -> True, PlotMarkers -> Automatic, > AxesLabel -> {"n", "error (s)"}, PlotLegends -> {"trapezoid", "Simpson"}] > ``` > *Look for:* straight parallel lines on log-log axes. What the guided run should produce: ![convergence of trapezoid and Simpson on the Molniya quarter arc](Media/orbital_convergence.png) > [!warning] Verify against ground truth (required) > Record these checked numbers; "it ran" is not verification. > - Warm-up: integrand $\tfrac49, 1, 4$ at $e=0.5$ (helper at $0$); Earth perihelion/aphelion $147.06/152.14$ million km; comet ratio $19$; ISS $T=5554.7$ s, $T/4=1388.672$ s, $h=51{,}981.8$ km²/s. > - Molniya geometry: $r(0)=6905.6$ km, $r(\pi)=46{,}214.4$ km. > - Circular quarter orbit (ISS): trapezoid $=1388.672$ s $=T/4$ exactly (error $\sim10^{-13}$, machine level). > - Molniya quarter arc: exact $t^\*=1645.4287$ s from Kepler's equation; at $n=320$ the trapezoid error is $6.2\times10^{-3}$ s and the Simpson error $5.5\times10^{-8}$ s. > - Measured orders: error ratios per doubling of $n$ are $4.00$ (trapezoid) and $16.0$ (Simpson), i.e., orders $2.00$ and $4.00$. > - Full Molniya period: $T=2\pi\sqrt{a^3/GM}=11.966$ h, and your full-orbit quadrature must reproduce it. > [!note] A surprise worth seeing: the trapezoid rule on a full orbit > Integrate the full orbit $[0,2\pi]$ with the plain trapezoid rule and compare to $T=11.966$ h: the error is $2.2$ s at $n=16$, $9\times10^{-6}$ s at $n=32$, and $7\times10^{-12}$ s at $n=64$, collapsing far faster than $O(h^2)$. The integrand is smooth and periodic, and for periodic integrands over a full period the trapezoid rule converges spectacularly (the [Euler-Maclaurin](https://en.wikipedia.org/wiki/Euler%E2%80%93Maclaurin_formula) correction terms cancel). The $O(h^2)$ story is for arcs with genuine endpoints, which is why the convergence study above uses a quarter arc. ### Generalization: the Molniya orbit, where Kepler's second law bites (Step 6) [Molniya orbits](https://en.wikipedia.org/wiki/Molniya_orbit) ($a=26{,}560$ km, $e=0.74$, $T\approx12$ h) are used by communications satellites serving high latitudes precisely because the satellite crawls near apogee. Quantify it: compute the transit time for the same $60°$ arc centered at perigee and at apogee. **Scenario A (equal arcs, unequal times).** *Prediction (before running):* the transit-time ratio should exceed the speed ratio. The satellite at apogee is both slower and farther out (the same $d\theta$ subtends a longer path), so the two factors compound. > [!example]- Step 6 code, Python (reference) > > **Step 6: perigee versus apogee dwell.** The same $60°$ arc, timed at both ends of the orbit; speeds from $v=h/r$. > ```python > arc = np.pi/3 > t_fast = transit_simp(-arc/2, arc/2, a_m, e_m, 1000) > t_slow = transit_simp(np.pi - arc/2, np.pi + arc/2, a_m, e_m, 1000) > print(t_fast/60, t_slow/60, t_slow/t_fast) # 12.51 435.40 34.80 > h_m = h_of(a_m, e_m) > print(h_m/r_of(0, a_m, e_m), h_m/r_of(np.pi, a_m, e_m)) # speeds 10.02, 1.50 km/s > same("speed ratio equals (1+e)/(1-e)", > (h_m/r_of(0, a_m, e_m))/(h_m/r_of(np.pi, a_m, e_m)), (1 + e_m)/(1 - e_m)) > ``` > *Look for:* the same $60°$ arc takes $12.51$ minutes at perigee and $435.40$ minutes ($7.26$ h, well over half the $12$ h period) at apogee, a factor of $34.80$; speeds $10.02$ and $1.50$ km/s with ratio exactly $(1+e)/(1-e)=6.69$ (conservation of angular momentum, $v=h/r$). > [!example]- Step 6 code, MATLAB > > **Step 6: perigee versus apogee dwell.** > ```matlab > arc = pi/3; n = 1000; h_m = h_of(a_m, e_m); > lims = [-arc/2, arc/2; pi-arc/2, pi+arc/2]; % perigee arc; apogee arc > t_arc = zeros(2, 1); > for j = 1:2 > th = linspace(lims(j,1), lims(j,2), n+1); dth = (lims(j,2) - lims(j,1))/n; > y = r_of(th, a_m, e_m).^2 / h_m; > t_arc(j) = dth/3*(y(1) + 4*sum(y(2:2:end-1)) + 2*sum(y(3:2:end-2)) + y(end)); > end > fprintf('%.2f %.2f %.2f\n', t_arc(1)/60, t_arc(2)/60, t_arc(2)/t_arc(1)) % 12.51 435.40 34.80 > fprintf('%.2f %.2f\n', h_m/r_of(0, a_m, e_m), h_m/r_of(pi, a_m, e_m)) % 10.02 1.50 km/s > same('speed ratio equals (1+e)/(1-e)', ... > (h_m/r_of(0, a_m, e_m))/(h_m/r_of(pi, a_m, e_m)), (1 + e_m)/(1 - e_m)) > ``` > *Look for:* the same numbers as the Python reference. > [!example]- Step 6 code, R > > **Step 6: perigee versus apogee dwell.** > ```r > arc <- pi/3 > t_fast <- transit_simp(-arc/2, arc/2, a_m, e_m, 1000) > t_slow <- transit_simp(pi - arc/2, pi + arc/2, a_m, e_m, 1000) > cat(t_fast/60, t_slow/60, t_slow/t_fast, "\n") # 12.51 435.40 34.80 > h_m <- h_of(a_m, e_m) > cat(h_m/r_of(0, a_m, e_m), h_m/r_of(pi, a_m, e_m), "\n") # 10.02 1.50 km/s > same("speed ratio equals (1+e)/(1-e)", > (h_m/r_of(0, a_m, e_m))/(h_m/r_of(pi, a_m, e_m)), (1 + e_m)/(1 - e_m)) > ``` > *Look for:* the same numbers as the Python reference. > [!example]- Step 6 code, Mathematica (numerical only) > > **Step 6: perigee versus apogee dwell.** > ```wolfram > arc = Pi/3.; > tFast = transitSimp[-arc/2, arc/2, am, em, 1000]; > tSlow = transitSimp[Pi - arc/2, Pi + arc/2, am, em, 1000]; > Print[{tFast/60, tSlow/60, tSlow/tFast}] (* {12.51, 435.40, 34.80} *) > hm = hOf[am, em]; > Print[{hm/rOf[0., am, em], hm/rOf[N[Pi], am, em]}] (* {10.02, 1.50} km/s *) > same["speed ratio equals (1+e)/(1-e)", > (hm/rOf[0., am, em])/(hm/rOf[N[Pi], am, em]), (1 + em)/(1 - em)] > ``` > *Look for:* the same numbers as the Python reference. **Scenario B (eccentricity drives the cost of accuracy).** Fix the quarter arc $[0,\pi/2]$ and re-run the Step 4 convergence study for $e=0,\ 0.3,\ 0.6,\ 0.9$ (exact values from Kepler's equation each time). Find the $n$ each eccentricity needs for a $10^{-3}$ s error. *Prediction:* higher $e$ concentrates the integrand's variation near perigee, so the constant in the $O(h^2)$/$O(h^4)$ error grows with $e$ and the required $n$ climbs; the order itself should not change. > [!warning] Verify the generalization (required) > - The $60°$ arcs: $12.51$ min at perigee, $435.40$ min at apogee, ratio $34.80$ (Simpson at $n=1000$ matches the Kepler-exact values to $\sim10^{-9}$). > - Speeds $10.02$ and $1.50$ km/s; their ratio must equal $(1+e)/(1-e)=6.6923$ through the helper. > - Scenario B: the measured orders stay $2.00$ and $4.00$ at every $e$; only the constants (and hence the required $n$) grow. ### Your orbits: three eccentricity personalities Now run the study on three real orbits you choose yourself, picked in advance to have three different eccentricity personalities. Look up each orbit's $a$ and $e$ (cite your source), use the right central body ($GM_{\text{Earth}}=398{,}600$ km³/s²; $GM_{\text{Sun}}=1.327\times10^{11}$ km³/s²), and read the difficulty dial $\bigl(\tfrac{1+e}{1-e}\bigr)^2$ *before* running anything: the order of operations matters every time, predict from the dial first, then compute. 1. **Orbit A, near-circular ($e\le0.05$).** A GPS satellite or a LEO orbit. Prediction to test: the integrand is nearly flat, tiny $n$ already reaches tight tolerances, and the perigee/apogee dwell ratio is close to $1$. (A nice pairing: GPS satellites sit at nearly the Molniya $a\approx26{,}560$ km, so they share the $12$ h period; everything that differs from your Molniya runs is eccentricity.) 2. **Orbit B, working eccentric ($0.2\le e\le0.7$).** A geostationary transfer orbit, a Tundra orbit, or a Molniya cousin. Prediction to test: this is the regime where Simpson's extra order visibly pays; state the expected dwell ratio from $e$ before timing the arcs. 3. **Orbit C, extreme ($e\ge0.9$).** A comet ([Halley](https://en.wikipedia.org/wiki/Halley%27s_comet), $e\approx0.967$) or a sungrazer. Prediction to test: the dial reads in the thousands, the required $n$ for a fixed tolerance climbs dramatically, and the dwell ratio becomes astronomical; check whether the *orders* (not the constants) survive. For each orbit, report: the dial value; a quarter-arc transit time against its Kepler-equation exact value (through the helper); the $n$ each rule needs for a $10^{-3}$ s error; and the $60°$ perigee/apogee dwell ratio. Close with a cross-orbit comparison: rank the three by the cost of accurate quadrature and defend the ranking from the convergence numbers and the dial, not from the pictures. #### Student task loop for implementation, analysis and reflection 1. **Predict** the difficulty ranking and dwell ratios from $e$ and the dial, before running. 2. **Implement** the workflow (Steps 2 to 5) on the orbit. 3. **Compare** the transit times against the Kepler-exact values and the required $n$ against your prediction. 4. **Interpret** what the numbers mean for the physical orbit and for the numerical method. ## Reflection Framework Address these (2-3 focused questions per category, no more): ### Convergence and cost - Your convergence data give straight lines on log-log axes. What do the slopes mean, and why does Simpson's line fall so much faster? - When is Simpson's extra bookkeeping worth it over the trapezoid rule, and when is it overkill? ### Geometry and physics - How do your perigee/apogee transit times express Kepler's second law without ever mentioning area? - Why does a Molniya satellite's ground station see it hang in the sky for hours per orbit? ### Mathematical insights - The integrand's max/min ratio $\bigl(\tfrac{1+e}{1-e}\bigr)^2$ acts like a difficulty dial (a conditioning statement for the quadrature problem). How did it show up in Scenario B and across your three orbits? - The full-orbit trapezoid result beat its own error bound by orders of magnitude. What property of the integrand (not the rule) explains that? ## (Optional) Mathematical Extensions These go beyond the project's stated level; attempt them if interested. - **Kepler's equation as root finding:** given $t$, solve $E-e\sin E=t\sqrt{GM/a^3}$ for $E$ with Newton's method, joining this project to the root-finding thread. - **Derive the eccentric-anomaly substitution** that produced the exact $t(\theta)$, and verify it reproduces your ground-truth values. - **Adaptive quadrature:** subdivide only where the integrand varies (near perigee) and compare function-evaluation counts with the uniform grids. - **Periodic superconvergence:** read about the [Euler-Maclaurin formula](https://en.wikipedia.org/wiki/Euler%E2%80%93Maclaurin_formula) and explain the full-orbit trapezoid result quantitatively. ## (Optional) Real-World Context ### Applications - **[International Space Station operations](https://en.wikipedia.org/wiki/International_Space_Station):** pass predictions, crew scheduling, and resupply rendezvous all start from transit-time integrals like yours. - **[Molniya communications satellites](https://en.wikipedia.org/wiki/Molniya_orbit):** the apogee dwell you computed is the entire design rationale for high-latitude coverage. - **[GPS](https://en.wikipedia.org/wiki/Global_Positioning_System):** positioning accuracy depends on continuously updated orbital parameters; the underlying propagation is numerical integration. - **[Space debris tracking](https://en.wikipedia.org/wiki/Space_debris):** collision screening propagates thousands of orbits numerically every day. ### Technical challenges - **[Orbit propagation](https://en.wikipedia.org/wiki/Orbit_modeling):** real orbits feel drag, Earth's oblateness, and third bodies; the exact Kepler formula dies, quadrature and ODE integrators survive. - **[Kepler's equation](https://en.wikipedia.org/wiki/Kepler%27s_equation):** even the "exact" route needs numerics in reverse (position from time is transcendental and is solved by root finding, our other course thread). ### Why it matters A Molniya ground station gets roughly $7$ hours of usable dwell from the $60°$ apogee arc you timed, from an orbit that laps Earth in $12$; that single number, computed by quadrature, is why the constellation works. The same transit integrals, at tighter tolerances, schedule every ISS resupply docking.