A unified system for predicting planetary precession fluctuations for all 7 planets (Mercury through Neptune) based on the Holistic Universe Model. For the theoretical derivation of these formulas (Fibonacci hierarchy, resonance loops, term categories), see Formula Derivation.
This system predicts how planetary precession rates deviate from their baseline values over time. The model uses a physical-beat feature matrix (~2421 terms per planet) whose feature frequencies are derived entirely from public/input/model-parameters.json — no hardcoded H_DIV_X constants. It captures the interactions between:
- Earth's perihelion longitude and its rate of change (ERD)
- Planet perihelion longitudes relative to Earth
- Earth's obliquity and eccentricity variations
- The six physical periods per planet (ecliptic perihelion, ICRF perihelion, obliquity cycle, ascending node, axial precession, wobble) plus all pairwise internal and Earth-cross beat frequencies — derived from the Earth Fundamental Cycle (H) via
planet_beats.py - Sideband structure: high-harmonic carriers (fundamentals and ICRF internal-beats) amplitude-modulated by Earth's ecliptic perihelion rate (see GROUP J, K, L below)
All 7 planets use the same feature structure with planet-specific coefficients trained via ridge regression. When JSON periods change, feature frequencies update automatically, but coefficients must be retrained.
The fundamental cycle of H years from which all other periods are derived (see Constants Reference for current value):
| Period | Description |
|---|---|
| H | Master cycle |
| H/3 | Inclination cycle |
| H/8 | Obliquity cycle |
| H/13 | Axial precession |
| H/16 | Earth perihelion |
Each planet's precession period is derived from H using Fibonacci-related fractions:
| Planet | Period Formula | Baseline Rate |
|---|---|---|
| Mercury | H × 8/11 | +531.4 "/cy |
| Venus | H × 8/6 | −289.9 "/cy (retrograde) |
| Mars | H × 8/35 | +1691.0 "/cy |
| Jupiter | H / 5 | +1932.5 "/cy |
| Saturn | H / 8 | −3092.0 "/cy (retrograde) |
| Uranus | H / 3 | +1159.5 "/cy |
| Neptune | H × 2 | +193.3 "/cy |
All planets are calibrated to their official J2000 (year 2000) perihelion longitudes. The perihelion at any year is calculated as:
θ(year) = θ₀ + 360° × (year - 2000) / period
tools/lib/python/ # Python physics library
├── constants_scripts.py # Shared constants (reads from constants.js via load_constants)
├── planet_beats.py # Fundamental periods + beats per planet (derived from JSON)
├── predictive_formula_physical.py # Primary: ~2421-term physical-beat feature builder
├── predictive_formula.py # Legacy 429-term builder + shared helpers (calc_erd, etc.)
├── observed_formula.py # Observed formula: feature builders (225/328 terms)
├── predict_precession.py # Prediction API (Earth + Planets output)
├── validate_precession.py # Validation against observed data
├── coefficients/ # ML coefficient output
│ ├── *_coeffs_physical.py # Physical-beat coefficients (7 files, ~2421 terms) — active
│ └── *_coeffs.py # Observed coefficients (7 files, 225/328 terms)
└── PREDICTIVE_FORMULA_GUIDE.mdx # This guide
tools/fit/ # Fitting scripts
├── python/
│ ├── train_precession_physical.py # Physical-beat predictive training (ridge)
│ ├── train_observed.py # Observed training (SVD least-squares)
│ ├── fit_perihelion_harmonics.py # Earth perihelion harmonics
│ ├── load_constants.py # Python bridge to constants.js
│ └── greedy_features_physical.py # Residual-analysis / candidate features
└── README.md # Dependency chain & refit guide
scripts/archive/ # Legacy 429-term system (superseded)
├── train_precession.py
├── greedy_features.py
└── *_coeffs_unified.py
data/
└── 01-holistic-year-objects-data.xlsx # Training data (Excel, multiple sheets)
From the command line:
cd tools/lib/python
python predict_precession.py 2000Output shows predictions for all planets at the specified year.
from predict_precession import predict, predict_total, get_all_predictions
# Get fluctuation for a single planet
fluctuation = predict(2000, 'mars') # Returns -113.93 "/cy
# Get total precession rate (baseline + fluctuation)
total = predict_total(2000, 'mars') # Returns +1568.07 "/cy
# Get predictions for all planets
results = get_all_predictions(2000)
for planet, data in results.items():
print(f"{planet}: {data['total']:.2f} \"/cy")from predict_precession import get_earth_parameters
earth = get_earth_parameters(2000)
print(f"Perihelion: {earth['perihelion']:.6f}°")
print(f"Obliquity: {earth['obliquity']:.6f}°")
print(f"Eccentricity: {earth['eccentricity']:.10f}")
print(f"ERD: {earth['erd']:.10f} °/yr")from predictive_formula import calc_planet_obliquity
# Dynamic obliquity anchored to J2000 axial tilt
for planet in ['Mercury', 'Venus', 'Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune']:
obliq = calc_planet_obliquity(planet, 2000)
print(f"{planet}: {obliq:.4f}°")Obliquity oscillates with Fibonacci-predicted periods (see doc 37). Venus and Neptune return static values (no obliquity cycle).
from predictive_formula import calc_obliquity
# Primary obliquity — 12-harmonic formula, RMSE ~0.005 arcsec
obliq = calc_obliquity(2000) # 23.439640° at J2000The mean (SOLSTICE_OBLIQUITY_MEAN_FITTED = 23.4534°) is a data-derived constant stored in fitted-coefficients.json. It represents the solstice-instant mean obliquity, which differs from the Pythagorean time-average by ~1.5" due to the instantaneous measurement geometry at solstice. The 12 harmonic coefficients are fitted from simulation data. See doc 14 §4 for the full derivation.
For H/3 vs H/8 component decomposition (dashboard charts), use computeObliquityIntegrals() in the orbital engine.
All four cardinal points (SS, WS, VE, AE) are supported via an optional cp_type parameter (default 'SS' for backward compatibility):
from predictive_formula import calc_solstice_ra, calc_solstice_jd, calc_solstice_year_length
# Where: RA of any cardinal point (fully derived, zero fitted constants)
ra_ss = calc_solstice_ra(2000, 'SS') # 90.0169° (summer solstice ≈ 6h)
ra_ws = calc_solstice_ra(2000, 'WS') # 270.0169° (winter solstice ≈ 18h)
ra_ve = calc_solstice_ra(2000, 'VE') # 0.0169° (vernal equinox ≈ 0h)
ra_ae = calc_solstice_ra(2000, 'AE') # 180.0169° (autumnal equinox ≈ 12h)
# When: Julian Day of any cardinal point (anchored at observed J2000 value)
jd_ss = calc_solstice_jd(2000, 'SS') # 2451716.575000 (Jun 21, 01:48 UTC)
jd_ws = calc_solstice_jd(2000, 'WS') # 2451900.067346 (Dec 21, 01:37 UTC)
# How long: time between consecutive events of the same type
yr_ss = calc_solstice_year_length(2000, 'SS') # 365.24160 days (shortest, −51s)
yr_ws = calc_solstice_year_length(2000, 'WS') # 365.24274 days (longest, +47s)The RA formula uses only H/3 and H/8 harmonics with amplitude A/sin(ε) — a pure geometric projection of the scene graph. The JD formula uses 12 harmonics per cardinal point (5 Fibonacci fundamentals + 7 overtones that are sums of Fibonacci: H/6=3+3, H/11=3+8, H/32=16+16, etc.). The H/16 phase rotates 90° between consecutive cardinal points, encoding the equation of center sweeping through the orbit.
RMSE: RA 0.089° for all types; JD: SS 2.7min, VE 3.0min, AE 5.0min, WS 5.3min over full H. See doc 14 for full derivation.
If you have updated observed data, you can retrain the coefficients:
Before training, verify the Earth perihelion formula produces correct values at J2000:
cd tools/lib/python
python -c "from predictive_formula import calc_earth_perihelion; print(f'{calc_earth_perihelion(2000):.6f}°')"Expected: ~102.94° (within ~0.01° of the Excel Earth Perihelion ICRF value at year 2000). The training script also prints this value automatically.
The training data is in data/01-holistic-year-objects-data.xlsx (Excel format). The step size should divide H (335,317) evenly — currently 23-year steps (14,579 points). Data comes from:
Sheet 'Holistic_objects_PerihelionPlan':
Model Year— Calendar year (must include year 2000){Planet} Precession Fluctuation— Observed fluctuations (″/century) for Mercury, Venus, Mars, Jupiter, Saturn, Uranus, Neptune
Earth orbital parameters (obliquity, eccentricity, perihelion) are computed analytically by tools/lib/python/predictive_formula.py, not read from Excel.
python tools/fit/python/train_precession_physical.py --writeThis will:
- Load observed data from Excel (downsampled by
stepYears=23for speed) - Build the ~2421-term physical-beat feature matrix for each data point
- Fit ridge regression (α=0.01) for each planet
- Write new coefficient files to
tools/lib/python/coefficients/*_coeffs_physical.py - Update
public/input/fitted-coefficients.json(PREDICT_COEFFS_PHYSICALkey) - Display RMSE and R² metrics
Per-planet mode (~7× faster for iteration):
python tools/fit/python/train_precession_physical.py --planet venus --write- Tweakpane (
src/script.js):node tools/fit/export-to-script.js --write— replaces thePREDICT_COEFFSblock in-place. - Holistic website:
node tools/fit/export-to-holistic.js --write— writescoefficients.ts(and verifiesplanets.ts::buildFeaturesterm count matches). - Node verification (
tools/lib/orbital-engine.js): picks up coefficients automatically viatools/lib/constants/fitted-coefficients.js, which reads the JSON.
cd tools/lib/python
python validate_precession.pyThis compares predictions against observed data and shows:
- RMSE, Mean Error, Max Error for each planet
- R² scores (should be > 0.998 for good fit)
- Sample comparisons at key years
The ~2421-term physical-beat feature matrix is organized into 8 groups. All frequencies derive from planet_beats.py::fundamental_periods(planet) — which reads model-parameters.json at import time.
| Group | Terms | Description |
|---|---|---|
| A | 49 | Earth reference: angle harmonics (δ, Σ), obliquity/eccentricity terms, ERD (linear/quadratic/cubic + cross-terms), 3δ/4δ, sum-angle×ERD, constant |
| B | 10 | Planet angle cross-terms: n=3,4 θ_P; obliq/ecc × 2δ; cos(θ_E)×cos/sin(θ_P) |
| C+D+E | 2·7·N_template | Physical periods & beats × 7 harmonics (sin/cos). Template per planet = 6 fundamentals + 15 internal-beat pairs (sum+diff) + 36 Earth-cross-beat pairs (sum+diff). HARMONIC_ORDERS = (1, 2, 4, 6, 8, 12, 16). Venus: 107 template entries → 1498 terms. |
| F | 4·N_fund | Fundamental period × angle: cos(φ)·cos(δ), sin(φ)·sin(δ), cos(φ)·cos(2δ), sin(φ)·sin(2δ) per fundamental |
| I | 4·N_fund | Fundamental × 2δ: sin(φ)·cos(2δ), sin(φ)·sin(2δ), cos(φ)·cos(2δ), cos(φ)·sin(2δ) per fundamental |
| J | 64·N_fund | Fundamental × Earth-cycle sidebands (n=1 carrier). 4 Earth cycles {ecl, obliq, icrf, axial} × 4 harmonics {1,2,3,4} × 4 combinations {sin·sin, sin·cos, cos·sin, cos·cos} per fundamental |
| L | 48·N_fund | High-harmonic fundamental × Earth ecl sidebands. Carrier n ∈ {6,10,12,16} × sideband k ∈ {1,2,3} × 4 combos. Added after residual analysis revealed Venus H/165 (= 12·V_icrf) sideband structure. |
| K | 144 | ICRF internal-beat carriers × Earth ecl sidebands. Carrier = n·(T_ecl ± T_icrf) for n ∈ {2,4,6,8,10,12} × sideband k ∈ {1,2,3} × 4 combos × 2 signs. Added after diagnosing Venus H/116 peak as 8·(V_ecl+V_icrf). |
Removed groups (were present in earlier versions; dropped after 0.00% impact diagnostic):
G (ERD × period, 4·N_fund)andH (ERD × period × angle, 4·N_fund)— 48 features total, zero contribution to any planet's RMSE.
The GROUP K and L feature classes were identified by running greedy_features_physical.py on Venus residuals and decomposing peaks in 8H/N form:
- Venus residual peak at H/116 = 8H/928 decomposed as
8 × (T_ecl + T_icrf)— an 8th-harmonic Venus ICRF internal beat. Sidebands at H/84, H/100, H/132, H/148, H/164 followed the pattern8H/(928 ± k·128)where 128 = 8H/128 = Earth ecliptic perihelion rate. This became GROUP K. - After GROUP K, residual peaks at H/165, H/181, H/197 decomposed as
12·T_icrf ± k·E_ecl— high-harmonic fundamental × Earth-ecl sidebands. This became GROUP L.
This methodology (express residual peaks as 8H/N integers, find the physical decomposition, add matching feature classes) is the recommended way to extend the matrix.
Saturn's precession period (H/8) equals the obliquity cycle, creating strong coupling with Earth's obliquity/eccentricity variations. The Group A obliq/ecc × diff-angle terms and the C+D+E beat combinations involving H/8 are critical.
Current model performance (~2421 terms per planet, 14,580 data points, 23-year steps, ridge α=0.01):
| Planet | RMSE ("/cy) | R² | Feature count |
|---|---|---|---|
| Mercury | 0.0974 | 0.999999 | 2421 |
| Venus | 0.7278 | 0.999998 | 2421 |
| Mars | 0.0980 | 0.999999 | 2435 |
| Jupiter | 0.0975 | 0.999999 | 2407 |
| Saturn | 0.0977 | 1.000000 | 2407 |
| Uranus | 0.1027 | 0.999998 | 2407 |
| Neptune | 0.0978 | 0.999985 | 2393 |
Feature count differs per planet because frozen periods (Uranus/Neptune ~200 Myr axial ≈ treated as infinite → merged with ICRF) and null beats (e.g. ecl+wobble when they coincide) drop template entries.
To add features to the model:
- Run
python tools/fit/python/greedy_features_physical.py --planet venus(or another planet) to find candidate terms with highest residual correlation. Convert residual peaks to 8H/N form and decompose physically (see "How the residual structure was discovered" above). - Edit
tools/lib/python/predictive_formula_physical.py: add a new GROUP block tobuild_features_physical()and update the matching count infeature_count(). Keep Python and JS in sync — the ordering is what makes coefficients portable. - Retrain all planets with
python tools/fit/python/train_precession_physical.py --write. - Port the new group to the JS side:
src/script.js::buildPredictiveFeatures(tweakpane),tools/lib/orbital-engine.js::buildPredictiveFeatures(Node verification), and the Holistic website'ssrc/lib/orbital/planets.ts::buildFeatures. Thetools/fit/export-to-*.jsscripts sync coefficients but not the feature builder structure.
-
Add planet to
model-parameters.json(withperihelionEclipticFraction,axialPrecessionFraction,ascendingNodeCyclesIn8H, andobliquityCycleFractionif applicable) andastro-reference.json(withlongitudePerihelion). -
Add column to
data/01-holistic-year-objects-data.xlsx:NewPlanet Precession Fluctuation. -
Register the planet in
train_precession_physical.py::PLANETS_FLUCT_COLSandpredictive_formula_physical.py::_TEMPLATES. -
Retrain:
python tools/fit/python/train_precession_physical.py --planet newplanet --write(generatestools/lib/python/coefficients/newplanet_coeffs_physical.pyand updates JSON). -
Import in
tools/lib/python/predict_precession.py.
Total precession rate = Baseline + Fluctuation
- Baseline: Mean precession rate derived from planet's period (360°/period per year, converted to "/cy)
- Fluctuation: Time-varying deviation predicted by the model
Key calculated values for each prediction:
- Earth Perihelion (θ_E): Sum of mean motion + Fourier harmonics (
PERI_HARMONICS_RAW) - ERD: Analytical derivative of Earth perihelion (rate deviation)
- Planet Perihelion (θ_P): Linear from J2000 value, period =
|perihelionEclipticYears| - Difference Angle (δ): θ_E - θ_P (key interaction term)
- Obliquity / Eccentricity: Time-varying Earth parameters
- Six physical periods per planet (ecl, icrf, obliq, asc, axial, wobble) + all pairwise internal beats (15 pairs × sum/diff) + all pairwise Earth-cross beats (36 pairs × sum/diff), each at 7 harmonic orders (1, 2, 4, 6, 8, 12, 16)
- Sidebands: high-harmonic carriers × Earth ecliptic perihelion (Group J: n=1 carrier × 4 Earth cycles; Group L: high-n fundamental × Earth ecl; Group K: high-n ICRF internal beat × Earth ecl) — these capture the amplitude-modulation structure that dominates Venus residuals.
The model captures how these parameters interact across multiple timescales to produce the observed fluctuation patterns.
- Holistic Universe Model documentation
- J2000 perihelion values from NASA JPL Horizons
- Training data spans one full Earth Fundamental Cycle (H years) of calculated values (14,579 data points, 23-year steps)