#!/usr/bin/env python3
"""
test_350 -- THE BIRTH OF A FOLD OF THE GAUSS MAP: THE CUBIC NORMAL FORM, THE EXACT EXCESS OF ABSOLUTE CURVATURE, AND
            THE delta^{3/2} LAW

For MB §6.2 (the companion of MA). The family is the embedded radial graphs X = r(chi) n, r = 1 + eps cos 2chi, chi the
polar angle on S^3. The outward normal of the meridian (r cos chi, r sin chi) makes the angle beta(chi) with the axis; by
SO(3)-symmetry beta determines the Gauss map. delta = eps - 1/5, x = chi - pi/2.

PRE-REGISTERED:
  P1 (the Gauss map, for every r) beta = chi - arctan(r'/r); beta' = (r^2 + 2r'^2 - r r'')/(r^2 + r'^2) = k_mer ds/dchi;
     k_sph = sin(beta)/(r sin chi); hence det K dV = 4 pi sin^2(beta) beta' dchi after integrating over the orbit
     spheres. For r(pi - chi) = r(chi): beta(pi - chi) = pi - beta(chi). (Symbolic, r arbitrary; the symmetry for the
     family.)
  P2 (the normal form, exactly) beta(pi/2 + x) - pi/2 = a1 x + a3 x^3 + O(x^5), with no even powers,
     a1 = (1 - 5eps)/(1 - eps), a3 = 8 eps (6eps^2 + eps + 1)/(3 (1 - eps)^3); so a1 = -(25/4) delta - (125/16) delta^2
     + O(delta^3) and a3 = 3/2 at eps = 1/5. Whitney's cusp conditions hold for (x, delta) -> (beta, delta) at the origin:
     d beta/dx = d^2 beta/dx^2 = 0, d^3 beta/dx^3 = 9, d^2 beta/dx d delta = -25/4.
  P3 (the excess, exactly) For eps > 1/5 let a be the fold point in (0, pi/2) (beta'(a) = 0). Then
     J_abs - 2 pi^2 = 8 pi [G(beta(a)) - G(beta(pi - a))], G(b) = b/2 - sin(2b)/4, equal to 2 int_{det K < 0} |det K| dV
     and to the direct quadrature of 4 pi int sin^2 beta |beta'| dchi minus 2 pi^2, to 1e-25 at eps = 1/4, 3/10, 2/5;
     J_abs/2pi^2 = 1.1237308, 1.3137532, 1.7261913 there (the values in the author's archive, to 1e-6); and J_abs = 2 pi^2 (to 1e-25) at
     eps = 1/10 and 1/5, where beta' >= 0.
  P4 (the scaling) At delta = 1e-3, 1e-4, 1e-5, with x0 = pi/2 - a and Dbeta = beta(a) - beta(pi - a):
     x0 / ((5/sqrt18) delta^{1/2}), Dbeta / ((125/(9 sqrt2)) delta^{3/2}) and
     (J_abs/2pi^2 - 1) / ((500/(9 sqrt2 pi)) delta^{3/2}) each lie within 3 delta of 1.

FIRST RUN, AND THE CHANGES:
  P1 failed in simplification only. sqrt(1 + r'^2/r^2) = A/r needs r > 0, which sympy cannot assume of an undefined
  function, so the identities for k_sph and det K dV did not reduce. Re-posed: those two are checked on the jet
  (r, r', r'') = (R0, R1, R2) with R0 > 0, after differentiation. The claim is unchanged.
  P4 failed as registered. The x0 ratio at delta = 1e-3 is 0.99681, outside 3 delta: its relative correction is about
  -3.19 delta, the other two about -2.64 delta, and the bound had been set from the J ratio alone. Re-posed: each ratio
  within 4 delta of 1, and each deviation falls by a factor between 9 and 11 per decade of delta (an O(delta)
  correction, as the law asserts). P2 and P3 passed.
  Before the bundle's final run, a reference to the archive's working notes was reworded; no code changed.
"""
import sys
import sympy as sp
import mpmath as mp

FAILURES = []


def check(name, ok, detail=""):
    print(f"  [{'PASS' if ok else 'FAIL'}] {name}" + (f"  -- {detail}" if detail else ""))
    if not ok:
        FAILURES.append(name)


