"""Cross-code check of the caustic benchmark: the same cold Gaussian sphere at rest evolved with pbhgr's polar-areal
Einstein-Vlasov code (spherical_ev.py: polar slicing, zero shift, independent implementation) and sampled at equal
central proper time. Usage: python3 v1_caustic_polar.py [C=2M/R0] [K] [r_out/R0] [t_end/tC0]"""
import sys, time, numpy as np
sys.path.insert(0, '.')
from pbhgr.spherical_ev import Grid, Run, RunConfig, Particles, solve_metric
from pbhgr.cosmo_ev import CosmoGrid
from pbhgr.cosmo_id import sample_shells_quiet
C = float(sys.argv[1]) if len(sys.argv) > 1 else 0.02
K = int(sys.argv[2]) if len(sys.argv) > 2 else 4000
rout = float(sys.argv[3]) if len(sys.argv) > 3 else 4.0
t_end_fac = float(sys.argv[4]) if len(sys.argv) > 4 else 1.2
R0 = 1.0; x_out = 4.0 * R0
M = C * R0 / 2; rho0 = M / (np.pi**1.5 * R0**3)
Rg = np.linspace(0, x_out, 40001)
prof = dict(R=Rg, dN_dR=4 * np.pi * Rg**2 * rho0 * np.exp(-Rg**2 / R0**2), v=np.zeros_like(Rg))
cgrid = CosmoGrid(x_out=x_out, K=2000, stretch=4.0) # identical shells to the areal-CMC benchmark
x, P0, Lsq, N = sample_shells_quiet(prof, cgrid.x_edge, 4)
keep = x < rout * R0
x, N = x[keep], N[keep]
tC0 = np.pi / 2 / np.sqrt(8 * np.pi * rho0 / 3)
Rlist = (0.01, 0.03, 0.1, 0.3)
tlist = np.array([0.9, 0.98, 1.0, 1.01, 1.02, 1.05, 1.1, 1.2, 1.3]) * tC0
tlist = tlist[tlist <= t_end_fac * tC0 + 1e-9]
P = Particles(x.copy(), np.zeros_like(x), np.zeros_like(x), N.copy(), np.ones_like(x))
grid = Grid(r_out=rout * R0, K=K)
met0 = solve_metric(grid, P)
print(f"polar-areal: C={C}, K={K} (dr={grid.dr:.1e}), r_out={rout} R0, shells={len(x)} (sum N={N.sum():.5f}), M={met0.M:.5f}, e^mu(0)={np.exp(met0.mu_edge[0]):.5f}, max 2m/r={np.max(2*met0.m_edge[1:]/grid.r_edge[1:]):.4f}, t_C(0)={tC0:.3f}", flush=True)
run = Run(grid, P, RunConfig(t_end=tlist[0], cfl=0.4, diag_every=10**9, horizon_threshold=0.995))
def tau_core(): return float(np.median(P.tau[P.r < 0.05]))
rate = float(np.exp(met0.mu_edge[0])); t0 = time.time(); out = {}
for tt in tlist:
while True:
rem = (tt - tau_core()) / rate
if rem < 1e-7 * tC0: break
t_a, tau_a = run.t, tau_core()
run.cfg.t_end = run.t + rem
run.run(None)
if run.outcome == "BH_approach": break
if run.t > t_a: rate = max((tau_core() - tau_a) / (run.t - t_a), 0.05)
met = solve_metric(grid, P)
comp = 2 * met.m_edge[1:] / grid.r_edge[1:]
out[tt] = ([float(np.interp(Rf, grid.r_edge, met.m_edge)) for Rf in Rlist], float(comp.max()), float(np.exp(met.mu_edge[0])), met.n_bad, run.t)
print(f"tau/tC0={tt/tC0:5.2f} (t/tC0={run.t/tC0:.3f}): m(<0.01, 0.03, 0.1, 0.3)/M = " + " ".join(f"{v/M:7.4f}" for v in out[tt][0]) + f" max 2m/r={out[tt][1]:.3f} e^mu(0)={out[tt][2]:.4f} bad={met.n_bad} steps={run.step} ({time.time()-t0:.0f}s)", flush=True)
if run.outcome == "BH_approach":
print("polar-areal gauge breakdown (2m/r -> 1): stopping"); break