Skip to Content
📄 ESSRT — Read the paper
ReferenceFormulas

Formulas and Sources

The model’s published quantities fall into three classes. Most are derived by the engines — the N-body integration of the planets and Earth’s orbit, and the lunisolar precession channel — and have no closed form: the package, the API, the MCP server and the @essrt/model-values registry serve them per epoch, and the Orbital Calculator evaluates them for any year. A second class are identities: frame arithmetic between quantities the engines produce, exact at every epoch. The third class is the frozen era device — the bounded harmonic bases the certified time-domain machinery runs on (Ledger 3 in Mathematical Foundations). This page states which is which, gives the formulas where formulas exist, and points to the page that owns each subject.

This page establishes:

Reproducing a value. Every number on this site carries a registry key; the registry is the published @essrt/model-values package, generated from the engine. The functions named below are the @essrt/physics package’s createModel() surface, the same code the 3D simulation runs; the API and the MCP server expose the same calls.


1. Where each quantity comes from

QuantitySourceObtain itJ2000
Obliquity ε(t)the obliquity hybrid — the spin axis integrated under one averaged precession equation on the engine’s own orbit-plane motion (Obliquity)earth.obliquityDeg(year)23.4393°
Eccentricity e(t)the N-body series — Earth’s eccentricity vector, a sum of the planets’ secular modes anchored at J2000 (Eccentricity)earth.eccentricity(year)0.01670244 (IAU mean element 0.01671022)
Longitude of perihelion ϖ(t)the same series, the argument of the vector (Precession)earth.perihelionLongitudeDeg(year)102.947°
Inclination and node on the invariable planethe N-body chain’s secular skeleton (Moon & Planets)earth.inclinationDeg(year), earth.ascendingNodeDeg(year)—
The planets’ elements of datethe N-body chain — anchors, secular modes and periodic terms in a governed artifactplanets—
Axial precession T_pthe composed lunisolar rate — spin from the recession history, the lunar torque on the Moon’s distance (Earth’s Clock)lunisolar25,771 yr
Apsidal period T_apsthe movement’s own apsidal rate of date — the beat of the one-family anomalistic and sidereal years, T_anom·T_sid/(T_anom − T_sid), in Julian years; a J2000 reading of the wandering apsidal motionregistry inclPrecYears~111,635 yr
Perihelion-of-date periodthe identity 1/T_peri = 1/T_p + 1/T_aps (§2)registry periPrecYears~20,938 yr
Obliquity beatthe hybrid’s dominant period, the beat of T_p against the nodal mode s₃registry obliqCycleYears~41,098 yr
Ecliptic (nodal) periodthe chain’s nodal moderegistry eclPrecYears~68,751 yr
Year lengths of datethe one-source year-length family — the series’ mean motion with the D6 sidereal-year channel (Days & Years)yearLengths.tropicalYearSecondsAtYear(year), siderealYearSecondsAtYear, anomalisticYearSecondsAtYear365.2422036 · 365.2563644 · 365.2596324 d
Day lengthsthe length-of-day stack — the tidal mean with the GIA and cycle terms, plus the identities of §3lengths.lodRealSeconds(year), siderealDaySeconds, stellarDaySeconds86,400.001880 · 86,164.092407 · 86,164.100782 s on LOD_real of J2000 (on the SI day of 86,400 s the IAU constants are 86,164.090531 · 86,164.098904 s)
ΔTthe ΔT stack — the integral of the day-length history plus the cycle termsepoch.deltaTSecondsAtYear(year)55.64 s anchor
Cardinal instants (the seasons)the apparent Sun’s crossings (§7)cardinal.jd(year, 'VE' | 'SS' | 'AE' | 'WS')—
The apparent Sunthe certified solar longitude (§7)eclipse.sunApparentLonDegAtJD(jd)—
The climate formulathe fitted per-regime amplitudes on the engine’s own lines (§6)the coefficient file—

2. The precession identities

Three of the model’s periods are identities between year lengths — the “coin rotation” relations — and one is an identity between periods. They are frame arithmetic and hold at every epoch; their values are the year lengths of date put through them, so at J2000 they are readings, not laws.

Tp=YsidYsid−Ytrop,Tperi=YanomYanom−Ytrop,Taps=YanomYanom−YsidT_p = \frac{Y_{\text{sid}}}{Y_{\text{sid}} - Y_{\text{trop}}}, \qquad T_{\text{peri}} = \frac{Y_{\text{anom}}}{Y_{\text{anom}} - Y_{\text{trop}}}, \qquad T_{\text{aps}} = \frac{Y_{\text{anom}}}{Y_{\text{anom}} - Y_{\text{sid}}}

with the years in the same unit (SI seconds of date). The three periods obey the addition identity

1Tperi=1Tp+1Taps\frac{1}{T_{\text{peri}}} = \frac{1}{T_p} + \frac{1}{T_{\text{aps}}}

— the perihelion-of-date period is the meeting frequency of the precession clock and the inertial apsidal motion, which turn in opposite senses (Precession).

