"""Verification of every identity and number quoted in the Lambert note (lambert-problem.qmd).

Part A is symbolic (sympy); Part B is numerical (both worked examples, the admissibility table,
monotonicity of the clock, multiple revolutions, the loss of precision near e = 1, the pinch at 180 deg).
Run from the sfm-slides directory:   python3 codes/lambert-verify.py        (about a minute)
"""
import warnings
import numpy as np
import sympy as sp
import mpmath as mp
from scipy.optimize import brentq
from _lambertp import mu, conic_through, p_window, time_between, tof_of_p, velocity

warnings.filterwarnings("ignore", category=RuntimeWarning)
ok = lambda label, cond: print(f"  [{'ok' if cond else 'FAIL'}] {label}")

# =================================================================================== A. symbolic
print("A. symbolic")
p, r1, r2 = sp.symbols("p r1 r2", positive=True)
C, S = sp.symbols("C S", real=True)                       # cos(dtheta), sin(dtheta); S^2 = 1 - C^2
A, B = p/r1 - 1, p/r2 - 1
ec = A
es = (A*C - B)/S
ok("e cos(f1 + dtheta) = B  (angle addition)", sp.simplify(ec*C - es*S - B) == 0)
ok("e sin(f1 + dtheta) = (A - B cos dtheta)/sin dtheta",
   sp.simplify(sp.expand(S*(es*C + ec*S) - (A - B*C)).subs(S**2, 1 - C**2)) == 0)

k, ell, m = r1*r2*(1 - C), r1 + r2, r1*r2*(1 + C)
c2 = r1**2 + r2**2 - 2*r1*r2*C
e2S2 = A**2 + B**2 - 2*A*B*C                              # e^2 sin^2(dtheta)
ok("e^2 sin^2 = A^2 + B^2 - 2AB cos", sp.simplify(ec**2 + es**2 - e2S2/S**2).subs(S**2, 1 - C**2) == 0 or
   sp.simplify(sp.expand((ec**2 + es**2)*S**2 - e2S2).subs(S**2, 1 - C**2)) == 0)
poly = sp.expand(e2S2*r1**2*r2**2 - (c2*p**2 - 2*k*ell*p + 2*k*r1*r2))
ok("r1^2 r2^2 (A^2+B^2-2AB cos) = c^2 p^2 - 2 k l p + 2 k r1 r2", sp.simplify(poly) == 0)
ok("k m = r1^2 r2^2 sin^2", sp.simplify((k*m - r1**2*r2**2*(1 - C**2))) == 0)
ok("c^2 = l^2 - 2m", sp.simplify(c2 - (ell**2 - 2*m)) == 0)
ok("e^2 - 1 = (c^2 p^2 - 2 k l p + k^2)/(r1 r2 sin)^2",
   sp.simplify(sp.expand((e2S2*r1**2*r2**2 - r1**2*r2**2*S**2) - (c2*p**2 - 2*k*ell*p + k**2)).subs(S**2, 1 - C**2)) == 0)
sq = sp.symbols("q", positive=True)                       # q = sqrt(2m)
pi_s, pii_s = k/(ell + sq), k/(ell - sq)
ok("p_i + p_ii = 2 k l/c^2 and p_i p_ii = k^2/c^2   (with c^2 = l^2 - q^2)",
   sp.simplify(pi_s + pii_s - 2*k*ell/(ell**2 - sq**2)) == 0 and sp.simplify(pi_s*pii_s - k**2/(ell**2 - sq**2)) == 0)
ok("e_min^2 = 1 - 2k/c^2 = (r1 - r2)^2/c^2",
   sp.simplify(1 - 2*k/c2 - (r1 - r2)**2/c2) == 0)

# half-angle relations (Section 4)
e_, f_ = sp.symbols("e f", positive=True)
cosE = (e_ + sp.cos(f_))/(1 + e_*sp.cos(f_))
a_ = sp.symbols("a", positive=True)
ok("elliptic: r = a(1 - e cos E) = p/(1 + e cos f)  with  cos E = (e + cos f)/(1 + e cos f), p = a(1-e^2)",
   sp.simplify(a_*(1 - e_*cosE) - a_*(1 - e_**2)/(1 + e_*sp.cos(f_))) == 0)
