#!/usr/bin/env python3
"""
SCOPE NOTE, 2026-10-13 (comps 201 and 216).  Every number below is evaluated at d_0 =
2.3 nm, the Casimir bound of its day, and is kept as found: (k d_0)^2 = 2.1e-24 and
r_break = 38 micrometres for a solar mass, the floor M = 2.9e18 kg.  comp 201 derived
d_0 = sqrt(2) hbar c / m_h = 2.2 am, and comp 216 re-evaluated this computation on that
length: (k d_0)^2 = 1.9e-42, r_break = 37 pm, the floor 2.8e9 kg.  The conclusion, that
the derivative expansion holds at every astrophysical horizon, is unchanged and stronger.

PROVENANCE: PROOF

Computation 154 -- The derivative expansion does not fail at an astrophysical horizon
=====================================================================================
Sec 8.1 declines to claim that the projected field equations (eq:einstein-derived)
cover the strong-field regime, and Sec 19.x ("The black-hole interior") records that
this is the one step standing between PST and the Wald Noether-charge entropy.  This
computation asks what the declining actually costs, by evaluating the book's OWN
validity criterion -- the higher-curvature terms of eq:projected-action are
"suppressed by (k d_0)^2 and negligible at all scales >> d_0" -- on the Schwarzschild
geometry.  It is the book's criterion, not a new one.

The answer is that "strong field" and "outside the expansion" are not the same region,
and at a horizon they are twenty-four orders of magnitude apart.

CHECKS
------
  (1) SCHWARZSCHILD IS RICCI-FLAT, exactly and symbolically.  R_mn = 0, R = 0, and the
      only non-vanishing curvature invariant is the Kretschmann scalar
      K = 12 r_s^2 / r^6.  Computed here from the metric via Christoffels and the
      Riemann tensor, not quoted.

  (2) THE EXPANSION PARAMETER AT THE HORIZON.  Take the local curvature scale to be
      k(r) = K(r)^(1/4), so k(r_s) = 12^(1/4)/r_s.  The book's suppression factor is
      then (k d_0)^2 = sqrt(12) (d_0/r_s)^2, which at r = r_s is 2.1e-24 for one solar
      mass and smaller for every heavier hole.  The horizon is not a strong-curvature
      region in the only sense the expansion cares about; it is strong *potential*,
      which the covariant expansion does not see.

  (3) WHERE THE EXPANSION DOES FAIL.  Setting k d_0 = 1 gives
      r_break = (12^(1/4) d_0)^(2/3) r_s^(1/3), i.e. r_break/r_s = (12^(1/4) d_0/r_s)^(2/3).
      For one solar mass that is 38 micrometres, 1.3e-8 of the horizon radius, enclosing
      2.1e-24 of the interior volume.  The region outside the projected description is
      the neighbourhood of the singularity, not the horizon -- which is exactly what
      Sec 8.1 already says is open ("no result on the central singularity").

  (4) THE MASS FLOOR, and its robustness to which PST length controls the expansion.
      k d_0 = 1 AT the horizon requires r_s = 12^(1/4) d_0 = 4.3 nm, i.e. M = 2.9e18 kg
      = 3.3e-12 solar masses.  Above that floor the horizon is inside the expansion's
      domain.  If the controlling length were the grain l_* = hbar c / M_* rather than
      the coherence length d_0, the floor drops to 1.9e8 kg: the conclusion does not
      depend on which of PST's two lengths is the cutoff, only on both being far below
      any astrophysical r_s.

  (5) THE NAMED CORRECTIONS VANISH IDENTICALLY, so the truncation is not merely small
      but exact here.  eq:projected-action names O(R_mn R^mn, R^2).  Both their
      contributions to the field equations and their contributions to Wald's
      Noether-charge integrand are proportional to R or to R_mn, which check (1) shows
      vanish.  Verified symbolically on the Schwarzschild solution.  Consequences:
      Schwarzschild solves the truncated equations exactly, and the Wald entropy of the
      truncated Lagrangian is exactly A/4G with no (d_0/r_s)^2 correction at all.
      PST's derived cosmological constant spoils exact Ricci-flatness -- the exterior is
      Schwarzschild-de Sitter, R = 4 Lambda -- so the check also carries the resulting
      correction: Lambda r_s^2 = 1e-45 for a solar mass, and the relative shift in the
      Wald entropy is 1e-119.

  (6) THE ONE TERM NOT NAMED is the Riemann-squared invariant, which a Seeley-DeWitt
      a_2 coefficient generally carries.  In four dimensions it is the Gauss-Bonnet
      combination up to the two terms already killed by (5), and its Wald contribution
      is computed here to be 64 pi^2 lambda -- INDEPENDENT OF THE MASS.  It is an
      additive constant, not an area term, so the area law survives it; against
      S_BH/k_B = 1.05e77 for one solar mass an O(1) shift is a relative 1e-77.

WHAT THIS DOES NOT ESTABLISH.  Three things, stated so the result is not oversold.
It does not derive the global causal structure: that r = r_s is an event rather than
merely an apparent horizon rests on the interior being censored, exactly as it does in
general relativity, and PST inherits that question unchanged rather than answering it.
It does not touch the substrate-level identification of Sec 8.1 -- that a horizon IS
the projection of the null-threshold set is still an identification, and nothing here
makes it a theorem.  And it says nothing about the statistical entropy; Wald's
construction contains no microstate count, which is precisely why it is available.

A "pass": (1) R_mn and R symbolically zero and K = 12 r_s^2/r^6; (2) the horizon
suppression factor below 1e-20 for every stellar-mass and heavier hole; (3) r_break
inside 1e-7 of r_s; (5) both named corrections symbolically zero on the solution;
(6) the Gauss-Bonnet Wald contribution independent of r_s.
"""

