#!/usr/bin/env python3
"""
PROVENANCE: NUMERICAL

PST Computation 139 — the d^-6 Casimir exponent derived from the projection kernel
==================================================================================
Advances the Casimir prediction (sec:casimir) from a dimensional-analysis /
parity ARGUMENT to a kernel-level DERIVATION, and quantifies its
(non-)distinctiveness.

The PST vacuum modifies the plate-mode spectrum by a form factor g(k d_0) set
by the substrate projection kernel (Comps 68, 73). The scalar Casimir pressure
between plates at separation d is

    P(d) = -C int_0^inf dk k^3 g(k d_0) / (e^{2 k d} - 1),   g(0)=1,

with g even and smooth (g'(0)=0: no preferred orientation / parity). Expanding
g(x) = 1 + a2 x^2 + a4 x^4 + ..., the leading correction is

    dP(d) = -C a2 d_0^2 int_0^inf dk k^5/(e^{2kd}-1)
          = -C a2 d_0^2 * Gamma(6) zeta(6) / (2 d)^6    ~   d_0^2 d^-6.

So the correction pressure scales as d^-6 for ANY smooth even kernel: the
EXPONENT is kernel-independent (hence NOT PST-distinctive; any single-length
smooth even cutoff gives d^-6), while the COEFFICIENT a2 (proportional to the
kernel's second moment) is kernel-specific and carries xi = 90/pi^2 for the
Gaussian.

Checks:
 (1) numerically integrate dP(d) for several smooth even kernels; fit the
     log-log slope and confirm exponent = -6 to <1% for all of them;
 (2) confirm the coefficient (xi-proxy |dP| d^6 / d0^2) is d-independent but
     kernel-specific (differs across kernels), matching Comp 68;
 (3) check the analytic leading term Gamma(6) zeta(6)/2^6 * a2 against numerics;
 (4) show a NON-smooth even kernel (exp(-|x|), a cusp at 0 with a |x| term)
     does NOT give -6, confirming smoothness+parity are the load-bearing inputs.
"""
import numpy as np
from scipy import integrate, special

np.seterr(over="ignore")   # e^{2kd} overflows at large k; integrand -> 0 there anyway

d0 = 1e-3

smooth_even = {
    "gaussian":       (lambda x: np.exp(-0.5 * x**2),      -0.5),   # a2 = g''(0)/2
    "sech^2":         (lambda x: 1.0 / np.cosh(x)**2,      -1.0),
    "lorentzian":     (lambda x: 1.0 / (1.0 + x**2),       -1.0),
    "sq-lorentzian":  (lambda x: 1.0 / (1.0 + x**2)**2,    -2.0),
}
nonsmooth = {"exp(-|x|)": lambda x: np.exp(-np.abs(x))}   # even but cusp -> |x| term


def dP(d, g):
    f = lambda k: k**3 * (g(k * d0) - 1.0) / np.expm1(2 * k * d)
    val, _ = integrate.quad(f, 0, np.inf, limit=400)
    return val


ds = np.array([0.05, 0.07, 0.1, 0.15, 0.2, 0.3])
analytic_scale = special.gamma(6) * special.zeta(6) / 2**6   # Gamma(6)zeta(6)/2^6

print(f"{'kernel':16} {'fit exponent':>12} {'coeff/d0^2':>14} {'analytic a2*C6':>15} {'ratio':>7}")
slopes = []
for name, (g, a2) in smooth_even.items():
    vals = np.array([dP(d, g) for d in ds])
    slope, _ = np.polyfit(np.log(ds), np.log(np.abs(vals)), 1)
    slopes.append(slope)
    coeff = np.mean(np.abs(vals) * ds**6 / d0**2)          # xi-proxy
    analytic = abs(a2) * analytic_scale                     # |dP| ~ |a2| C6 d0^2 d^-6
    print(f"{name:16} {slope:12.4f} {coeff:14.4f} {analytic:15.4f} {coeff/analytic:7.3f}")

# non-smooth kernel: different exponent
for name, g in nonsmooth.items():
    vals = np.array([dP(d, g) for d in ds])
    slope, _ = np.polyfit(np.log(ds), np.log(np.abs(vals)), 1)
    print(f"{name:16} {slope:12.4f}   (non-smooth: cusp gives a |x| term, exponent != -6)")

exp_ok = all(abs(s + 6.0) < 0.06 for s in slopes)           # all smooth even -> -6
coeffs = []
for name, (g, a2) in smooth_even.items():
    v = np.array([dP(d, g) for d in ds])
    coeffs.append(np.mean(np.abs(v) * ds**6 / d0**2))
coeff_varies = (max(coeffs) / min(coeffs)) > 1.5            # coefficient is kernel-specific

print(f"\nexponent = -6 for all smooth even kernels: {exp_ok}")
print(f"coefficient kernel-specific (varies by {max(coeffs)/min(coeffs):.1f}x): {coeff_varies}")
print("\nRESULT:", "PASS  (d^-6 exponent derived + kernel-independent; coefficient kernel-specific)"
      if (exp_ok and coeff_varies) else "FAIL")