"""Reckoning 002 (STAND-IN DATA!): Kepler's vicarious hypothesis rebuilt on a MODERN EPHEMERIS (JPL DE441 via Horizons), NOT on Brahe's observations.

Model (longitudes as seen from the Sun S): Mars moves on a circle of radius 1 about centre C; S is displaced e_s from C toward perihelion; the equant E
is e_q from C on the far side (Ptolemy's 'bisected' form has e_s = e_q). Mars's equant angle advances uniformly, w = w0 + n (t - t0). Mars is the point
P of the circle on the ray from E at angle w; its heliocentric longitude is  lon = phi + atan2(P_y, P_x + e_s)  and its distance from S is |P - S|.

What this script does
  A. fit all 12 oppositions (phi, e_s, e_q, w0, n free)  -> how well can a circle+equant represent the OPPOSITIONS?
  B. fit only FOUR oppositions (n fixed from A, 'the mean motion known from the long opposition table') -> residuals at the other eight
  C. the bisected variant e_s = e_q = e (free), fit to all 12 -> Kepler's '8 arcminutes'?
  D. DISTANCES: compare the Sun-Mars distance each model predicts with the ephemeris at the oppositions (Kepler's real reason to reject the unbisected fit)
  E. noise: repeat A and C with 1.5' Gaussian noise on the longitudes (Brahe's accuracy is about 1-2') and count how often C's misfit stands out

Run (from reckoning/):  python3 kepler_vicarious.py     (needs data/kepler_ephem.csv from kepler_fetch.py, and scipy)
"""
import math, itertools, numpy as np
from scipy.optimize import least_squares
from kepler_oppositions import opp, rows, jd

T0 = opp[0][0]
opp_t = np.array([o[0] for o in opp]); opp_l = np.radians([o[1] for o in opp]); opp_r = np.array([o[2] for o in opp])
ARC = 180 * 60 / math.pi
n0 = 2 * math.pi / 686.98   # first guess only; replaced by the fit in A


def wrap(x):
    return (x + math.pi) % (2 * math.pi) - math.pi


def geom(par, t, nn):
    """returns (heliocentric longitude, distance from S) for parameters (phi, es, eq, w0)."""
    phi, es, eq, w0 = par
    w = w0 + nn * (t - T0)
    ux, uy = np.cos(w), np.sin(w)
    Eu = eq * ux
    tt = -Eu + np.sqrt(np.maximum(Eu * Eu - eq * eq + 1.0, 0))
    Px, Py = eq + tt * ux, tt * uy
    return phi + np.arctan2(Py, Px + es), np.hypot(Px + es, Py)


def fit_free(idx, target, bisected=False, fixed_n=None):
    best = None
    for phi0 in np.radians(np.arange(0, 360, 30)):
        for w00 in np.radians(np.arange(0, 360, 90)):
            if bisected:
                def unpack(p):
                    return [p[0], p[1], p[1], p[2]], (p[3] if fixed_n is None else fixed_n)
                x0, lo, hi = [phi0, .09, w00, n0], [-10, -.3, -10, n0 * .99], [10, .3, 10, n0 * 1.01]
                if fixed_n is not None: x0, lo, hi = x0[:3], lo[:3], hi[:3]
            else:
                def unpack(p):
                    return list(p[:4]), (p[4] if fixed_n is None else fixed_n)
                x0, lo, hi = [phi0, .09, .09, w00, n0], [-10, -.3, -.3, -10, n0 * .99], [10, .3, .3, 10, n0 * 1.01]
                if fixed_n is not None: x0, lo, hi = x0[:4], lo[:4], hi[:4]
            xs = [1, .1, .1, 1, 1e-4] if not bisected else [1, .1, 1, 1e-4]
            xs = xs[:len(x0)]
            r = least_squares(lambda p: wrap(geom(unpack(p)[0], opp_t[idx], unpack(p)[1])[0] - target), x0, bounds=(lo, hi), x_scale=xs)
            if best is None or r.cost < best[0].cost:
                best = (r, unpack)
    r, unpack = best
    par, nn = unpack(r.x)
    return par, nn, r.cost


def lon_res(par, nn, idx=None):
    t = opp_t if idx is None else opp_t[idx]
    l = opp_l if idx is None else opp_l[idx]
    return wrap(geom(par, t, nn)[0] - l) * ARC


def dist_res(par, nn):
    return (geom(par, opp_t, nn)[1] / opp_r) - 1.0  # relative error, model unit radius=1 scaled below


