reference_window.py

"""Estimate a sampled agreement window between numerical references.

Varies the precision of one reference, compares an independent
implementation, and masks unsupported pointwise comparisons using existing
reference utilities. See Lesson 2G.

Usage
-----
python -m examples.reference_window

Model
-----
Non-hierarchical Pythagorean three-body state with spatial perturbation 0.01;
401 samples on [0, 20]. IAS15 with epsilon 1e-9 and 1e-12; external REBOUND
BS with rtol 1e-12 and atol 1e-14 (distinct from the core BS). Agreement norm
uses position scale 5, velocity scale sqrt(12/5) and threshold 1e-8.

Notes
-----
Prints crossing brackets (or a censored no-crossing result), finite pointwise
sample counts and energy drift; writes no files. Windows are sampled and
threshold-dependent. They are not rigorous error bounds or Lyapunov times.
Each pair is checked separately. Requires REBOUND.
"""
import numpy as np
from astrosimbench.systems import pythagorean_three_body
from astrosimbench.references import rebound_reference,validity_horizon,pointwise_error
from astrosimbench.diagnostics import invariant_history


def main():
    r,v,m,_=pythagorean_three_body(spatial_perturbation=.01)
    times=np.linspace(0,20,401)
    a=rebound_reference(r,v,m,times,epsilon=1e-9)
    b=rebound_reference(r,v,m,times,epsilon=1e-12)
    c=rebound_reference(r,v,m,times,method='bs',rtol=1e-12,atol=1e-14)
    for label,check in [('IAS15 tolerance change',b),('IAS15 vs BS',c)]:
        w=validity_horizon(a,check,position_scale=5.,velocity_scale=np.sqrt(12/5),threshold=1e-8)
        print(label,':',f'no crossing through t={times[-1]} (censored)' if w.censored else f'crossing bracket [{w.last_valid_sample}, {w.t_valid}]')
        print('pointwise valid samples:',np.count_nonzero(np.isfinite(pointwise_error(a,check,w))))
    for out in (a,b,c):
        history=invariant_history(out,m)
        print(out.method, 'max absolute energy drift:',np.max(np.abs(history['energy_drift'])))
    print('After the earliest observed crossing: use invariants and converged ensemble/distribution diagnostics, not trajectory ranking.')

if __name__=='__main__': main()