from __future__ import annotations
import math
import sympy as sp

# ---------------------------------------------------------------- constants
G_N   = 6.67430e-11          # CODATA 2018, m^3 kg^-1 s^-2
C_L   = 2.99792458e8         # m/s, exact
HBAR  = 1.054571817e-34      # J s
K_B   = 1.380649e-23         # J/K, exact
M_SUN = 1.98892e30           # kg
D0    = 2.3e-9               # m   -- coherence length, diagnostic bound eq:d0bound
L_STAR = 1.5e-19             # m   -- grain hbar c / M_*, M_* = 4 pi m_h sqrt(2/3)
LAMBDA = 1.1e-52             # m^-2 -- observed cosmological constant, Table of symbols
L_P2   = HBAR * G_N / C_L**3 # m^2  -- Planck area, the G of the geometric units below

def r_s(M):
    return 2.0 * G_N * M / C_L**2


# ------------------------------------------------- (1) symbolic Schwarzschild
def schwarzschild_curvature():
    """Christoffels -> Riemann -> Ricci, R, Kretschmann.  Nothing quoted."""
    t, r, th, ph, rs = sp.symbols('t r theta phi r_s', positive=True)
    x = [t, r, th, ph]
    f = 1 - rs / r
    g = sp.diag(-f, 1 / f, r**2, r**2 * sp.sin(th)**2)
    ginv = g.inv()
    n = 4

    Gam = [[[sp.simplify(sum(ginv[a, d] * (sp.diff(g[d, b], x[c])
                                           + sp.diff(g[d, c], x[b])
                                           - sp.diff(g[b, c], x[d]))
                             for d in range(n)) / 2)
             for c in range(n)] for b in range(n)] for a in range(n)]

    def riem_mixed(a, b, c, d):                      # R^a_{bcd}
        e = sp.diff(Gam[a][b][d], x[c]) - sp.diff(Gam[a][b][c], x[d])
        e += sum(Gam[a][c][s] * Gam[s][b][d] - Gam[a][d][s] * Gam[s][b][c]
                 for s in range(n))
        return sp.simplify(e)

    Rm = [[[[riem_mixed(a, b, c, d) for d in range(n)] for c in range(n)]
           for b in range(n)] for a in range(n)]
    Ric = sp.Matrix(n, n, lambda b, d: sp.simplify(sum(Rm[a][b][a][d] for a in range(n))))
    Rsc = sp.simplify(sum(ginv[b, d] * Ric[b, d] for b in range(n) for d in range(n)))

    # the metric is diagonal, so lowering and raising are factor-wise
    Rl = [[[[sp.simplify(g[a, a] * Rm[a][b][c][d]) for d in range(n)] for c in range(n)]
           for b in range(n)] for a in range(n)]
    K = 0
    for a in range(n):
        for b in range(n):
            for c in range(n):
                for d in range(n):
                    if Rl[a][b][c][d] != 0:
                        K += (ginv[a, a] * ginv[b, b] * ginv[c, c] * ginv[d, d]
                              * Rl[a][b][c][d]**2)
    K = sp.simplify(K)
    return dict(sym=(t, r, th, ph, rs), g=g, ginv=ginv, Ric=Ric, R=Rsc, Rl=Rl, K=K)


