"""Reconstruction 001: Michelson-Morley 1887, from the printed tables.
Method (stated before looking at results):
 1. data: 17 readings/row (marks 16,1..16; first and last are the same mark after one full turn), divisions*0.02 = fringes.
 2. remove linear drift using the end points (x17-x1 over 16 steps); 16 points remain, theta_k = 22.5deg * k.
 3. least-squares/DFT harmonic amplitudes h=1..8 in fringes. A rotation signal from ether drift has period 180deg = harmonic 2.
 4. noise floor: the other harmonics (3..8) carry no predicted signal; their RMS is the empirical noise scale.
 5. expected amplitude of harmonic 2 for v = orbital speed: A = D v^2/(lambda c^2), D = 11 m, lambda = 5.9e-7 m (sodium)  [peak-to-peak 2A is the 'displacement on 90 deg rotation' the paper quotes as 0.4].
Run: python3 mm1887.py
"""
import json, math
d = json.load(open('data/mm1887.json'))
c, v = 299792458.0, 29780.0
D, lam = 11.0, 5.9e-7
A_exp = D * v * v / (lam * c * c)
print('expected harmonic-2 amplitude (orbital v=29.78 km/s): %.3f fringe (peak-to-peak %.3f; paper says 0.4)' % (A_exp, 2 * A_exp))

def detrend(x):
    drift = (x[16] - x[0]) / 16.0
    return [x[k] - drift * k for k in range(16)]

def harm(y, h):
    n = len(y); m = sum(y) / n
    a = sum((y[k] - m) * math.cos(2 * math.pi * h * k / n) for k in range(n)) * 2 / n
    b = sum((y[k] - m) * math.sin(2 * math.pi * h * k / n) for k in range(n)) * 2 / n
    return math.hypot(a, b), math.degrees(math.atan2(b, a))

def report(label, x_div):
    y = detrend([u * 0.02 for u in x_div])
    amps = [harm(y, h)[0] for h in range(1, 9)]
    noise = math.sqrt(sum(a * a for a in amps[2:]) / len(amps[2:]))
    print('%-14s h1=%.3f h2=%.3f | noise(rms h3-8)=%.3f | h2/expected=%.3f' % (label, amps[0], amps[1], noise, amps[1] / A_exp))
    return amps[1], noise

res = {}
for sess in ('noon', 'evening'):
    rows = d[sess]; days = list(rows)
    out = [report(sess + ' ' + k, rows[k]) for k in days]
    mean = [sum(rows[k][i] for k in days) / len(days) for i in range(17)]
    res[sess] = report(sess + ' MEAN', mean)
h2 = [res['noon'][0], res['evening'][0]]
print('\nBest case for a signal: noon and evening mean harmonic-2 = %.3f and %.3f fringe' % tuple(h2))
for a in h2:
    print('  -> implied ether-wind speed if ALL of it were signal: %.1f km/s (orbital 29.8; paper says < ~5-7.5)' % (v * math.sqrt(a / A_exp) / 1000))

# --- phase consistency: a real drift fixed to the lab frame at a given hour would give the same phase of harmonic 2 on different days
print('\nPhase of harmonic 2 (degrees, of the 180-deg cycle) by day:')
for sess in ('noon', 'evening'):
    ph = []
    for k, x in d[sess].items():
        y = detrend([u * 0.02 for u in x]); a, p = harm(y, 2); ph.append(p); print('  %-8s %-8s amp %.3f phase %6.1f' % (sess, k, a, p))
    # resultant length R of unit phasors (R=1: perfectly consistent; ~0.5-0.6 expected for 3 random phases)
    R = math.hypot(sum(math.cos(math.radians(p)) for p in ph), sum(math.sin(math.radians(p)) for p in ph)) / len(ph)
    print('  %s: phase concentration R = %.2f' % (sess, R))
import random
random.seed(1); trials = 200000
Rs = []
for _ in range(trials):
    ph = [random.uniform(0, 2 * math.pi) for _ in range(3)]
    Rs.append(math.hypot(sum(math.cos(p) for p in ph), sum(math.sin(p) for p in ph)) / 3)
for sess in ('noon', 'evening'):
    ph = [harm(detrend([u * 0.02 for u in x]), 2)[1] for x in d[sess].values()]
    R = math.hypot(sum(math.cos(math.radians(p)) for p in ph), sum(math.sin(math.radians(p)) for p in ph)) / 3
    print('  %s: P(R >= %.2f | random phases, 3 days) = %.3f' % (sess, R, sum(r >= R for r in Rs) / trials))