tE2 = (1 - cosE)/(1 + cosE)
ok("tan^2(E/2) = (1-e)/(1+e) tan^2(f/2)   [tan^2(f/2) = (1-cos f)/(1+cos f)]",
   sp.simplify(tE2 - (1 - e_)/(1 + e_)*(1 - sp.cos(f_))/(1 + sp.cos(f_))) == 0)
eh = sp.symbols("e_h", positive=True)                     # e > 1
cosh_H = (eh + sp.cos(f_))/(1 + eh*sp.cos(f_))
ok("hyperbolic: |a|(e cosh H - 1) = p/(1+e cos f), p = |a|(e^2-1)",
   sp.simplify(a_*(eh*cosh_H - 1) - a_*(eh**2 - 1)/(1 + eh*sp.cos(f_))) == 0)
ok("tanh^2(H/2) = (e-1)/(e+1) tan^2(f/2)",
   sp.simplify((cosh_H - 1)/(cosh_H + 1) - (eh - 1)/(eh + 1)*(1 - sp.cos(f_))/(1 + sp.cos(f_))) == 0)

# 180 degrees
ok("dtheta = pi: A + B = 0  =>  p = 2 r1 r2/(r1 + r2)", sp.simplify(sp.solve(sp.Eq(A + B, 0), p)[0] - 2*r1*r2/(r1 + r2)) == 0)

# Lagrange coefficients (Section 8): in the frame at r1
mu_s = sp.symbols("mu", positive=True)
th = sp.symbols("theta", real=True)                        # dtheta, cos = C, sin = S
R1 = sp.Matrix([r1, 0]); R2 = r2*sp.Matrix([C, S])
s0 = sp.sqrt(mu_s/p)
V1 = s0*sp.Matrix([es, 1 + ec])                            # (v_r, v_theta) at r1; r1 is along +x
esb, ecb = (A - B*C)/S, B
vr2, vt2 = s0*esb, s0*(1 + ecb)
V2 = sp.Matrix([vr2*C - vt2*S, vr2*S + vt2*C])
kap = 1 - C
F = 1 - r2*kap/p; G = r1*r2*S/sp.sqrt(mu_s*p); Gd = 1 - r1*kap/p
Fd = (F*Gd - 1)/G
sub = lambda x: sp.simplify(sp.expand(x.subs(S, sp.sqrt(1 - C**2))))
ok("r2 = F r1 + G v1", sub(R2 - F*R1 - G*V1) == sp.zeros(2, 1))
ok("v2 = Fdot r1 + Gdot v1", sub(V2 - Fd*R1 - Gd*V1) == sp.zeros(2, 1))
ok("v_r1/v_theta1 sin = cos - r1/r2 + (r1/p)(1 - cos)", sub(es*(r1/p)*S - (C - r1/r2 + (r1/p)*kap)) == 0)
ok("v_r2/v_theta2 sin = r2/r1 - cos - (r2/p)(1 - cos)", sub(esb*(r2/p)*S - (r2/r1 - C - (r2/p)*kap)) == 0)
ok("F Gdot - Fdot G = 1 (by construction); Fdot = sqrt(mu/p) tan(dtheta/2)[(1-cos)/p - 1/r1 - 1/r2]",
   sub(Fd - s0*(1 - C)/S*(kap/p - 1/r1 - 1/r2)) == 0)
ok("F Gdot - Fdot G computed from the definitions equals 1 with the independent Fdot",
   sub(F*Gd - s0*(1 - C)/S*(kap/p - 1/r1 - 1/r2)*G - 1) == 0)
ok("independent check: [r2, v2] = [r1, v1] (angular momentum)",
   sub((R2[0]*V2[1] - R2[1]*V2[0]) - (R1[0]*V1[1] - R1[1]*V1[0])) == 0)

# =================================================================================== B. numerical
print("B. numerical")
MU = mu
R1_, R2_ = 6678.137, 42164.0