def main():
    print("=" * 100)
    print("  Computation 154 -- the derivative expansion does not fail at an astrophysical horizon")
    print("=" * 100)
    print()

    # ---- (1) -------------------------------------------------------------
    print("  (1) SCHWARZSCHILD CURVATURE, computed symbolically from the metric")
    S = schwarzschild_curvature()
    t, r, th, ph, rs = S['sym']
    ricci_zero  = all(sp.simplify(e) == 0 for e in S['Ric'])
    scalar_zero = sp.simplify(S['R']) == 0
    K_expected  = sp.simplify(S['K'] - 12 * rs**2 / r**6) == 0
    print(f"      R_mn == 0 (all sixteen components) : {ricci_zero}")
    print(f"      R    == 0                          : {scalar_zero}")
    print(f"      K    == 12 r_s^2 / r^6             : {K_expected}      [K = {sp.simplify(S['K'])}]")
    check1 = ricci_zero and scalar_zero and K_expected
    print()

    # ---- (2) -------------------------------------------------------------
    print("  (2) THE BOOK'S SUPPRESSION FACTOR (k d_0)^2 EVALUATED AT THE HORIZON, k = K^(1/4)")
    print(f"      {'object':<24}{'M / kg':>12}{'r_s / m':>14}{'(k d_0)^2':>14}")
    ladder = [("Moon-mass PBH",      7.35e22),
              ("Earth-mass PBH",     5.97e24),
              ("1 solar mass",       M_SUN),
              ("10 solar masses",    10 * M_SUN),
              ("Cyg X-1, 21 M_sun",  21 * M_SUN),
              ("Sgr A*, 4.3e6 M_sun", 4.3e6 * M_SUN),
              ("M87*, 6.5e9 M_sun",  6.5e9 * M_SUN)]
    horizon_factors = {}
    for name, M in ladder:
        R = r_s(M)
        kd2 = math.sqrt(12.0) * (D0 / R)**2          # (K^(1/4) d_0)^2 at r = r_s
        horizon_factors[name] = kd2
        print(f"      {name:<24}{M:12.3e}{R:14.4e}{kd2:14.3e}")
    check2 = all(v < 1e-20 for k, v in horizon_factors.items() if "PBH" not in k)
    print(f"      -> below 1e-20 for every stellar-mass and heavier hole: {check2}")
    print()

    # ---- (3) -------------------------------------------------------------
    print("  (3) WHERE THE EXPANSION ACTUALLY FAILS: r_break = (12^(1/4) d_0)^(2/3) r_s^(1/3)")
    print(f"      {'object':<24}{'r_s / m':>14}{'r_break / m':>14}{'r_break/r_s':>14}{'vol. frac.':>13}")
    worst = 0.0
    for name, M in [("1 solar mass", M_SUN), ("10 solar masses", 10 * M_SUN),
                    ("Sgr A*", 4.3e6 * M_SUN), ("M87*", 6.5e9 * M_SUN)]:
        R = r_s(M)
        rb = (12.0**0.25 * D0)**(2.0 / 3.0) * R**(1.0 / 3.0)
        worst = max(worst, rb / R)
        print(f"      {name:<24}{R:14.4e}{rb:14.4e}{rb / R:14.3e}{(rb / R)**3:13.2e}")
    check3 = worst < 1e-7
    print(f"      -> the excluded region is the singularity's neighbourhood, not the horizon: {check3}")
    print()

    # ---- (4) -------------------------------------------------------------
    print("  (4) MASS FLOOR, k d_0 = 1 AT the horizon, and robustness to the choice of length")
    for L, label in [(D0, "coherence length d_0 = 2.3 nm"),
                     (L_STAR, "grain  l_* = 1.5e-19 m     ")]:
        R = 12.0**0.25 * L
        M = R * C_L**2 / (2 * G_N)
        T_H = HBAR * C_L**3 / (8 * math.pi * G_N * M * K_B)
        print(f"      {label}  ->  r_s = {R:.3e} m,  M = {M:.3e} kg = {M / M_SUN:.2e} M_sun,"
              f"  T_H = {T_H:.2e} K")
    floor_d0 = 12.0**0.25 * D0 * C_L**2 / (2 * G_N)
    check4 = floor_d0 < 1e19 and floor_d0 / M_SUN < 1e-10
    print(f"      -> both floors lie far below any astrophysical mass; the conclusion does")
    print(f"         not depend on which PST length is the cutoff: {check4}")
    print()

    # ---- (5) -------------------------------------------------------------
    print("  (5) THE NAMED CORRECTIONS O(R_mn R^mn, R^2) VANISH IDENTICALLY ON THE SOLUTION")
    Ric, Rsc, ginv, g = S['Ric'], S['R'], S['ginv'], S['g']
    # field-equation contributions, standard variations, evaluated on Schwarzschild
    #   from  R^2        :  2 grad_m grad_n R - 2 g_mn box R - (1/2) g_mn R^2 + 2 R R_mn
    #   from  R_ab R^ab  :  every term carries a factor R_mn or R (box R_mn, grad grad R, ...)
    ric2 = sp.simplify(sum(ginv[a, a] * ginv[b, b] * Ric[a, b]**2 for a in range(4) for b in range(4)))
    print(f"      R           = {sp.simplify(Rsc)}   -> every term of the R^2 variation carries R")
    print(f"      R_ab R^ab   = {ric2}   -> every term of the Ricci-squared variation carries R_ab")
    print(f"      sum |R_ab|  = {sp.simplify(sum(sp.Abs(e) for e in Ric))}")
    print(f"      Wald integrand dL/dR_abcd: the R^2 piece is 2R g..g.. and the Ricci-squared")
    print(f"      piece is built from R^bd g^ac, so both are zero by the two lines above.")
    # A critic's objection: PST derives a cosmological constant, so the exterior is
    # Schwarzschild-de Sitter and is NOT Ricci-flat -- R = 4 Lambda, R_mn = Lambda g_mn.
    # Both corrections then survive, at these magnitudes.
    lam_rs2 = LAMBDA * r_s(M_SUN)**2
    rel_wald = 128 * math.pi * L_P2 * LAMBDA          # dS/S from an R^2 term, alpha ~ 1
    print(f"      with the derived Lambda = {LAMBDA:.1e} m^-2 the exterior is Schwarzschild-de Sitter,")
    print(f"      not Ricci-flat: Lambda r_s^2 = {lam_rs2:.1e} for one solar mass, and the resulting")
    print(f"      relative Wald correction 128 pi alpha G Lambda = {rel_wald:.1e} with alpha ~ 1.")
    check5 = (sp.simplify(Rsc) == 0) and (ric2 == 0) and rel_wald < 1e-100
    print(f"      -> Schwarzschild solves the TRUNCATED equations exactly, and S_Wald = A/4G")
    print(f"         with no (d_0/r_s)^2 correction: {check5}")
    print()

    # ---- (6) -------------------------------------------------------------
    print("  (6) THE UNNAMED Riemann^2 TERM: Gauss-Bonnet, and its Wald contribution is topological")
    # On a Ricci-flat background dL_GB/dR_abcd = 2 R^abcd.  With the horizon binormal
    # eps_ab (eps_tr = 1, since |g_tt g_rr| = 1 for Schwarzschild; eps_ab eps^ab = -2),
    #     R^abcd eps_ab eps_cd = 4 R^{trtr} = 4 (g^tt g^rr)^2 R_trtr = 4 R_trtr .
    Rl = S['Rl']
    # convention self-test: the same formula on L = R/(16 pi G) must return A/4G.
    #   dL/dR_abcd = (1/16 pi G) * (1/2)(g^ac g^bd - g^ad g^bc);  contracted with the
    #   binormal this is -2/(16 pi G), and S = -2 pi * (-2/(16 pi G)) * A = A/(4G).
    eh_contract = sp.Rational(1, 2) * (-2 - 2)                # eps.eps = -2 on each pairing
    eh_entropy_coeff = sp.simplify(-2 * sp.pi * eh_contract / (16 * sp.pi))   # times A/G
    print(f"      convention self-test: Wald applied to L = R/(16 pi G) gives S = "
          f"{eh_entropy_coeff} * A/G   (must be 1/4)")
    R_trtr = sp.simplify(Rl[0][1][0][1])
    norm = sp.simplify(2 * ginv[0, 0] * ginv[1, 1])           # eps_ab eps^ab, must be -2
    contraction = sp.simplify(4 * (ginv[0, 0] * ginv[1, 1])**2 * R_trtr)
    at_horizon = sp.simplify(contraction.subs(r, rs))
    # S_GB = -2 pi lambda * Int (2 R^abcd eps_ab eps_cd) sqrt(h) d^2x , area 4 pi r_s^2
    S_GB = sp.simplify(-2 * sp.pi * 2 * at_horizon * 4 * sp.pi * rs**2)
    print(f"      R_trtr                       = {R_trtr}")
    print(f"      binormal normalisation eps.eps = {norm}   (must be -2)")
    print(f"      R^abcd eps_ab eps_cd at r=r_s  = {at_horizon}")
    print(f"      S_GB / lambda                  = {S_GB}      <- no r_s, hence no area")
    S_BH = 4 * math.pi * G_N * M_SUN**2 / (HBAR * C_L)
    print(f"      S_BH/k_B for one solar mass    = {S_BH:.3e}; an O(1) additive shift is {1 / S_BH:.1e} of it")
    check6 = (sp.simplify(norm + 2) == 0) and (sp.diff(S_GB, rs) == 0) \
        and (eh_entropy_coeff == sp.Rational(1, 4))
    print(f"      -> additive constant, not an area term; the area law is untouched: {check6}")
    print()

    # ---- assessment -------------------------------------------------------
    print("=" * 100)
    print("  ASSESSMENT")
    print("=" * 100)
    print(f"  (1) Schwarzschild Ricci-flat, K = 12 r_s^2/r^6, symbolically         : {check1}")
    print(f"  (2) (k d_0)^2 < 1e-20 at the horizon for every astrophysical hole    : {check2}")
    print(f"  (3) breakdown radius inside 1e-7 of r_s (the singularity, not the horizon) : {check3}")
    print(f"  (4) mass floor 2.9e18 kg (d_0) / 1.9e8 kg (l_*), both far below any hole   : {check4}")
    print(f"  (5) the named higher-curvature corrections vanish identically on the solution : {check5}")
    print(f"  (6) the Gauss-Bonnet Wald contribution is mass-independent           : {check6}")
    print()
    ok = check1 and check2 and check3 and check4 and check5 and check6
    print("  RESULT: " + ("NARROWS Sec 8.1 -- the projected field equations are valid AT the horizon by "
                          "twenty-four orders of margin; what they do not cover is the singularity"
                          if ok else "MISMATCH -- revise"))


if __name__ == "__main__":
    main()