PeriodJ2000What it is
T_p, axial precession25,771 yrthe clock: the composed lunisolar precession period
T_aps, apsidal~111,635 yrthe movement’s own apsidal rate of date (the anomalistic–sidereal beat) — a J2000 reading; the apsidal motion wanders (N-body Exploration)
T_peri, perihelion-of-date~20,938 yrthe beat, 1/T_p + 1/T_aps
T_aps / T_p4.332the ratio the founding observation reads at J2000

3. The day-length identities

Given the tropical year of date Y_trop(s) in SI seconds and the mean solar day D = 86,400 s:

Dsid=Ytrop(s)Ytrop(s)/86400+1,Dstar=Dsid(1+cos⁡εTp⋅R),R=Ytrop(s)DsidD_{\text{sid}} = \frac{Y_{\text{trop}}(s)}{Y_{\text{trop}}(s)/86400 + 1}, \qquad D_{\text{star}} = D_{\text{sid}}\left(1 + \frac{\cos\varepsilon}{T_p \cdot R}\right), \quad R = \frac{Y_{\text{trop}}(s)}{D_{\text{sid}}}

The sidereal day is one rotation against the equinox — shorter than the solar day because Earth turns about one extra degree a day to keep up with its orbit. The stellar day is one rotation against the fixed stars, longer than the sidereal day by the equinox’s westward precession, projected onto the equator by cos ε (precession in right ascension, m = p cos ε). At J2000 the difference is 8.37 ms. On the model’s physical day of J2000 (LOD_real, 86,400.001880 s) the two days are 86,164.092407 s and 86,164.100782 s; on the SI day of exactly 86,400 s the IAU constants are 86,164.090531 s and 86,164.098904 s — the same two days 1.88 ms apart, the day’s excess. The measured mean solar day the day bases close on is 86,400.000427 s (Days & Years).


4. The frozen era device — the harmonic bases

The certified time-domain machinery — the cardinal-point instants, the year-length combs and the ΔT stack’s cycle terms — runs on bounded harmonic bases: Fourier series on fixed fractions of the fitted anchor (Earth’s Clock), with coefficients fitted by the published pipeline against the model’s own N-body movement and the Espenak–Meeus ΔT canon. No polynomial term is admitted, because polynomials compound at deep time; the bases stay bounded at every epoch. This is Ledger 3: four fitted constants and the bases’ coefficients, none of which enters the published Earth orbit or the published clock.

The tropical-year comb, for example, is twelve lines around the mean tropical year; its three largest sit at the device’s obliquity line (41,915 yr, about 2 s of amplitude), its apsidal line (111,772 yr, under 1 s) and its perihelion-of-date line (20,957 yr, about half a second), with the nine remaining lines each under a tenth of a second. The device’s lines sit within 0.1 to 2.5 % of the published periods of the same names (the obliquity beat ~41,098 yr, the apsidal period ~111,635 yr, the perihelion-of-date period ~20,938 yr); the two are different objects — one a fitted basis, the other the engine’s dynamics — and the model keeps them apart by name.

The coefficients live in the public repository’s public/input/fitted-coefficients.json, the single source of truth for fitted values, and the pipeline that writes them is documented in tools/fit/README.md.


5. Planetary precession as seen from Earth

A planet’s perihelion rate read from Earth differs from its ecliptic rate by the reference frame: the equatorial projection of the perihelion motion, with Earth’s obliquity of date the only Earth input. This is a typed identity, not Mercury’s relativistic anomaly, which the model carries as a derived first-order post-Newtonian term (Mercury Precession).

α˙=λ˙ cos⁡εcos⁡2λ+sin⁡2λcos⁡2ε  +  −sin⁡λcos⁡λsin⁡εcos⁡2λ+sin⁡2λcos⁡2ε ε˙\dot\alpha = \dot\lambda\,\frac{\cos\varepsilon}{\cos^2\lambda + \sin^2\lambda\cos^2\varepsilon} \;+\; \frac{-\sin\lambda\cos\lambda\sin\varepsilon}{\cos^2\lambda + \sin^2\lambda\cos^2\varepsilon}\,\dot\varepsilon

Where:

  • λ˙\dot\lambda = the planet’s sidereal perihelion rate. The identity is a frame projection and holds whichever rate is fed into it; what is fed in here is the N-body chain’s own 1800–2100 window with the relativistic term off, so that the relativistic advance stays an explicit, separately stated addition rather than one already folded into the projected number
  • λ\lambda = the planet’s J2000 perihelion longitude advanced at that rate
  • ε\varepsilon, ε˙\dot\varepsilon = the model’s obliquity of date and its rate

The rate fed in is the sidereal one, so the table gives that rather than a period: a perihelion period on this page would otherwise be read against the equinox of date, which is a different quantity and about ten times faster.

