Arc Length of an Ellipse
The arc length along an ellipse with semi-major axis $a$ and parameter $m = e^2$ (where $e$ is the eccentricity) from the end of the major axis to eccentric angle $\varphi$ is
$$s(\varphi) = a\,E(\varphi \mid m) = a\!\int_0^{\varphi}\!\sqrt{1 - m\sin^2\theta}\,\mathrm{d}\theta.$$The complete perimeter is $4a\,E(m) = 4a\,E(\pi/2\mid m)$. Using Carlson forms the evaluation is exact for all $m\in[0,1)$:
% addpath('src') % if not yet on the path a = 1.5; % semi-major axis [AU] ecc = 0.6; % eccentricity m = ecc.^2; % parameter m = e² phi = linspace(0, pi/2, 200); % E(phi|m) — incomplete elliptic integral of the second kind [~, E_phi] = elliptic12(phi, m .* ones(size(phi))); arc_length = a .* E_phi; % s(φ) = a · E(φ|m) % Complete perimeter: 4a·E(m) [~, E_complete] = ellipke(m); perimeter = 4 .* a .* E_complete;
The Solar System — Ellipses as a First Approximation
Each planet travels on an ellipse with the Sun at one focus — at least to first approximation. The orbital elements (semi-major axis $a$, eccentricity $e$, inclination $i$, longitude of ascending node $\Omega$, argument of perihelion $\omega$) are well-known from centuries of observation. The interactive scene below uses the J2000 mean elements from the JPL Planetary Ephemeris (Standish & Williams, 1992) to place each planet on its correct orbit.
But this is only an approximation: the mutual gravitational attraction of the planets slowly perturbs the orbital elements over time. The most famous consequence is the precession of Mercury's perihelion, described in the next section.
% J2000 mean elements for Earth (Standish & Williams 1992) a = 1.00000018; % semi-major axis [AU] ecc = 0.01673163; % eccentricity inc = 0; % inclination [rad] (≈0 for Earth) Om = -5.1126*pi/180; % longitude of ascending node ob = 102.930*pi/180; % longitude of perihelion L0 = 100.467*pi/180; % mean longitude at J2000 n = 35999.373*pi/(180*36525); % mean motion [rad/day] % Position at day d from J2000 d = 0:365; M = mod(L0 + n.*d - ob, 2*pi); % mean anomaly E = M + ecc.*sin(M); % eccentric anomaly (1st-order) % refine with ellipj or Newton-Raphson for better accuracy nu = 2.*atan2(sqrt(1+ecc).*sin(E/2), sqrt(1-ecc).*cos(E/2)); r = a .* (1 - ecc.*cos(E)); % heliocentric distance [AU] x = r .* cos(nu + ob - Om); % ecliptic x [AU] y = r .* sin(nu + ob - Om);
Mercury's Perihelion Precession — Model vs Reality
In a purely Keplerian two-body system, Mercury's orbit would be a fixed ellipse, repeating identically every 87.97 days. In reality, the gravitational pull of the other planets causes the orbital elements to drift slowly — the most visible effect being a steady precession of the perihelion (the point of closest approach to the Sun).
By the late 19th century astronomers had measured a total precession of $\approx 5600$ arcseconds per century. Newtonian mechanics, accounting for all known planetary perturbations, predicted $\approx 5557$ "/cy. The residual $\boldsymbol{43}$" /cy remained unexplained until Einstein's General Relativity in 1915. It was one of GR's first quantitative predictions.
The figure below shows Mercury's orbit under three models, separated by a time $T$ controlled by the slider:
- Keplerian: fixed elements, perfect ellipse, perihelion never moves.
- Newtonian perturbations: 5557 "/cy precession from planetary tugs.
- Observed (= Newtonian + GR): full 5600 "/cy — matches reality.
% Mercury mean elements (J2000) and secular rates [deg/century] a = 0.38710; ecc = 0.20564; ob0 = 77.456; % longitude of perihelion at J2000 [deg] rate_kep = 0; % Keplerian: no precession rate_newt = 531/3600; % Newtonian perturbations [deg/cy] rate_obs = 574/3600; % Observed (Newtonian + GR) [deg/cy] T_cy = 1000; % time elapsed [centuries] % Perihelion longitude under each model ob_kep = ob0 + rate_kep .* T_cy; % stays fixed ob_newt = ob0 + rate_newt .* T_cy; % 1.475° after 100 centuries ob_obs = ob0 + rate_obs .* T_cy; % 1.594° after 100 centuries % GR residual after T_cy centuries [arcsec] delta_arcsec = (rate_obs - rate_newt) .* T_cy .* 3600; fprintf('GR residual: %.1f arcsec (= %.3f deg)\n', delta_arcsec, delta_arcsec/3600);
Keplerian Orbits and Arc Length
In polar coordinates centred on the Sun (one focus), the orbit is $r(\theta) = a(1-e^2)/(1+e\cos\theta)$. The arc length from perihelion to true anomaly $\theta$ is
$$L(\theta) = a\!\int_0^{E_{\mathrm{anom}}(\theta)}\!\sqrt{1 - e^2\sin^2 u}\,\mathrm{d}u = a\,E\!\left(E_{\mathrm{anom}}(\theta)\,\middle|\,e^2\right),$$where $E_{\mathrm{anom}} = 2\arctan\!\bigl(\sqrt{(1-e)/(1+e)}\,\tan(\theta/2)\bigr)$ is the eccentric anomaly.
a = 1.524; % Mars semi-major axis [AU] ecc = 0.0934; % Mars eccentricity m = ecc^2; theta = linspace(0, 2*pi, 361); % Eccentric anomaly from true anomaly E_anom = 2 .* atan2(sqrt(1-ecc) .* sin(theta/2), ... sqrt(1+ecc) .* cos(theta/2)); E_anom(E_anom < 0) = E_anom(E_anom < 0) + 2*pi; % Arc length L(θ) = a · E(E_anom | m) [~, E_phi] = elliptic12(E_anom, m .* ones(size(E_anom))); L = a .* E_phi; % Full orbital perimeter = 4a·E(m) [~, Ec] = ellipke(m); perimeter = 4 * a * Ec; fprintf('Mars perimeter: %.4f AU\n', perimeter);
Kepler's Equation
Time is encoded in Kepler's equation $M = E - e\sin E$, where $M = n(t-t_0)$ is the mean anomaly (linear in time) and $E$ is the eccentric anomaly. There is no closed-form inverse; we solve numerically via Newton–Raphson:
$$E_{k+1} = E_k + \frac{M - E_k + e\sin E_k}{1 - e\cos E_k}.$$The combined pipeline $t \!\to\! M \!\to\! E \!\to\! \theta \!\to\! L$ gives arc length as a function of time, with $E(\,\cdot\mid e^2)$ evaluated at each step.
a = 1.524; ecc = 0.0934; m = ecc^2; T = 686.97; % orbital period [days] t = linspace(0, T, 500); % Step 1: mean anomaly M (linear in time) M = 2*pi .* t ./ T; % Step 2: solve Kepler's equation M = E − e sin E (Newton–Raphson) E = M; for k = 1:50 E = E + (M - E + ecc.*sin(E)) ./ (1 - ecc.*cos(E)); end % Step 3: arc length L(t) = a · E(E_anom | m) using elliptic12 [~, E_phi] = elliptic12(E, m .* ones(size(E))); L = a .* E_phi; % Step 4: true anomaly for position on orbit nu = 2 .* atan2(sqrt(1+ecc).*sin(E/2), sqrt(1-ecc).*cos(E/2)); r = a .* (1 - ecc.*cos(E)); plot(t, L); xlabel('Days from perihelion'); ylabel('Arc length [AU]');
Associate Integrals $B$, $D$, $J$ — Cancellation-Free Arithmetic
The standard integrals $F$, $E$, $\Pi$ are linear combinations of three associates:
$$B(\varphi|m)=\!\int_0^\varphi\!\frac{\cos^2\!\theta}{\Delta}\mathrm{d}\theta,\quad D(\varphi|m)=\!\int_0^\varphi\!\frac{\sin^2\!\theta}{\Delta}\mathrm{d}\theta,\quad \Delta=\sqrt{1-m\sin^2\!\theta},$$ $$F=B+D,\qquad E=B+(1-m)D,\qquad \Pi=B+D+nJ.$$Computing $E$ as $F - mD$ loses significant digits when $m\approx 1$ because $F\approx mD$ (large cancellation). The form $B+(1-m)D$ keeps all digits since both terms are individually smooth and positive.
phi = linspace(0, pi/2, 200); m = 0.97; % near-singular case where cancellation matters most n = 0.5; % characteristic for J % Incomplete associate integrals B(φ|m), D(φ|m), J(φ,n|m) [B, D, J] = ellipticBDJ(phi, m, n); % Recover F, E, Π without cancellation F = B + D; E = B + (1 - m) .* D; % ← stable; avoids F - m·D subtraction Pi = B + D + n .* J; % Complete versions at φ = π/2 [Bc, Dc, ~, Jc] = ellipticBD(m, n); [K, Ek] = ellipke(m); fprintf('K = B+D = %.15g (ellipke: %.15g)\n', Bc+Dc, K); fprintf('E = B+mcD = %.15g (ellipke: %.15g)\n', Bc+(1-m)*Dc, Ek);
Carlson Duplication: Quadratic Convergence
The $R_F(x,y,z)$ algorithm replaces $(x,y,z)$ by their Landen-transformed mean at each step. After $n$ steps the relative error is $O(4^{-n})$ — each iteration roughly doubles the number of correct decimal digits. Convergence is declared when the scaled deviations $|X|, |Y|, |Z| < 3.2\times10^{-3}$, at which point a 7th-order polynomial in $E_2,E_3$ (elementary symmetric functions of $X,Y,Z$) yields the full double-precision answer.
% Direct Carlson symmetric forms x = 0.5; y = 1.0; z = 2.0; RF = carlsonRF(x, y, z); RD = carlsonRD(x, y, z); % Recover F(φ|m) and E(φ|m) — DLMF 19.25.5 / 19.25.7 phi = pi/3; m = 0.7; s = sin(phi); c = cos(phi); d2 = 1 - m*s^2; F_carl = s .* carlsonRF(c^2, d2, 1); E_carl = s .* carlsonRF(c^2, d2, 1) - (m/3).*s^3.* carlsonRD(c^2, d2, 1); % Cross-check against elliptic12 [F_ref, E_ref] = elliptic12(phi, m); fprintf('F: carlson=%.15g ref=%.15g diff=%.2e\n', F_carl, F_ref, abs(F_carl-F_ref)); fprintf('E: carlson=%.15g ref=%.15g diff=%.2e\n', E_carl, E_ref, abs(E_carl-E_ref)); % R_D encodes complete D(m) = (1/3)·R_D(0, 1-m, 1) m_v = [0.3, 0.7, 0.9]; D_complete = (1/3) .* carlsonRD(0, 1-m_v, 1); [~, D_ref] = ellipticBD(m_v);
Browser Validation Against Octave
The table below compares the JavaScript Carlson implementation running in this
browser against reference values computed by Octave using the elliptic12
function from this library. Differences should be at or below machine epsilon.
phi_tab = [pi/6, pi/4, pi/3, pi/2, pi/6, pi/4, pi/3, pi/2]; m_tab = [0.3, 0.3, 0.3, 0.3, 0.7, 0.7, 0.7, 0.7]; % Legendre incomplete integrals F and E [F, E] = elliptic12(phi_tab, m_tab); % Carlson forms give identical results s = sin(phi_tab); c = cos(phi_tab); d2 = 1 - m_tab.*s.^2; F_carl = s .* carlsonRF(c.^2, d2, ones(size(d2))); E_carl = F_carl - (m_tab/3).*s.^3.* carlsonRD(c.^2, d2, ones(size(d2))); fprintf('max|F diff| = %.2e\n', max(abs(F - F_carl))); fprintf('max|E diff| = %.2e\n', max(abs(E - E_carl))); % Also check associate integrals [B, D] = ellipticBDJ(phi_tab, m_tab); assert(max(abs(B + D - F)) < 1e-13, 'F = B+D failed'); assert(max(abs(B + (1-m_tab).*D - E)) < 1e-13, 'E = B+(1-m)D failed');
| φ | m | F (Octave ref) | F (JS) | |ΔF| | E (Octave ref) | E (JS) | |ΔE| |
|---|
Green rows: $|\Delta| < 10^{-12}$. Red rows would indicate a bug (should not appear).