#!/usr/bin/env python3
"""Cold-limit convergence: nearly cold cloud (L_rms=0.05) at three (K, N) resolutions.
Reports ADM drift, max compactness and the time at which 2m/r first exceeds 0.5."""
import json, sys, time
import numpy as np
from pbhgr.spherical_ev import Grid, Run, RunConfig, solve_metric
from pbhgr.initial_data import shell_family

L_rms = float(sys.argv[1]) if len(sys.argv) > 1 else 0.05
nu, flatness = 0.15, 1.0
M = 1.0; R = 2 * M / nu; t_dyn = np.sqrt(R**3 / M)
rows = []
for K, n in ((300, 20000), (600, 40000), (1200, 80000)):
    rng = np.random.default_rng(11)
    g = Grid(r_out=6 * R, K=K)
    P = shell_family(M, R, n, rng, flatness=flatness, L_rms=L_rms)
    for _ in range(4): P.N *= M / solve_metric(g, P).M
    run = Run(g, P, RunConfig(t_end=3.0 * t_dyn, cfl=0.4, diag_every=25, horizon_threshold=0.995))
    t0 = time.time(); out = run.run()
    dg = run.diag; comp = np.array([d["max_2m_over_r"] for d in dg]); tt = np.array([d["t"] for d in dg]) / t_dyn
    Mt = np.array([d["M_total"] for d in dg])
    cross = tt[np.argmax(comp > 0.5)] if np.any(comp > 0.5) else None
    row = dict(K=K, n=n, outcome=out, max_comp=float(comp.max()), t_cross_0p5=cross,
               adm_drift_max=float(np.max(np.abs(Mt - Mt[0]))), adm_drift_at_end=float(Mt[-1] - Mt[0]),
               t_end=float(tt[-1]), r_min=float(np.min([d["r_min_particle"] for d in dg])), wall=time.time() - t0)
    rows.append(row); print(json.dumps(row), flush=True)
json.dump(rows, open(f"runs/verify_cold_L{L_rms}.json", "w"), indent=1)