if __name__ == '__main__':
    allidx = list(range(12))
    # ---- A
    pA, nA, cA = fit_free(allidx, opp_l)
    rA = lon_res(pA, nA)
    print('A. all 12 oppositions, circle+equant (phi,e_s,e_q,w0,n free):')
    print('   period %.3f d (modern sidereal period 686.98 d); e_s=%.5f e_q=%.5f (Kepler quoted 0.11332, 0.07232); mean (e_s+e_q)/2=%.5f' % (2 * math.pi / nA, abs(pA[1]), abs(pA[2]), (abs(pA[1]) + abs(pA[2])) / 2))
    print('   rms residual %.2f arcmin, max %.2f arcmin' % (math.sqrt(np.mean(rA ** 2)), np.max(np.abs(rA))))
    # ---- B
    lon = np.degrees(opp_l) % 360; best = None
    for idx in itertools.combinations(range(12), 4):
        s = sorted(lon[list(idx)]); g = min((s[(i + 1) % 4] - s[i]) % 360 for i in range(4))
        if best is None or g > best[0]: best = (g, list(idx))
    idx = best[1]
    pB, nB, cB = fit_free(idx, opp_l[idx], fixed_n=nA)
    oth = [i for i in range(12) if i not in idx]
    rB = lon_res(pB, nB, oth)
    print('B. four most uniform oppositions %s (longitudes %s), n fixed from A:' % (idx, [round(lon[i]) for i in idx]))
    print('   e_s=%.5f e_q=%.5f; residuals at the other 8 oppositions (arcmin):' % (abs(pB[1]), abs(pB[2])), np.round(rB, 2).tolist(), ' max %.2f' % np.max(np.abs(rB)))
    # ---- C
    pC, nC, cC = fit_free(allidx, opp_l, bisected=True)
    rC = lon_res(pC, nC)
    print('C. bisected (e_s = e_q = e), all 12 oppositions:  e=%.5f  rms %.2f arcmin, max %.2f arcmin' % (abs(pC[1]), math.sqrt(np.mean(rC ** 2)), np.max(np.abs(rC))))
    print('   residuals (arcmin):', np.round(rC, 2).tolist())
    pC4, nC4, _ = fit_free(idx, opp_l[idx], bisected=True, fixed_n=nA)
    rC4 = lon_res(pC4, nC4, oth)
    print("C'. bisected, fitted to the same FOUR oppositions only (n from A): e=%.5f; residuals at the other 8 (arcmin):" % abs(pC4[1]), np.round(rC4, 2).tolist(), ' max %.2f' % np.max(np.abs(rC4)))
    # ---- D  distances: model radius R is set by one overall scale (best fit) since longitudes do not fix it
    def scaled_dist_err(par, nn):
        d = geom(par, opp_t, nn)[1]; s = np.sum(d * opp_r) / np.sum(d * d); return (s * d) / opp_r - 1.0
    dA, dC = scaled_dist_err(pA, nA), scaled_dist_err(pC, nC)
    print('D. Sun-Mars distance at the 12 oppositions, model vs ephemeris (relative error, one best scale each):')
    print('   unbisected (A): rms %.4f, max |err| %.4f ;  bisected (C): rms %.4f, max |err| %.4f' % (math.sqrt(np.mean(dA ** 2)), np.max(np.abs(dA)), math.sqrt(np.mean(dC ** 2)), np.max(np.abs(dC))))
    # ---- E  noise
    rng = np.random.default_rng(3); trials = 40; sig = 1.5; outA, outC = [], []
    for _ in range(trials):
        tgt = opp_l + rng.normal(0, sig / ARC, 12)
        pa, na, _ = fit_free(allidx, tgt); pc, nc, _ = fit_free(allidx, tgt, bisected=True)
        ra = wrap(geom(pa, opp_t, na)[0] - tgt) * ARC; rc = wrap(geom(pc, opp_t, nc)[0] - tgt) * ARC
        outA.append(math.sqrt(np.mean(ra ** 2))); outC.append(math.sqrt(np.mean(rc ** 2)))
    print('E. with %.1f arcmin Gaussian noise on the longitudes (%d trials): rms residual unbisected %.2f (sd %.2f), bisected %.2f (sd %.2f); bisected worse in %d/%d trials' %
          (sig, trials, np.mean(outA), np.std(outA), np.mean(outC), np.std(outC), sum(c > a for a, c in zip(outA, outC)), trials))