# ------------------------------------------------------------------------------------------------------------- P1
print("P1  the Gauss map")
ch, e, x, d = sp.symbols("chi epsilon x delta", real=True)
rf = sp.Function("r")(ch)
r1, r2 = sp.diff(rf, ch), sp.diff(rf, ch, 2)
beta_g = ch - sp.atan(r1 / rf)
A = sp.sqrt(rf**2 + r1**2)
kmer = (rf**2 + 2 * r1**2 - rf * r2) / A**3
ksph = (rf - r1 * sp.cot(ch)) / (rf * A)
ok1 = sp.simplify(sp.diff(beta_g, ch) - (rf**2 + 2 * r1**2 - rf * r2) / (rf**2 + r1**2)) == 0
ok1 &= sp.simplify(sp.diff(beta_g, ch) - kmer * A) == 0
# the next two on the jet (r, r', r'') = (R0, R1, R2) with R0 > 0, after differentiating: sqrt(1 + r'^2/r^2) = A/r
# needs r > 0, which sympy cannot assume of an undefined function
R0 = sp.Symbol("R0", positive=True)
R1, R2 = sp.symbols("R1 R2", real=True)
jet = lambda expr: expr.subs({r2: R2}).subs({r1: R1}).subs({rf: R0})
sinb = jet(sp.expand_trig(sp.sin(beta_g)))
ok1 &= sp.simplify(sinb / (R0 * sp.sin(ch)) - jet(ksph)) == 0
# det K dV over the orbit spheres: k_mer k_sph^2 (r sin chi)^2 A dchi * 4 pi = 4 pi sin^2 beta beta' dchi
ok1 &= sp.simplify(jet(kmer * ksph**2 * (rf * sp.sin(ch))**2 * A) - sinb**2 * jet(sp.diff(beta_g, ch))) == 0
r = 1 + e * sp.cos(2 * ch)
beta = ch - sp.atan(sp.diff(r, ch) / r)
ok1 &= sp.simplify(beta.subs(ch, sp.pi - ch) - (sp.pi - beta)) == 0
check("beta' = (r^2+2r'^2-rr'')/(r^2+r'^2) = k_mer ds/dchi; k_sph = sin beta/(r sin chi); det K dV = 4pi sin^2 beta dbeta; "
      "beta(pi - chi) = pi - beta(chi)", ok1)

# ------------------------------------------------------------------------------------------------------------- P2
print("P2  the cubic normal form")
ser = sp.expand(sp.series((beta - sp.pi / 2).subs(ch, sp.pi / 2 + x), x, 0, 6).removeO())
a1, a3 = sp.simplify(ser.coeff(x, 1)), sp.simplify(ser.coeff(x, 3))
ok2 = all(sp.simplify(ser.coeff(x, k)) == 0 for k in (0, 2, 4))
ok2 &= sp.simplify(a1 - (1 - 5 * e) / (1 - e)) == 0
ok2 &= sp.simplify(a3 - 8 * e * (6 * e**2 + e + 1) / (3 * (1 - e)**3)) == 0
ok2 &= sp.expand(sp.series(a1.subs(e, sp.Rational(1, 5) + d), d, 0, 3).removeO()
                 - (-sp.Rational(25, 4) * d - sp.Rational(125, 16) * d**2)) == 0
ok2 &= a3.subs(e, sp.Rational(1, 5)) == sp.Rational(3, 2)
bd = (beta - sp.pi / 2).subs({ch: sp.pi / 2 + x, e: sp.Rational(1, 5) + d})
at0 = {x: 0, d: 0}
ok2 &= sp.simplify(sp.diff(bd, x).subs(at0)) == 0 and sp.simplify(sp.diff(bd, x, 2).subs(at0)) == 0
ok2 &= sp.simplify(sp.diff(bd, x, 3).subs(at0)) == 9 and sp.simplify(sp.diff(bd, x, d).subs(at0)) == sp.Rational(-25, 4)
check("beta - pi/2 = a1 x + a3 x^3 + O(x^5), a1 = (1-5eps)/(1-eps), a3 = 8eps(6eps^2+eps+1)/(3(1-eps)^3); cusp conditions",
      ok2, f"a1 = {sp.factor(a1)}, a3 = {sp.factor(a3)}")

# ------------------------------------------------------------------------------------------------------------- P3
print("P3  the excess of absolute curvature, exactly")
mp.mp.dps = 40
G = lambda b: b / 2 - mp.sin(2 * b) / 4


def family(ev):
    rr = lambda c: 1 + ev * mp.cos(2 * c)
    r_1 = lambda c: -2 * ev * mp.sin(2 * c)
    r_2 = lambda c: -4 * ev * mp.cos(2 * c)
    bp = lambda c: (rr(c)**2 + 2 * r_1(c)**2 - rr(c) * r_2(c)) / (rr(c)**2 + r_1(c)**2)
    bb = lambda c: c - mp.atan(r_1(c) / rr(c))
    return bp, bb


def fold_point(bp, lo, hi):
    for _ in range(160):                      # bp(lo) > 0 > bp(hi)
        mid = (lo + hi) / 2
        lo, hi = (mid, hi) if bp(mid) > 0 else (lo, mid)
    return (lo + hi) / 2


ok3 = True
rows = []
want = {"1/4": mp.mpf("1.1237308"), "3/10": mp.mpf("1.3137532"), "2/5": mp.mpf("1.7261913")}
for lab, ev in (("1/4", mp.mpf(1) / 4), ("3/10", mp.mpf(3) / 10), ("2/5", mp.mpf(2) / 5)):
    bp, bb = family(ev)
    a = fold_point(bp, mp.mpf("0.5"), mp.pi / 2 - mp.mpf("1e-15"))
    exact = 8 * mp.pi * (G(bb(a)) - G(bb(mp.pi - a)))
    neg = 4 * mp.pi * mp.quad(lambda c: mp.sin(bb(c))**2 * abs(bp(c)), [a, mp.pi / 2, mp.pi - a])
    tot = 4 * mp.pi * mp.quad(lambda c: mp.sin(bb(c))**2 * abs(bp(c)), [0, a, mp.pi / 2, mp.pi - a, mp.pi])
    ok3 &= abs(exact - 2 * neg) < mp.mpf("1e-25") and abs(exact - (tot - 2 * mp.pi**2)) < mp.mpf("1e-25")
    ok3 &= abs(tot / (2 * mp.pi**2) - want[lab]) < mp.mpf("1e-6")
    rows.append(f"eps={lab}: fold at chi = {mp.nstr(a, 10)}, J_abs/2pi^2 = {mp.nstr(tot / (2 * mp.pi**2), 10)}")
for ev in (mp.mpf(1) / 10, mp.mpf(1) / 5):
    bp, bb = family(ev)
    tot = 4 * mp.pi * mp.quad(lambda c: mp.sin(bb(c))**2 * abs(bp(c)), [0, mp.pi / 2, mp.pi])
    ok3 &= abs(tot - 2 * mp.pi**2) < mp.mpf("1e-25") and bp(mp.pi / 2) >= -mp.mpf("1e-35")
check("J_abs - 2pi^2 = 8pi[G(beta(a)) - G(beta(pi-a))] = 2 int_{det K<0}|det K| = quadrature; the archive's values; "
      "2pi^2 at eps <= 1/5", ok3, "; ".join(rows))

# ------------------------------------------------------------------------------------------------------------- P4
print("P4  the scaling laws")
cx = 5 / mp.sqrt(18)
cb = 125 / (9 * mp.sqrt(2))
cJ = 500 / (9 * mp.sqrt(2) * mp.pi)
ok4 = True
rows = []
devs = []
for dv in (mp.mpf("1e-3"), mp.mpf("1e-4"), mp.mpf("1e-5")):
    bp, bb = family(mp.mpf(1) / 5 + dv)
    xh = cx * mp.sqrt(dv)
    a = fold_point(bp, mp.pi / 2 - 2 * xh, mp.pi / 2 - xh / 2)
    x0 = mp.pi / 2 - a
    Db = bb(a) - bb(mp.pi - a)
    Jx = 8 * mp.pi * (G(bb(a)) - G(bb(mp.pi - a))) / (2 * mp.pi**2)
    q = [x0 / (cx * mp.sqrt(dv)), Db / (cb * dv**1.5), Jx / (cJ * dv**1.5)]
    ok4 &= all(abs(v - 1) < 4 * dv for v in q)
    devs.append([abs(v - 1) for v in q])
    rows.append(f"delta={mp.nstr(dv, 2)}: " + ", ".join(mp.nstr(v, 10) for v in q))
for k in range(2):
    ok4 &= all(9 < devs[k][j] / devs[k + 1][j] < 11 for j in range(3))
check("x0 ~ (5/sqrt18) delta^1/2, Dbeta ~ (125/(9sqrt2)) delta^3/2, J_abs/2pi^2 - 1 ~ (500/(9sqrt2 pi)) delta^3/2, "
      "each within 4 delta, deviations falling tenfold per decade", ok4, "; ".join(rows))

print()
if FAILURES:
    print(f"test_350: {len(FAILURES)} FAILED")
    sys.exit(1)
print("test_350: ALL PASS")