def window(r1_, r2_, dth):
    k_, l_, m_ = r1_*r2_*(1 - np.cos(dth)), r1_ + r2_, r1_*r2_*(1 + np.cos(dth))
    return k_, l_, m_, *p_window(r1_, r2_, dth)

def secant(dt, p0, p1, r1_, r2_, dth, nmax=14):
    rows = [(p0, tof_of_p(p0, r1_, r2_, dth)), (p1, tof_of_p(p1, r1_, r2_, dth))]
    for _ in range(nmax):
        (pa, ta), (pb, tb) = rows[-2], rows[-1]
        pn = pb - (tb - dt)*(pb - pa)/(tb - ta)
        rows.append((pn, tof_of_p(pn, r1_, r2_, dth)))
        if abs(rows[-1][1] - dt) < 1e-6:
            break
    return rows

def lagr(pp, r1_, r2_, dth):
    kap_ = 1 - np.cos(dth)
    f_ = 1 - r2_*kap_/pp; g_ = r1_*r2_*np.sin(dth)/np.sqrt(MU*pp); gd_ = 1 - r1_*kap_/pp
    fd_ = np.sqrt(MU/pp)*np.tan(dth/2)*(kap_/pp - 1/r1_ - 1/r2_)
    return f_, g_, gd_, fd_

dth1 = np.radians(175.0)
k1, l1, m1, pi1, pii1 = window(R1_, R2_, dth1)
print(f"  Example 1: k={k1:.5e}, l={l1:.3f}, m={m1:.5e}, p_i={pi1:.3f}, p_ii={pii1:.3f}")
print(f"    parabola at p_ii: {tof_of_p(pii1*(1 - 1e-12), R1_, R2_, dth1)/3600:.4f} h")
print("    secant from 11300, 11800 (dt = 4 h):  step, p, dt[h], e")
rows = secant(4*3600.0, 11300.0, 11800.0, R1_, R2_, dth1)
for i, (pp, tt) in enumerate(rows):
    print(f"      {i}  {pp:11.4f}  {tt/3600:8.4f}  {conic_through(pp, R1_, R2_, dth1)[0]:.5f}")
p1 = brentq(lambda x: tof_of_p(x, R1_, R2_, dth1) - 4*3600.0, pi1*(1 + 1e-9), pii1*(1 - 1e-9), xtol=1e-12)
e1, f1 = conic_through(p1, R1_, R2_, dth1)
A1, B1 = p1/R1_ - 1, p1/R2_ - 1
a1 = p1/(1 - e1**2); n1 = np.sqrt(MU/a1**3)
Ef = lambda f: 2*np.arctan2(np.sqrt(1 - e1)*np.sin(f/2), np.sqrt(1 + e1)*np.cos(f/2))
E1, E2 = Ef(f1), Ef(f1 + dth1); M1, M2 = E1 - e1*np.sin(E1), E2 - e1*np.sin(E2)
print(f"    p={p1:.3f}, A={A1:.6f}, B={B1:.6f}, e sin f1={e1*np.sin(f1):.6f}, e={e1:.5f}, f1={np.degrees(f1):.3f}, f2={np.degrees(f1 + dth1):.3f} deg")
print(f"    a={a1:.2f}, n={n1:.6e}, E1={E1:.6f}, E2={E2:.6f}, M1={M1:.6f}, M2={M2:.6f}, (M2-M1)/n={(M2 - M1)/n1:.2f} s")
v1 = velocity(p1, e1, f1); v2 = velocity(p1, e1, f1 + dth1)
g1 = np.degrees(np.arctan2(v1[0], v1[1])); g2 = np.degrees(np.arctan2(v2[0], v2[1]))
print(f"    v1=(vr,vt)=({v1[0]:.5f}, {v1[1]:.5f}), v2=({v2[0]:.5f}, {v2[1]:.5f}) km/s; flight-path angles {g1:.2f}, {g2:.2f} deg")
vc1, vc2 = np.sqrt(MU/R1_), np.sqrt(MU/R2_)
dv1 = np.hypot(v1[0], v1[1] - vc1); dv2 = np.hypot(v2[0], v2[1] - vc2)
ah = 0.5*(R1_ + R2_); dvH1 = vc1*(np.sqrt(R2_/ah) - 1); dvH2 = vc2*(1 - np.sqrt(R1_/ah))
print(f"    dv1={dv1:.5f}, dv2={dv2:.5f}, total={dv1 + dv2:.5f} km/s;  Hohmann {dvH1 + dvH2:.5f} km/s in {np.pi*np.sqrt(ah**3/MU)/3600:.3f} h")
Av, Bv = A1, B1
vr1_c = np.sqrt(MU/p1)*(Av*np.cos(dth1) - Bv)/np.sin(dth1); vr2_c = np.sqrt(MU/p1)*(Av - Bv*np.cos(dth1))/np.sin(dth1)
vt1_c, vt2_c = np.sqrt(MU*p1)/R1_, np.sqrt(MU*p1)/R2_
ok("closed-form end velocities (no e, no f) equal the element-based ones",
   np.allclose([vr1_c, vt1_c, vr2_c, vt2_c], [v1[0], v1[1], v2[0], v2[1]], atol=1e-9))