PlanetNewtonian sidereal rate λ̇ (″/cy)θ₀ (J2000)
Mercury529.077.455°
Venus14.1131.865°
Mars1,596.5336.099°
Jupiter488.5375.558°
Saturn-1,576.7449.657°
Uranus806.6170.531°
Neptune9,101.7397.442°

Output at J2000 — the Earth-frame right-ascension rate, the ecliptic baseline projected into the simulation’s equatorial frame plus the obliquity-rate term (not the ecliptic rate observers publish):

PlanetEarth-frame RA rate (″/century)Baseline, eclipticEarth-frame excess
Mercury~+576+529.0~+47
Venus~+414.1~-10
Mars~+1,497+1,596.5~-100
Jupiter~+458+488.5~-30
Saturn~-1,718-1,576.7~-142
Uranus~+740+806.6~-66
Neptune~+8,872+9,101.7~-229

For the physical interpretation — why these are reference-frame effects, calculable entirely from Earth’s motion — see Mercury Precession: Perihelion Precession Across the Solar System. The two kinds of planetary “perihelion rate” — the long-term mean and the window trend — are separated on the N-body Exploration page.


6. The climate formula

The climate formula predicts paleoclimate proxy values (benthic δ¹⁸O, atmospheric δ¹³C, ice-core CO₂) at any age tt inside a fitted window — the formula plotted in the Earth Climate Analysis panel of the 3D Simulation :

Cnorm(t)  =  c0  +  ∑i=128[aicos⁡ ⁣2πtTi+bisin⁡ ⁣2πtTi]  +  ∑j∈TL3γj⋅H(t−tj),C(t)  =  Cnorm(t)⋅σy  +  μy  +  (m⋅t+b)C_{\text{norm}}(t) \;=\; c_0 \;+\; \sum_{i=1}^{28} \left[a_i \cos\!\frac{2\pi t}{T_i} + b_i \sin\!\frac{2\pi t}{T_i}\right] \;+\; \sum_{j \in T_{L3}} \gamma_j \cdot \mathcal{H}(t - t_j), \qquad C(t) \;=\; C_{\text{norm}}(t) \cdot \sigma_y \;+\; \mu_y \;+\; \bigl(m \cdot t + b\bigr)

The 28 line periods TiT_i are the engine’s own secular beats — five climatic-precession lines p+gᵢ, ten obliquity lines p+sᵢ, ten planetary eccentricity beats gᵢ−gⱼ, and the 405.6-kyr g₂−g₅ metronome with its harmonics — read from the model’s N-body modes and precession clock and fixed before any fit; the amplitudes and the Heaviside steps are fitted per regime. The formula is valid only inside the regime window it was fitted on. The line list, the per-regime coefficients, the validity windows and the fit statistics live in public/input/climate-formula-coefficients.json and the Python reference implementation scripts/milankovitch_climate_formula.py; the scientific narrative and the per-regime results are on Climate Formula.


7. The apparent Sun

The apparent solar longitude driving the eclipse chain, the cardinal-point instants and the visible simulator is assembled with zero fitted solar constants:

λ(t)=L0+ntrop⋅(t−t0)+D(t)+EoC(e(t), M(t))\lambda(t) = L_0 + n_{trop} \cdot (t - t_0) + D(t) + \mathrm{EoC}\big(e(t),\, M(t)\big)
  • L0L_0 = 280.46646° (J2000 anchor); ntropn_{trop} = 360° per mean tropical year.
  • D(t)D(t) — the year-length drift integral: the closed-form integral of the tropical-year comb (§4), whose two largest lines derive from Acot⁡εA\cot\varepsilon with AA = 0.63607° and stand in the exact ratio of their periods, each scaled by the derived obliquity-torque factor.
  • e(t)e(t), ϖ(t)\varpi(t) — Earth’s elements from the N-body series plus the derived mean-element offset of the era (Δϖ = +70.08″, Δe = +5.81e-6).
  • EoC\mathrm{EoC} — the classical Kepler equation of centre, no fitted coefficients.

Verified at 1.03″ RMS vs JPL in the modern window (the Meeus Ch. 25 reference on the same instrument: 1.02″). Full presentation: The Derived Sun.


8. Verification

CheckModelReference
Obliquity at J200023.4393°IAU 2006: 23.439279°
Eccentricity at J20000.01670244 (secular element)IAU mean element 0.01671022
Longitude of perihelion at J2000102.947°JPL mean element
Obliquity at −10,00024.1589°La2004: 24.1592°, Chapront: 24.3053°
Obliquity over the last 100,000 years1.2″ rmsagainst La2004
Obliquity over the last 1,000,000 years13.5″ rmsagainst La2004
The apparent Sun, modern window1.03″ rmsagainst JPL

The Orbital Calculator computes every quantity of §1 for any year, with the J2000 reference values side by side; the Interactive 3D Simulation shows the same state visually, and its Framework Verification panel plots each quantity against the reference solutions. The research record behind the derivations is doc 109  and the calculation map doc 110  in the public repository.


Return to Mathematical Foundations or explore the Appendix for reference data.

Last updated on: