#!/usr/bin/env python3
"""LUMEN Second Edition -- Substrate chapter verifier.

Reproduces, by direct numerical integration of the cubic phase density
    rho(x) = 16 pi^3 x^3 + 3 pi^2 x^2 + 2 pi x   on x in [0,1],
the fine-structure inverse, the three layer fractions, the self-lensing energy E_self, the
moment ratio MU that sets the mass hierarchy, and the non-collapsibility of the three layers.

Method: trapezoidal integration on a 1,000,001-point grid (numpy only); cross-checked against
the closed forms 4pi^3+pi^2+pi, etc. No network, no corpus dependency, deterministic.
Independently pinned in the First Edition by verify_P001.py, verify_P003.py, verify_P004.py."""
import sys, os
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import numpy as np
from numpy import pi
from se_verify_common import Report

trap = getattr(np, "trapezoid", np.trapz)
x = np.linspace(0, 1, 1_000_001)
I = lambda f: float(trap(f, x))
rho  = 16*pi**3*x**3 + 3*pi**2*x**2 + 2*pi*x
rhop = 48*pi**3*x**2 + 6*pi**2*x + 2*pi

R = Report("Substrate",
    "the fine-structure inverse, the three layer fractions, and the self-lensing energy, from one cubic density",
    "trapezoidal integration of rho(x) on [0,1] (numpy), cross-checked against closed forms",
    sources=["CODATA-2022 (alpha^-1 = 137.035999177)"],
    pins=["verify_P001.py", "verify_P003.py", "verify_P004.py"])

ai = I(rho)
R.check("integral of rho  =  alpha^-1", ai, 137.035999177, source="CODATA-2022", tol=1e-4,
        note="seed identity 4pi^3+pi^2+pi; the +2.2 ppm offset is the flagged seed gap")
bulk, bnd, edge = I(16*pi**3*x**3), I(3*pi**2*x**2), I(2*pi*x)
R.check("bulk layer  4pi^3", bulk, 4*pi**3, source="closed form", tol=1e-3)
R.check("boundary layer  pi^2", bnd, pi**2, source="closed form", tol=1e-3)
R.check("edge layer  pi", edge, pi, source="closed form", tol=1e-3)
R.check("bulk fraction", 100*bulk/ai, 90.51, unit="%", source="geometric", tol=2e-3)
R.check("boundary fraction", 100*bnd/ai, 7.20, unit="%", source="geometric", tol=2e-3)
R.check("edge fraction", 100*edge/ai, 2.29, unit="%", source="geometric", tol=2e-3)
E_self = I(rhop**2) / (2*ai**2)
R.check("self-lensing energy E_self", E_self, 13.177, source="P04 / verify_P116", tol=1e-3,
        note="int(rho')^2 / 2(int rho)^2")
mu0, mu1 = I(rho), I(x*rho)
R.check("mass-hierarchy ratio MU = mu1/mu0", mu1/mu0, 0.79334, source="P03 / verify_P116", tol=2e-3)
# layer independence: deleting the edge term moves E_self out of the oscillation arena [4pi, E_self]
E2 = I((48*pi**3*x**2+6*pi**2*x)**2)/(2*I(16*pi**3*x**3+3*pi**2*x**2)**2)
R.check("layers non-collapsible", abs(E2-E_self) > 0.3, kind="fact",
        note="removing the edge moves E_self by %.2f (out of the arena)" % abs(E2-E_self))
R.emit()