fL, gL, gdL, fdL = lagr(p1, R1_, R2_, dth1)
print(f"    F={fL:.5f}, G={gL:.3f} s, Gdot={gdL:.5f}, Fdot={fdL:.4e};  F Gdot - Fdot G = {fL*gdL - fdL*gL:.9f}")
r1v = np.array([R1_, 0.0]); r2v = R2_*np.array([np.cos(dth1), np.sin(dth1)])
w1 = (r2v - fL*r1v)/gL; w2 = (gdL*r2v - r1v)/gL
print(f"    r2 = {np.round(r2v, 3).tolist()} km;  v1 = {np.round(w1, 5).tolist()};  v2 = {np.round(w2, 5).tolist()} km/s")
ok("Lagrange velocities equal the element-based ones (v1)", np.allclose(w1, [v1[0], v1[1]], atol=1e-9))
rot = np.array([[np.cos(dth1), -np.sin(dth1)], [np.sin(dth1), np.cos(dth1)]])
ok("... and v2", np.allclose(w2, rot @ np.array(v2), atol=1e-9))
print(f"    1 - cos(dtheta) = {1 - np.cos(dth1):.6f}; (r2 - F r1) = {np.round(r2v - fL*r1v, 3).tolist()}; (gd r2 - r1) = {np.round(gdL*r2v - r1v, 3).tolist()}")

# Lambert's theorem in Lagrange's form for the elliptic arc of Example 1
c_ = np.linalg.norm(r2v - r1v); s_ = 0.5*(R1_ + R2_ + c_)
al = 2*np.arcsin(np.sqrt(s_/(2*a1))); be = 2*np.arcsin(np.sqrt((s_ - c_)/(2*a1)))
tL = ((al - np.sin(al)) - (be - np.sin(be)))/n1
print(f"    Lambert's theorem: c={c_:.3f}, s={s_:.3f}, alpha={al:.6f}, beta={be:.6f}, dt={tL:.3f} s (target 14400)")

# ---- Example 2 (2 h): hyperbola
print("  Example 2 (dt = 2 h):  step, p, dt[h], e")
for i, (pp, tt) in enumerate(secant(2*3600.0, 11900.0, 12200.0, R1_, R2_, dth1)):
    print(f"      {i}  {pp:11.4f}  {tt/3600:8.4f}  {conic_through(pp, R1_, R2_, dth1)[0]:.5f}")
p2 = brentq(lambda x: tof_of_p(x, R1_, R2_, dth1) - 2*3600.0, pii1*(1 + 1e-9), 3*pii1, xtol=1e-12)
e2, f2a = conic_through(p2, R1_, R2_, dth1)
a2 = p2/(1 - e2**2); nh = np.sqrt(MU/(-a2)**3); finf = np.arccos(-1/e2)
Hf = lambda f: 2*np.arctanh(np.sqrt((e2 - 1)/(e2 + 1))*np.tan(f/2))
H1, H2 = Hf(f2a), Hf(f2a + dth1); Mh1, Mh2 = e2*np.sinh(H1) - H1, e2*np.sinh(H2) - H2
print(f"    p={p2:.3f}, e={e2:.5f}, f1={np.degrees(f2a):.3f}, f2={np.degrees(f2a + dth1):.3f}, f_inf={np.degrees(finf):.3f} deg")
print(f"    a={a2:.2f}, n_h={nh:.6e}, H1={H1:.6f}, H2={H2:.6f}, Mh1={Mh1:.6f}, Mh2={Mh2:.6f}, (Mh2-Mh1)/n_h={(Mh2 - Mh1)/nh:.2f} s")
u1 = velocity(p2, e2, f2a); u2 = velocity(p2, e2, f2a + dth1)
du1 = np.hypot(u1[0], u1[1] - vc1); du2 = np.hypot(u2[0], u2[1] - vc2)
print(f"    v1=({u1[0]:.5f}, {u1[1]:.5f}), v2=({u2[0]:.5f}, {u2[1]:.5f}); dv1={du1:.5f}, dv2={du2:.5f}, total={du1 + du2:.5f} km/s")

# ---- the parabola (Barker): continuity of the clock across e = 1
def barker(pp, r1_, r2_, dth, short=True):
    """Time on the parabola of semi-latus rectum pp between the two points (e = 1, |f| < pi)."""
    A_, B_ = pp/r1_ - 1, pp/r2_ - 1
    f1_ = np.arctan2((A_*np.cos(dth) - B_)/np.sin(dth), A_)
    f2_ = f1_ + dth
    Dn = lambda f: np.tan(f/2) + np.tan(f/2)**3/3
    return 0.5*np.sqrt(pp**3/MU)*(Dn(f2_) - Dn(f1_))
tb = barker(pii1, R1_, R2_, dth1)
print(f"  Barker's equation at p_ii (Example 1): {tb/3600:.5f} h; ellipse side {tof_of_p(pii1*(1 - 1e-8), R1_, R2_, dth1)/3600:.5f} h, hyperbola side {tof_of_p(pii1*(1 + 1e-8), R1_, R2_, dth1)/3600:.5f} h")
ok("the clock is continuous across the parabola (Barker = limit from both sides, 5e-6; the clock is steep here, dln(dt)/dln(p) = -22)",
   abs(tb - tof_of_p(pii1*(1 - 1e-8), R1_, R2_, dth1)) < 5e-6*tb and abs(tb - tof_of_p(pii1*(1 + 1e-8), R1_, R2_, dth1)) < 5e-6*tb)

# ---- which parabola is admissible for which way (Section 3 / the window figure)
dthL = np.radians(260.0); RAx, RBx = 7000.0, 11000.0
for label, dthx, p_ok, p_bad in (("short way, 100 deg", np.radians(100.0), "ii", "i"), ("long way, 260 deg", dthL, "i", "ii")):
    pi_x, pii_x = p_window(RAx, RBx, dthx)
    pk = {"i": pi_x, "ii": pii_x}
    A_, B_ = pk[p_ok]/RAx - 1, pk[p_ok]/RBx - 1
    f1x = np.arctan2((A_*np.cos(dthx) - B_)/np.sin(dthx), A_)
    f2x = f1x + dthx
    t_ok = barker(pk[p_ok], RAx, RBx, dthx)
    side = 1 - 1e-8 if p_ok == "ii" else 1 + 1e-8
    t_nb = tof_of_p(pk[p_ok]*side, RAx, RBx, dthx)
    sidebad = 1 + 1e-8 if p_bad == "i" else 1 - 1e-8
    t_bad = tof_of_p(pk[p_bad]*sidebad, RAx, RBx, dthx)
    print(f"    {label}: parabola at p_{p_ok}: f1={np.degrees(f1x):.1f}, f2={np.degrees(f2x):.1f} deg (< 180), Barker {t_ok/60:.3f} min, neighbour {t_nb/60:.3f} min;  next to p_{p_bad}: dt = {t_bad/3600:.2e} h")
    ok(f"{label}: p_{p_ok} is an admissible parabola (finite time, f2 < pi); dt diverges towards p_{p_bad}",
       f2x < np.pi and abs(t_ok - t_nb) < 1e-5*t_ok and t_bad > 1e3*t_ok)

# ---- Euler's parabolic time: (1/3) sqrt(2/mu) [ s^(3/2) -/+ (s - c)^(3/2) ]  (minus: short way, plus: long way)
def euler_par(r1_, r2_, dth):
    c_ = np.sqrt(r1_**2 + r2_**2 - 2*r1_*r2_*np.cos(dth)); s_ = 0.5*(r1_ + r2_ + c_)
    sgn = 1.0 if np.sin(dth) > 0 else -1.0
    return np.sqrt(2/MU)/3*(s_**1.5 - sgn*(s_ - c_)**1.5), c_, s_
good_par = True
for r1x, r2x, dx in ((R1_, R2_, 175.0), (R1_, R2_, 185.0), (7000.0, 11000.0, 100.0), (7000.0, 11000.0, 260.0), (5000.0, 30000.0, 40.0), (5000.0, 30000.0, 300.0)):
    dthx = np.radians(dx); pi_x, pii_x = p_window(r1x, r2x, dthx)
    pp_ = pii_x if dx < 180 else pi_x
    tb_ = barker(pp_, r1x, r2x, dthx); te_, c_x, s_x = euler_par(r1x, r2x, dthx)
    good_par &= abs(tb_ - te_) < 1e-9*tb_
    print(f"    Euler's parabolic time, r1={r1x:.0f}, r2={r2x:.0f}, dtheta={dx}: {te_/3600:.5f} h   (Barker at the admissible parabola: {tb_/3600:.5f} h)")
ok("Euler's parabolic time equals Barker's equation at the admissible parabola (both ways)", good_par)

# ---- admissibility table (Section 3): scan random geometries
print("  admissibility of p (zero revolutions), random geometries:")
rng = np.random.default_rng(1)
bad = 0; nsc = 0
for _ in range(400):
    r1r, r2r = rng.uniform(1, 10, 2)
    dth = rng.uniform(0.15, 2*np.pi - 0.15)
    if abs(dth - np.pi) < 0.08:
        continue
    pi_r, pii_r = p_window(r1r, r2r, dth)
    for pp in np.exp(rng.uniform(np.log(1e-2), np.log(1e3), 40)):
        if min(abs(pp - pi_r), abs(pp - pii_r)) < 1e-6*pp:
            continue
        e_, f1_ = conic_through(pp, r1r, r2r, dth)
        if e_ < 1:
            admiss = True
        else:
            finf_ = np.arccos(-1/e_)
            admiss = (-finf_ < f1_ < finf_) and (-finf_ < f1_ + dth < finf_)
        predicted = (pp > pi_r) if dth < np.pi else (pp < pii_r)
        bad += (admiss != predicted); nsc += 1
ok(f"window table matches the direct test (unwrapped anomalies inside (-f_inf, f_inf)) on {nsc} samples", bad == 0)

# ---- monotonicity of the clock (zero revolutions)
print("  monotonicity of dt(p), zero revolutions, random geometries:")
viol = 0; ncase = 0
for _ in range(120):
    r1r, r2r = rng.uniform(1, 10, 2)*6000
    dth = rng.uniform(0.2, 2*np.pi - 0.2)
    if abs(dth - np.pi) < 0.1:
        continue
    pi_r, pii_r = p_window(r1r, r2r, dth)
    if dth < np.pi:
        ps = np.concatenate([np.geomspace(pi_r*(1 + 1e-6), pii_r*(1 - 1e-6), 200), np.geomspace(pii_r*(1 + 1e-6), 1e3*pii_r, 200)])
        sign = -1
    else:
        ps = np.concatenate([np.geomspace(1e-4*pi_r, pi_r*(1 - 1e-6), 200), np.geomspace(pi_r*(1 + 1e-6), pii_r*(1 - 1e-6), 200)])
        sign = +1
    ts = np.array([tof_of_p(q, r1r, r2r, dth) for q in ps])
    good = np.isfinite(ts)
    d = np.diff(ts[good])*sign
    viol += int(np.any(d < -1e-9*np.abs(ts[good][1:]))); ncase += 1
ok(f"dt falls with p (short way) / rises with p (long way) in all {ncase} geometries", viol == 0)
pl = np.geomspace(1e2*pii1, 1e6*pii1, 5)
tl = np.array([tof_of_p(q, R1_, R2_, dth1) for q in pl])
sl = np.diff(np.log(tl))/np.diff(np.log(pl))
print(f"    short way, p -> infinity: dt = {tl[0]/3600:.3e}, {tl[-1]/3600:.3e} h; log-log slope {sl[-1]:.3f} (expect -1/2)")
dth_l = np.radians(185.0); pi_l, pii_l = p_window(R1_, R2_, dth_l)
pl2 = np.geomspace(1e-6*pi_l, 1e-3*pi_l, 4); tl2 = np.array([tof_of_p(q, R1_, R2_, dth_l) for q in pl2])
print(f"    long way, p -> 0: dt = {tl2[0]/3600:.3e} ... {tl2[-1]/3600:.3e} h; slope {np.diff(np.log(tl2))[-1]/np.diff(np.log(pl2))[-1]:.3f}")
print(f"    dt at p_i(1+1e-9) for the short way: {tof_of_p(pi1*(1 + 1e-9), R1_, R2_, dth1)/3600:.3e} h  (diverges as p -> p_i)")

# ---- multiple revolutions: the four-panel geometry
print("  multiple revolutions, r1=7000 km, r2=11000 km (r2/r1 = 11/7):")
RA, RB = 7000.0, 11000.0
T1c = 2*np.pi*np.sqrt(RA**3/MU)
def solve_all(deg, N, dt):
    dth = np.radians(deg); pi_r, pii_r = p_window(RA, RB, dth)
    ps = np.linspace(pi_r*(1 + 1e-9), pii_r*(1 - 1e-9), 30001)
    ts = np.array([tof_of_p(q, RA, RB, dth, N) for q in ps]) - dt
    out = []
    for i in range(len(ps) - 1):
        if np.isfinite(ts[i]) and np.isfinite(ts[i + 1]) and ts[i]*ts[i + 1] < 0:
            out.append(brentq(lambda x: tof_of_p(x, RA, RB, dth, N) - dt, ps[i], ps[i + 1], xtol=1e-9))
    return out
for deg in (100.0, 260.0):
    for N in (0, 1, 2):
        dth = np.radians(deg); pi_r, pii_r = p_window(RA, RB, dth)
        ps = np.linspace(pi_r*(1 + 1e-9), pii_r*(1 - 1e-9), 30001)
        ts = np.array([tof_of_p(q, RA, RB, dth, N) for q in ps])
        print(f"    dtheta={deg:5.1f}, N={N}: min dt = {np.nanmin(ts)/60:8.2f} min ({np.nanmin(ts)/T1c:.3f} T1) at p={ps[np.nanargmin(ts)]:.1f}")
DT = 2.0*T1c
print(f"    dt = 2 T1 = {DT/60:.2f} min; solutions:")
cnt = {}
for deg in (100.0, 260.0):
    for N in (0, 1):
        sols = solve_all(deg, N, DT); cnt[(deg, N)] = len(sols)
        for q in sols:
            e_, f1_ = conic_through(q, RA, RB, np.radians(deg))
            print(f"      dtheta={deg:5.1f} N={N}: p={q:9.2f}, e={e_:.4f}, a={q/(1 - e_**2):9.1f}, f1={np.degrees(f1_):7.2f} deg")
