"""Reckoning 002, follow-up (criticism by Tern, 2026-10-24): is the bisected misfit robust to SMOOTH, correlated errors, not only independent noise?
Two tests, on the same ephemeris stand-in data, all twelve oppositions:
 T1 (Tern's): add A*sin(k*L - phi) to the longitudes (k = 1, 2, 3; random phi) plus 1' white noise; fit the free-equant and bisected models; report how often the
    free model still fits better and by how much.
 T2 (adversarial): for k and amplitude A, CHOOSE phi to make the bisected model fit best; report the smallest A for which the bisected rms falls to <= 1.5'
    (the free model's rms on the clean data is 0.42'). That is the size of smooth systematic needed to rescue the bisected hypothesis.
Run (from reckoning/): python3 kepler_systematic.py        (needs scipy; ~2 minutes)
"""
import math, numpy as np
from kepler_vicarious import *
allidx = list(range(12)); L = opp_l
def rms_for(target, bisected):
    par, nn, cost = fit_free(allidx, target, bisected=bisected)
    r = wrap(geom(par, opp_t, nn)[0] - target) * ARC
    return math.sqrt(np.mean(r ** 2))
rng = np.random.default_rng(5)
print('T1: smooth systematic A*sin(kL - phi), random phi (30 draws each) + 1 arcmin white noise')
print('  k  A(arcmin)  free rms   bisected rms   bisected worse in')
for k in (1, 2, 3):
    for A in (2, 4, 6):
        f, b, worse = [], [], 0
        for _ in range(30):
            ph = rng.uniform(0, 2 * math.pi); tgt = L + (A * np.sin(k * L - ph) + rng.normal(0, 1, 12)) / ARC
            rf, rb = rms_for(tgt, False), rms_for(tgt, True); f.append(rf); b.append(rb); worse += rb > rf
        print('  %d  %5.0f      %6.2f     %6.2f        %2d/30' % (k, A, np.mean(f), np.mean(b), worse))
print('\nT2: smallest smooth systematic (best phase) that lets the BISECTED model fit to <= 1.5 arcmin rms')
for k in (1, 2, 3):
    need = None
    for A in (1, 2, 3, 4, 5, 6, 8, 10, 12):
        best = 1e9
        for ph in np.radians(np.arange(0, 360, 15)):
            best = min(best, rms_for(L + A * np.sin(k * L - ph) / ARC, True))
        if best <= 1.5: need = (A, best); break
    print('  k=%d: %s' % (k, ('A = %d arcmin (bisected rms %.2f)' % need) if need else 'not rescued up to A = 12 arcmin'))