INFLATION: RECOMPUTE A PRIMORDIAL SPECTRUM (v1)
Expected files were made with Python 3.9.6, numpy 2.0.2 and scipy 1.13.1.

1. Install requirements.txt in your own virtual environment.
2. Check the solver against closed-form answers first:
     python verify_benchmarks.py my-benchmarks.json
   Power-law inflation has an exact spectrum (Stewart & Lyth 1993, eqs. 31,
   41, 43). The script also checks the canonical Wronskian and moves the start
   ratio k/(aH) and the tolerance. It exits nonzero when a check fails.
3. Recompute one fixed model and compare it with the saved answer:
     python recompute_inflation_spectrum.py inputs/step-ripples.json my-step.json
     python compare_with_expected.py my-step.json expected/step-ripples.json
   Each preset takes about one second on a laptop.

UNITS AND CONVENTIONS
Reduced Planck mass M_pl = 1. Time is N = ln a (e-folds). The pivot
k* = 0.05 Mpc^-1 equals aH when Nstar e-folds remain before epsilon_H = 1.
Each mode starts in the Bunch-Davies state at k/(aH) = 100 and is read out at
k/(aH) = 1e-3. P_R = k^3 |R|^2 / (2 pi^2); P_t = 8 k^3 |h|^2 / (2 pi^2).
The amplitude field is V0 in M_pl^4. P_R scales exactly with V0 at fixed Nstar,
so short-reheating (Nstar = 46) gives A_s = 1.488e-9 with the same V0; rescale
V0 if you want A_s = exp(3.044)/1e10.
The step follows Adams, Cresswell & Easther (2001) eq. 1,
V = V_smooth [1 + c tanh((phi - phi_s)/d)], with phi_s and d given in M_pl.

WHAT THE FILES ARE
inputs/    five fixed models used by the Sandbox Physics Lab presets.
expected/  our own outputs on a 27-point grid from 2e-4 to 0.3 Mpc^-1.
These are model calculations, not observations. Compare them with published
constraints such as Planck 2018 X (arXiv:1807.06211v2), n_s = 0.9649 +/- 0.0042.

A FAILURE CASE WORTH REPEATING
An earlier version used a fixed absolute tolerance atol = 1e-30. For
Nstar = 46 the starting |R| ranges from about 2e-29 to 1e-24 across the grid
(2e-23 to 1e-18 for Nstar = 55), so that tolerance was not small. Over the 27-point grid the scalar power came out
up to 0.55% wrong and the Wronskian drifted by up to 0.88%, while every run
reported success. The Nstar = 55 plateau stayed below 1e-4. The current
code scales atol with each component's starting size. To see the failure,
replace the atol line in solve_mode with a constant 1e-30 and rerun
inputs/short-reheating.json; then check the Wronskian with wronskian_drift().

LIMITS
Single canonical field, Einstein gravity, Bunch-Davies vacuum, no reheating
dynamics after epsilon_H = 1. The spectra are primordial; no CMB transfer
function or likelihood is applied, so no fit to data is claimed.
Course: /stories/inflation-perturbations/recompute-a-primordial-spectrum/