ok("one solution for N = 0 (both ways), two for N = 1 (both ways)", [cnt[(100.0, 0)], cnt[(260.0, 0)], cnt[(100.0, 1)], cnt[(260.0, 1)]] == [1, 1, 2, 2])
ok("no N = 1 solution below the minimum time", len(solve_all(100.0, 1, 0.9*np.nanmin([tof_of_p(q, RA, RB, np.radians(100.0), 1) for q in np.linspace(p_window(RA, RB, np.radians(100.0))[0]*(1 + 1e-9), p_window(RA, RB, np.radians(100.0))[1]*(1 - 1e-9), 3000)]))) == 0)

# ---- precision loss near e = 1 (double against 50 digits)
print("  precision of the elliptic clock just below p_ii (Example 1 geometry):  1-e, relative error of dt, eps/(1-e)")
mp.mp.dps = 50
MPMU = mp.mpf("398600.4418")
def clock_mp(pp, r1_, r2_, dth):
    pp, r1_, r2_, dth = mp.mpf(pp), mp.mpf(r1_), mp.mpf(r2_), mp.mpf(dth)
    A_, B_ = pp/r1_ - 1, pp/r2_ - 1
    ecm, esm = A_, (A_*mp.cos(dth) - B_)/mp.sin(dth)
    em = mp.sqrt(ecm**2 + esm**2); f1m = mp.atan2(esm, ecm); f2m = f1m + dth
    am = pp/(1 - em*em)
    Em = lambda f: 2*mp.atan2(mp.sqrt(1 - em)*mp.sin(f/2), mp.sqrt(1 + em)*mp.cos(f/2))
    E1m, E2m = Em(f1m), Em(f2m)
    dE = (E2m - E1m) % (2*mp.pi)
    dM = dE - em*(mp.sin(E1m + dE) - mp.sin(E1m))
    return dM*mp.sqrt(am**3/MPMU), em
for target in (1.6e-5, 1.6e-8, 1.6e-11, 1.6e-13):
    lo, hi = mp.mpf(pii1)*(1 - mp.mpf("1e-3")), mp.mpf(pii1)
    for _ in range(200):
        mid = (lo + hi)/2
        if 1 - clock_mp(mid, R1_, R2_, dth1)[1] > target:
            lo = mid
        else:
            hi = mid
    pd = float(mid)
    td = tof_of_p(pd, R1_, R2_, dth1); tm, em_ = clock_mp(pd, R1_, R2_, dth1)
    print(f"      1-e = {1 - float(em_):.2e}:  relative error = {abs(td - float(tm))/float(tm):.1e},  eps/(1-e) = {np.finfo(float).eps/(1 - float(em_)):.1e}")

# ---- the pinch at 180 degrees
print("  window width near 180 degrees (Example 1 radii):  numeric against 4 (r1 r2)^{3/2} (pi - dtheta)/l^2")
for deg in (179.9, 179.0, 175.0):
    dth = np.radians(deg); _, _, _, a_i, a_ii = window(R1_, R2_, dth)
    print(f"      {deg}: p_ii - p_i = {a_ii - a_i:.3f} km;  formula {4*(R1_*R2_)**1.5*(np.pi - dth)/(R1_ + R2_)**2:.3f} km;  p_H = {2*R1_*R2_/(R1_ + R2_):.3f}")
dth = np.radians(179.9); pi_, pii_ = p_window(R1_, R2_, dth)
t_lo = tof_of_p(pi_ + 0.05*(pii_ - pi_), R1_, R2_, dth)/3600; t_hi = tof_of_p(pi_ + 0.95*(pii_ - pi_), R1_, R2_, dth)/3600
print(f"      at 179.9 deg the window ({pii_ - pi_:.2f} km) sweeps dt from {t_lo:.2f} h to {t_hi:.2f} h")
pa = brentq(lambda x: tof_of_p(x, R1_, R2_, dth) - 5*3600.0, pi_*(1 + 1e-9), pii_*(1 - 1e-9), xtol=1e-12)
pb = brentq(lambda x: tof_of_p(x, R1_, R2_, dth) - 1.8*3600.0, pii_*(1 + 1e-9), pii_ + 100.0, xtol=1e-12)
print(f"      at 179.9 deg: dt = 5 h needs p = {pa:.3f} km, dt = 1.8 h needs p = {pb:.3f} km: a change of {pb - pa:.2f} km")