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

PST Computation 137 — Equation of state w = -1 from steady-state non-dilution
=============================================================================
Witness computation for the cosmological equation-of-state claim
(sec:cosmology-eos): the instantiation stress-energy has w = -1.

Two independent exact routes, both verified here:

(A) Continuity route. In an FRW universe the covariant conservation of the
    stress-energy gives the continuity equation
        rho_dot + 3 H (rho + p) = 0,   with solution   rho(a) = rho_0 a^{-3(1+w)}.
    The PST instantiation density is non-diluting (rho_dot = 0), because the
    substrate source has no spatial extent to stretch. rho constant forces the
    exponent -3(1+w) = 0, i.e. w = -1. Verified symbolically and numerically.

(B) Conformal / vacuum-form route. A stress tensor proportional to the metric,
    T_{mu nu} = rho_Lambda g_{mu nu} (the form the substrate contributes,
    isotropic scalar projection), has mixed form T^mu_nu = rho_Lambda delta^mu_nu.
    Matching to a perfect fluid T^mu_nu = diag(rho, -p, -p, -p) gives p = -rho,
    hence w = -1. Verified by explicit index manipulation.

Exact symbolic + linear algebra, hence PROVENANCE PROOF. Note the load-bearing
INPUT is non-dilution (rho_dot = 0); this computation shows w = -1 FOLLOWS from
it, it does not derive non-dilution (that is the steady-state premise of
sec:cosmology-eos).
"""
import sympy as sp

# ---------- (A) continuity route ----------
a, rho0, H, t = sp.symbols("a rho0 H t", positive=True)
w = sp.symbols("w", real=True)          # w may be negative (indeed w = -1)
rho = rho0 * a ** (-3 * (1 + w))                      # general-w solution of continuity
# non-dilution: d rho / d a = 0 for all a  <=>  exponent zero
drho_da = sp.simplify(sp.diff(rho, a))
# exponent of a in rho:
exponent = -3 * (1 + w)
w_from_nondilution = sp.solve(sp.Eq(exponent, 0), w)   # exponent = 0
routeA = (w_from_nondilution == [-1])

# numeric sanity: density ratio rho(2)/rho(1) for several w (1 = non-diluting)
ratios = {ww: float((2.0) ** (-3 * (1 + ww))) for ww in (-1.0, -1.0 / 3, 0.0, 1.0 / 3)}

# ---------- (B) conformal / vacuum-form route ----------
# metric signature (+,-,-,-); T_{mu nu} = rhoL g_{mu nu}
rhoL = sp.symbols("rho_L", positive=True)
g = sp.diag(1, -1, -1, -1)                             # Minkowski local frame
T_lower = rhoL * g                                     # T_{mu nu}
T_mixed = g.inv() * T_lower                            # T^mu_nu = g^{mu a} T_{a nu}
# perfect fluid comparison: T^mu_nu = diag(rho, -p, -p, -p)
rho_val = T_mixed[0, 0]
p_val = -T_mixed[1, 1]
w_B = sp.simplify(p_val / rho_val)
routeB = (w_B == -1)

print("(A) continuity + non-dilution:")
print(f"    rho(a) = rho0 * a^(-3(1+w));  d rho/d a = 0  =>  w = {w_from_nondilution}")
print(f"    density ratio rho(2)/rho(1) by w: {{" +
      ", ".join(f'{k:+.3f}:{v:.4f}' for k, v in ratios.items()) + "}")
print(f"    (only w=-1 gives ratio 1.0000, i.e. non-diluting)   route A -> w=-1: {routeA}")
print("\n(B) T_{mu nu} = rho_L g_{mu nu}  =>  T^mu_nu = diag" +
      f"({T_mixed[0,0]}, {T_mixed[1,1]}, {T_mixed[2,2]}, {T_mixed[3,3]})")
print(f"    rho = {rho_val}, p = {p_val}, w = p/rho = {w_B}   route B -> w=-1: {routeB}")

ok = routeA and routeB and abs(ratios[-1.0] - 1.0) < 1e-12
print("\nRESULT:", "PASS  (w = -1 by both routes)" if ok else "FAIL")