#!/usr/bin/env python3
"""Step 6.  Collect the key numbers into data/results.json and print the
tables quoted in report.md (also written to data/report_tables.md)."""
import os

import numpy as np
import pandas as pd

from common import DATA, KMS_TO_PCMYR, load_json, save_json


def fmt(x, n=3):
    if x is None or (isinstance(x, float) and not np.isfinite(x)):
        return "n/a"
    return f"{x:.{n}g}"


def rule(dv):
    return "any speed (upper bound)" if np.isinf(dv) else f"dv < {dv:g} km/s"


def main():
    s1 = load_json("s1_density_kinematics.json")
    s2 = load_json("s2_census.json")
    s3 = load_json("s3_exposure.json")
    s4 = load_json("s4_settlement.json")
    lines = []
    P = lines.append

    P("## Table 1. Density, speed, rate\n")
    P("| quantity | value |\n|---|---|")
    P(f"| GCNS rows | {s1['gcns_rows']} |")
    P(f"| GCNS systems <50 pc (after merging) | {s1['n_systems_gcns_50pc_count']} |")
    P(f"| GCNS system density <50 pc | {s1['n_systems_gcns_50pc']:.4f} +- {s1['n_systems_gcns_50pc_err']:.4f} pc^-3 |")
    t = s1["tenpc"]
    P(f"| 10 pc sample stellar systems | {t['n_systems_stellar']} ({t['dens_systems_stellar']:.4f} +- {t['dens_systems_stellar_err']:.4f} pc^-3) |")
    P(f"| ... of which present in GCNS | {t['stellar_systems_with_gcns_member']} ({100 * t['gcns_system_completeness_10pc']:.1f}%) |")
    P(f"| GCNS sources within 10 pc / of which brown dwarfs | {t['gcns_sources_10pc']} / {t['gcns_sources_10pc_brown_dwarfs']} |")
    P(f"| GCNS systems within 10 pc / with a brown-dwarf brightest member | {t['gcns_systems_10pc']} / {t['gcns_systems_10pc_bd_primary']} |")
    P(f"| adopted system density | {s1['n_sys_adopted']:.4f} +- {s1['n_sys_adopted_err']:.4f} pc^-3 |")
    P(f"| mean speed relative to Sun, (4/pi)<vT> | {s1['mean_speed_adopted_kms']:.2f} +- {s1['mean_speed_adopted_err']:.2f} km/s |")
    P(f"| check: <|v|> of {s1['check_rv_subsample_N']} clean-RV stars | {s1['check_v3d_mean_kms']:.2f} +- {s1['check_v3d_mean_err']:.2f} (3D) vs {s1['check_4overpi_vT_kms']:.2f} +- {s1['check_4overpi_vT_err']:.2f} (4/pi <vT>) |")
    P(f"| RV-availability bias, 4<=MG<12 | with RV {s1['bias_MG4to12_v_rv']:.1f} (N={s1['bias_MG4to12_N_rv']}), without {s1['bias_MG4to12_v_norv']:.1f} +- {s1['bias_MG4to12_v_norv_err']:.1f} km/s (N={s1['bias_MG4to12_N_norv']}) |")
    P(f"| Gamma(1 pc), raw GCNS flux sum | {s1['gamma1pc_gcns_raw']:.2f} +- {s1['gamma1pc_gcns_raw_staterr']:.2f} /Myr |")
    P(f"| Gamma(1 pc), adopted | {s1['gamma1pc_adopted']:.2f} +- {s1['gamma1pc_adopted_err']:.2f} systems/Myr (half-range of the variants below; internal error {s1['gamma1pc_internal_err']:.2f}) |")
    P(f"| 10 pc sample alone: density x own mean 3D speed ({t['tenpc_own_speed_N']} systems, {t['tenpc_own_speed_mean_kms']:.1f} +- {t['tenpc_own_speed_err_kms']:.1f} km/s) | {t['gamma1pc_tenpc_own_speed']:.2f} +- {t['gamma1pc_tenpc_own_speed_err']:.2f} /Myr |")
    P(f"| model-free lower bound ({s2['lower_bound_window']}; N = {s2['lower_bound_N']}) | >= {s2['lower_bound_gamma1pc']:.2f} +- {s2['lower_bound_gamma1pc_err']:.2f} /Myr |")
    P(f"| velocity sample for simulations | {s1['velocity_sample_N']} stars (N_eff {s1['velocity_sample_Neff']:.0f}); {s1['velocity_sample_N_speed_gt_200kms']} faster than 200 km/s; sigma(U,V,W) = {s1['sigma_U_kms']:.1f}, {s1['sigma_V_kms']:.1f}, {s1['sigma_W_kms']:.1f} km/s |")

    P("\nRate variants (stellar systems), Gamma(1 pc) per Myr:\n")
    rv = pd.read_csv(os.path.join(DATA, "rate_variants.csv"))
    P("| variant | Gamma(1 pc) |\n|---|---|")
    for _, r in rv.iterrows():
        P(f"| {r.variant} | {r.gamma1pc:.2f} |")
    P(f"| (other definitions) every stellar object of the 10 pc sample counted separately | {s3['gamma1pc_all_stellar_objects_10pc']:.2f} |")
    P(f"| (other definitions) all 10 pc systems including brown-dwarf primaries | {s3['gamma1pc_all_systems_incl_BD_10pc']:.2f} |")

    P("\nModel-free lower bound, raw clean-RV counts in all of Gaia DR3 (no completeness model):\n")
    lb = pd.read_csv(os.path.join(DATA, "rate_lower_bound.csv"))
    P("| |t| < [Myr] | d < [pc] | N | N / (2 T d^2) [per Myr at 1 pc] |\n|---|---|---|---|")
    for _, r in lb.iterrows():
        P(f"| {r.T_Myr:g} | {r.d_pc:g} | {int(r.N_clean)} | {r.gamma1pc_lower:.2f} +- {r.gamma1pc_lower_err:.2f} |")

    P("\n## Table 2. Direct census (clean DR3 RV, now <100 pc, |t|<1 Myr)\n")
    cen = pd.read_csv(os.path.join(DATA, "census_obs_vs_pred.csv"))
    m = cen[(cen.T_Myr == 1.0) & (cen.ruwe_cut == 0) & np.isinf(cen.vmax)]
    P(f"parent systems <100 pc: {s2['parent_systems_100pc']}; with clean RV: {s2['parent_with_clean_rv']}\n")
    P("| d [pc] | N_obs | N_pred | obs/pred (68% Poisson) | RV completeness |\n|---|---|---|---|---|")
    for _, r in m.iterrows():
        P(f"| {r.d_pc:g} | {int(r.N_obs)} | {r.N_pred:.1f} | {r.ratio:.2f} ({r.ratio_lo:.2f}-{r.ratio_hi:.2f}) | {r.rv_completeness:.2f} |")
    P(f"\nslope alpha of N(<d) ~ d^alpha: {s2['slope_dph_alpha']:.2f} +- {s2['slope_dph_alpha_err']:.2f}")

    # derived consistency numbers quoted in the report
    from scipy.stats import poisson
    r05 = m[m.d_pc == 0.5].iloc[0]; r1 = m[m.d_pc == 1.0].iloc[0]; r5 = m[m.d_pc == 5.0].iloc[0]
    BJ18 = 19.7
    scale = BJ18 / s1["gamma1pc_adopted"]
    chk = dict(
        P_zero_within_0p5pc=float(poisson.pmf(0, r05.N_pred)),
        P_le_obs_within_1pc=float(poisson.cdf(int(r1.N_obs), r1.N_pred)),
        P_le_obs_within_5pc=float(poisson.cdf(int(r5.N_obs), r5.N_pred)),
        N_pred_5pc_if_rate_were_BJ18=float(r5.N_pred * scale),
        N_pred_1pc_if_rate_were_BJ18=float(r1.N_pred * scale),
        P_le_obs_within_1pc_if_BJ18=float(poisson.cdf(int(r1.N_obs), r1.N_pred * scale)),
        gamma1pc_from_10pc_stellar_systems=float(np.pi * t["dens_systems_stellar"]
                                                 * s1["mean_speed_adopted_kms"] * KMS_TO_PCMYR),
        gamma1pc_from_10pc_stellar_systems_err=float(np.pi * t["dens_systems_stellar_err"]
                                                     * s1["mean_speed_adopted_kms"] * KMS_TO_PCMYR),
        ratio_adopted_over_BJ18=float(1 / scale),
        lower_bound_over_BJ18=float(s2["lower_bound_gamma1pc"] / BJ18),
    )
    P("\nDerived checks:")
    for k, v in chk.items():
        P(f"- {k}: {fmt(v)}")
    for T in (0.5, 1.5):
        mm = cen[(cen.T_Myr == T) & (cen.ruwe_cut == 0) & np.isinf(cen.vmax)]
        P(f"\n|t| < {T:g} Myr: " + "; ".join(f"d={r.d_pc:g}: {int(r.N_obs)}/{r.N_pred:.1f}" for _, r in mm.iterrows()))
    mm = cen[(cen.T_Myr == 1.0) & (cen.ruwe_cut == 1.4) & np.isinf(cen.vmax)]
    P("\nRUWE < 1.4 only, |t| < 1 Myr: " + "; ".join(f"d={r.d_pc:g}: {int(r.N_obs)}/{r.N_pred:.1f}" for _, r in mm.iterrows()))
    mm = cen[(cen.T_Myr == 1.0) & (cen.ruwe_cut == 0) & (cen.vmax == 150)]
    P("\nv < 150 km/s only, |t| < 1 Myr: " + "; ".join(f"d={r.d_pc:g}: {int(r.N_obs)}/{r.N_pred:.1f}" for _, r in mm.iterrows()))
    P("\ndifferential chi-square of observed vs expected counts in the 7 distance bins: see results.json (checks.census_differential_chi2)")
    P(f"\nBailer-Jones 2022 list: {s2['bj22_matched_in_candidates']} of {s2['bj22_rows']} recovered; median |d_lin - d_BJ22| = "
      f"{s2['bj22_within2Myr_median_abs_dd_pc']:.3f} pc for the {s2['bj22_within2Myr_N']} within 2 Myr")

    P("\n## Table 3. Exposure (distinct systems passing within d)\n")
    P("F(T), rate averaged over the last T Myr / present rate (orbit integration): "
      + ", ".join(f"{k} Myr: {v:.3f}" for k, v in s3["F_window"].items()))
    P(f"\norbit-factor range over potential and disc variants (2.4 Gyr average): {s3['F_orbit_min']:.2f}-{s3['F_orbit_max']:.2f}; "
      f"Martinez-Barbosa et al. 2017 scaled to 1 pc: " + ", ".join(f"{x:.1f}" for x in s3["martinez_barbosa_2017_gamma1pc"])
      + " per Myr (ratios " + ", ".join(f"{x:.2f}" for x in s3["martinez_barbosa_2017_ratio_to_no_migration"]) + ")\n")
    exp = pd.read_csv(os.path.join(DATA, "exposure_table.csv"))
    P("| era | T [Myr] | F | " + " | ".join(f"d={d:g} pc" for d in sorted(exp.d_pc.unique())) + " |")
    P("|---|---|---|" + "---|" * exp.d_pc.nunique())
    for era, g in exp.groupby("era", sort=False):
        cells = [f"{fmt(r.N_central)} ({fmt(r.N_low, 2)}-{fmt(r.N_high, 2)})" for _, r in g.sort_values("d_pc").iterrows()]
        P(f"| {era} | {g.T_Myr.iloc[0]:g} | {g.F_central.iloc[0]:.2f} | " + " | ".join(cells) + " |")
    P(f"\npassages within 1500 au per Gyr: {s3['passages_within_1500au_per_Gyr']:.2f}; systems inside the 100 light-year radio sphere: {s3['radio_sphere_systems_30p7pc']:.0f}")
    P(f"fraction of the encounter flux at the Sun slower than 30 / 10 km/s: {s3['sun_flux_frac_v_lt_30']:.3f} / {s3['sun_flux_frac_v_lt_10']:.4f}")

    P("\n## Table 4. Settlement runs (settleable fraction f = 1; fronts along U)\n")
    sim = pd.read_csv(os.path.join(DATA, "settlement_runs.csv"))
    P("| rule | d_hop | t_delay [Myr] | n d^3 | T2 kinetic [kyr] | T2 sim [kyr] (min-max) | seeds reaching 1000 | t(1 kpc) [Myr] | front 250-1000 pc [km/s] | travel frac |")
    P("|---|---|---|---|---|---|---|---|---|---|")
    for _, r in sim.iterrows():
        P(f"| {rule(r.dv_max_kms)} | {r.d_hop_pc:g} | {r.t_delay_Myr:g} | {r.eta:.2f} | "
          f"{1e3 * r.T2_kinetic_Myr:.0f} | {1e3 * r.T2_sim_Myr:.0f} ({1e3 * r.T2_sim_min:.0f}-{1e3 * r.T2_sim_max:.0f}) | "
          f"{int(r.n_reached)}/{int(r.n_real)} | {fmt(r.t_arr_1000pc)} | {fmt(r.front_kms)} | {fmt(r.travel_frac, 2)} |")
    for dv in sorted(sim.dv_max_kms.unique()):
        sub = sim[sim.dv_max_kms == dv]
        dev = (sub.T2_sim_Myr / sub.T2_kinetic_Myr - 1).abs()
        chk[f"T2_sim_vs_kinetic_max_dev_{'any' if np.isinf(dv) else 'dv30'}"] = float(dev.max())
        chk[f"T2_sim_vs_kinetic_n_within_15pct_{'any' if np.isinf(dv) else 'dv30'}"] = int((dev < 0.15).sum())
        chk[f"T2_cells_{'any' if np.isinf(dv) else 'dv30'}"] = int(len(sub))
        chk[f"seeds_reaching_1000_{'any' if np.isinf(dv) else 'dv30'}"] = f"{int(sub.n_reached.sum())}/{int(sub.n_real.sum())}"
        P(f"\n{rule(dv)}: T2 sim {1e3 * sub.T2_sim_Myr.min():.0f}-{1e3 * sub.T2_sim_Myr.max():.0f} kyr; "
          f"|sim/kinetic - 1| < 15% in {int((dev < 0.15).sum())} of {len(sub)} cells (worst {100 * dev.max():.0f}%); "
          f"seeds reaching 1000 systems: {int(sub.n_reached.sum())}/{int(sub.n_real.sum())}")

    g = s4["galaxy"]
    P("\n## Galaxy extrapolation\n")
    for k, v in g.items():
        P(f"- {k}: {fmt(v)}")
    P(f"- hop below which 1e11 systems cannot be reached in 1 Gyr (no delay, f = 1, local density): "
      f"{s4['d_crit_1Gyr_pc_dv30']:.3f} pc = {s4['d_crit_1Gyr_au_dv30']:.0f} au (dv<30 km/s), "
      f"{s4['d_crit_1Gyr_pc_any']:.3f} pc = {s4['d_crit_1Gyr_au_any']:.0f} au (any speed)")
    # growth-limited time for hops as short as Oort-cloud-grazing passages
    # (T_grow scales as d^-2 from the 1 Gyr critical hop computed in s4)
    AU_PC = 1.0 / 206264.8
    for d_au in (1500.0, 10000.0):
        tg = (s4["d_crit_1Gyr_pc_any"] / (d_au * AU_PC)) ** 2
        tg30 = (s4["d_crit_1Gyr_pc_dv30"] / (d_au * AU_PC)) ** 2
        chk[f"T_grow_1e11_Gyr_hop_{d_au:.0f}au_any"] = float(tg)
        chk[f"T_grow_1e11_Gyr_hop_{d_au:.0f}au_dv30"] = float(tg30)
        P(f"- growth-limited time to 1e11 systems for d_hop = {d_au:.0f} au: {tg30:.3g} Gyr (dv<30 km/s), "
          f"{tg:.3g} Gyr (any speed)")

    # What the straight-line box leaves out (assumption: flat rotation curve,
    # Vc = 233 km/s at R0 = 8.2 kpc).  A star carrying the front at radial
    # speed U oscillates over +-U/kappa instead of moving in a straight line,
    # and differential rotation shears a settled band of radial width dR at
    # Omega*dR.
    VC, R0 = 233.0, 8.2
    omega = VC / R0                       # km/s/kpc
    kappa = np.sqrt(2.0) * omega
    u_front = g["front_kms_median_U_dv30"]
    amp = u_front / kappa                 # kpc
    wrap = 2 * np.pi * R0 / (omega * 2 * amp * KMS_TO_PCMYR * 1e-3)   # Myr
    chk.update(dict(galaxy_omega_kms_per_kpc=float(omega), galaxy_kappa_kms_per_kpc=float(kappa),
                    galaxy_epicycle_amplitude_kpc_front_star=float(amp),
                    galaxy_shear_speed_kms_across_epicycle=float(omega * 2 * amp),
                    galaxy_shear_wrap_time_Myr=float(wrap),
                    galaxy_T_cross_Myr_times_pi_over_2=float(g["T_cross_Myr_median_U_dv30"] * np.pi / 2)))
    P(f"- beyond the box (assumed flat rotation curve): epicyclic frequency {kappa:.1f} km/s/kpc; a front star with "
      f"U = {u_front:.0f} km/s has epicycle amplitude {amp:.1f} kpc; shear across that epicycle {omega * 2 * amp:.0f} km/s; "
      f"azimuthal wrap time {wrap:.0f} Myr; radial crossing time with the mean outward speed of an epicycle "
      f"(2/pi of the peak): {g['T_cross_Myr_median_U_dv30'] * np.pi / 2:.0f} Myr")

    P("\n## Front-speed tests (no delay, f = 1)\n")
    ft = pd.read_csv(os.path.join(DATA, "front_tests.csv"))
    P("| rule | d_hop | case | axis | width [pc] | N stars | stars with |vx|>200 km/s | front [km/s] | +x side | -x side | travel frac |")
    P("|---|---|---|---|---|---|---|---|---|---|---|")
    for _, r in ft.iterrows():
        P(f"| {rule(r.dv_max_kms)} | {r.d_hop_pc:g} | {r.case} | {r.axis} | {r.width_pc:g} | {int(r.slab_N)} | {int(r.n_vx_gt_200)} | "
          f"{fmt(r.front_kms)} | {fmt(r.front_plus_kms)} | {fmt(r.front_minus_kms)} | {fmt(r.travel_frac, 2)} |")

    # fronts with the 30 km/s limit by hop length and direction (all f = 1 runs)
    allf = pd.concat([sim[(sim.dv_max_kms == 30) & sim.front_kms.notna()][["d_hop_pc", "axis", "front_kms"]],
                      ft[(ft.dv_max_kms == 30) & ft.front_kms.notna()][["d_hop_pc", "axis", "front_kms"]]])
    P("\nFront speed with the 30 km/s limit [km/s], range over runs (number of runs):\n")
    axes = ["U (radial)", "V (rotation)", "W (vertical)"]
    P("| d_hop [pc] | " + " | ".join(axes) + " |\n|---|---|---|---|")
    for d in sorted(allf.d_hop_pc.unique()):
        cells = []
        for ax in axes:
            v = allf[(allf.d_hop_pc == d) & (allf.axis == ax)].front_kms
            cells.append("n/a" if len(v) == 0 else (f"{v.min():.0f}" if len(v) == 1 else f"{v.min():.0f}-{v.max():.0f}") + f" ({len(v)})")
        P(f"| {d:g} | " + " | ".join(cells) + " |")

    # differential chi-square of the census against the expectation
    dN_obs = np.diff(np.r_[0, m.N_obs.values]); dN_pred = np.diff(np.r_[0.0, m.N_pred.values])
    chk["census_differential_chi2"] = float(np.sum((dN_obs - dN_pred) ** 2 / dN_pred))
    chk["census_differential_bins"] = int(len(dN_obs))

    P("\n## Time-step convergence of the doubling time\n")
    tc = pd.read_csv(os.path.join(DATA, "timestep_convergence.csv"))
    P("| rule | d_hop | step | dt [kyr] | T2 sim [kyr] | seeds reaching 1000 |\n|---|---|---|---|---|---|")
    for _, r in tc.iterrows():
        P(f"| {r.rule} | {r.d_hop_pc:g} | {r.step} | {r.dt_kyr:.1f} | {r.T2_sim_kyr:.0f} | {int(r.n_reached)} |")

    P("\n## Rendezvous limit: kinetic growth versus the velocity change a ship can make (f = 1, no delay)\n")
    kd = pd.read_csv(os.path.join(DATA, "kinetic_vs_dvcap.csv"))
    P("| dv limit [km/s] | share of pair-encounter flux | kinetic speed [km/s] (half-size sample) | T2 at 0.5 / 1 / 2 / 5 pc [Myr] | time to 1e11 systems at 1 pc [Myr] |")
    P("|---|---|---|---|---|")
    for _, r in kd.iterrows():
        P(f"| {'none' if np.isinf(r.dv_max_kms) else f'{r.dv_max_kms:g}'} | {r.pair_flux_frac:.4f} | {r.kinetic_speed_kms:.3g} ({r.kinetic_speed_kms_small_sample:.3g}) | "
          f"{fmt(r['T2_Myr_d0.5'])} / {fmt(r['T2_Myr_d1'])} / {fmt(r['T2_Myr_d2'])} / {fmt(r['T2_Myr_d5'])} | {fmt(r['T_grow_1e11_Myr_d1'])} |")
    vr = pd.read_csv(os.path.join(DATA, "voyager_class_runs.csv"))
    P("\nSimulated check, 10 km/s ships with a 10 km/s rendezvous limit (200 pc cube, growth from 10 to 1000 systems):\n")
    P("| d_hop | N stars | T2 kinetic [kyr] | T2 sim [kyr] (min-max) | seeds reaching 1000 |\n|---|---|---|---|---|")
    for _, r in vr.iterrows():
        P(f"| {r.d_hop_pc:g} | {int(r.N)} | {1e3 * r.T2_kinetic_Myr:.0f} | {1e3 * r.T2_sim_Myr:.0f} ({1e3 * r.T2_sim_min:.0f}-{1e3 * r.T2_sim_max:.0f}) | "
          f"{int(r.n_reached)}/{int(r.n_real)} |")

    P("\n## Settleable fraction f (dv < 30 km/s)\n")
    fk = pd.read_csv(os.path.join(DATA, "f_scaling_kinetic.csv"))
    for td in sorted(fk.t_delay_Myr.unique()):
        P(f"\nKinetic time to 1e11 systems [Myr], delay {td:g} Myr (rows: d_hop pc; columns: f)\n")
        pv = fk[(fk.dv_max_kms == 30) & (fk.t_delay_Myr == td)].pivot(index="d_hop_pc", columns="f", values="T_grow_1e11_Myr")
        pv = pv[sorted(pv.columns, reverse=True)]
        P("| d_hop [pc] | " + " | ".join(f"f={c:g}" for c in pv.columns) + " |")
        P("|---|" + "---|" * len(pv.columns))
        for d, row in pv.iterrows():
            P(f"| {d:g} | " + " | ".join(fmt(v) for v in row.values) + " |")
    P("\nKinetic doubling time [Myr], no delay (rows: d_hop pc; columns: f)\n")
    pv = fk[(fk.dv_max_kms == 30) & (fk.t_delay_Myr == 0)].pivot(index="d_hop_pc", columns="f", values="T2_Myr")
    pv = pv[sorted(pv.columns, reverse=True)]
    P("| d_hop [pc] | " + " | ".join(f"f={c:g}" for c in pv.columns) + " |")
    P("|---|" + "---|" * len(pv.columns))
    for d, row in pv.iterrows():
        P(f"| {d:g} | " + " | ".join(fmt(v) for v in row.values) + " |")
    fr = pd.read_csv(os.path.join(DATA, "f_scaling_runs.csv"))
    P("\nSimulated checks at f < 1 (no delay):\n")
    P("| kind | d_hop | f | T2 kinetic [Myr] | T2 sim [Myr] (min-max) | seeds reaching 1000 | front [km/s] | t(1 kpc) [Myr] | travel frac |\n|---|---|---|---|---|---|---|---|---|")
    for _, r in fr.iterrows():
        if r.kind == "cube":
            P(f"| cube | {r.d_hop_pc:g} | {r.f:g} | {fmt(r.T2_kinetic_Myr)} | {fmt(r.T2_sim_Myr)} ({fmt(r.T2_sim_min)}-{fmt(r.T2_sim_max)}) | "
              f"{int(r.n_reached)}/{int(r.n_real)} | | | |")
        else:
            P(f"| slab | {r.d_hop_pc:g} | {r.f:g} | {fmt(r.T2_kinetic_Myr)} | | | {fmt(r.front_kms)} | {fmt(r.t_arr_1000pc)} | {fmt(r.travel_frac, 2)} |")

    P("\n## Origins that stall under the 30 km/s limit\n")
    for k, v in s4["origins"].items():
        P(f"- {k}: {fmt(v)}")
    for d in ("1", "2"):
        o = s4[f"origins_sim_d{d}"]
        P(f"- simulated random origins, d_hop = {d} pc, no delay: {o['n_1000_within_30Myr']}/{o['n']} reach 1000 systems within 30 Myr, "
          f"{o['n_1000_within_100Myr']}/{o['n']} within 100 Myr; median origin speed {fmt(o['median_seed_speed_started'])} km/s (started) "
          f"vs {fmt(o['median_seed_speed_stalled'])} km/s (stalled)")

    mf = pd.read_csv(os.path.join(DATA, "meanfield_growth.csv"))
    P("\n## Kinetic growth time for 1e11 systems, ln(1e11)/s [Myr] (f = 1; rows: d_hop pc; columns: t_delay Myr)\n")
    for dv in sorted(mf.dv_max_kms.unique()):
        P(f"\nrendezvous limit: {'none' if np.isinf(dv) else f'{dv:g} km/s'}\n")
        pv = mf[mf.dv_max_kms == dv].pivot(index="d_hop_pc", columns="t_delay_Myr", values="T_grow_1e11_Myr")
        P("| d_hop [pc] | " + " | ".join(f"{c:g}" for c in pv.columns) + " |")
        P("|---|" + "---|" * len(pv.columns))
        for d, row in pv.iterrows():
            P(f"| {d:g} | " + " | ".join(fmt(v) for v in row.values) + " |")

    P("\n## Table 5. 95% limits on x\n")
    con = pd.read_csv(os.path.join(DATA, "fermi_constraint.csv"))
    for label in ("any speed", "v<30 km/s", "v<10 km/s"):
        P(f"\npassers: {label}\n")
        pv = con[con.passers == label].pivot(index="persistence_Myr", columns="d_pc", values="x95")
        P("| persistence [Myr] | " + " | ".join(f"d={d:g} pc" for d in pv.columns) + " |")
        P("|---|" + "---|" * len(pv.columns))
        for L, row in pv.iterrows():
            P(f"| {L:g} | " + " | ".join(fmt(v, 2) for v in row.values) + " |")
    P("\npassers: any speed, conservative (low end of rate, era-independent low F)\n")
    pv = con[con.passers == "any speed"].pivot(index="persistence_Myr", columns="d_pc", values="x95_conservative")
    P("| persistence [Myr] | " + " | ".join(f"d={d:g} pc" for d in pv.columns) + " |")
    P("|---|" + "---|" * len(pv.columns))
    for L, row in pv.iterrows():
        P(f"| {L:g} | " + " | ".join(fmt(v, 2) for v in row.values) + " |")

    P(f"\nStatic percolation: d50 real {s4['static']['R40_d50_real']:.2f} pc, Poisson {s4['static']['R40_d50_poisson']:.2f} pc "
      f"(N = {s4['static']['R40_N']} systems within 40 pc); infinite-field threshold {s4['static']['d_c_infinite_poisson_adopted_density']:.2f} pc")
    P(f"Simulator validation: {s4['validation_contacts_sim']} contacts vs {s4['validation_contacts_expected']:.1f} expected "
      f"({s4['validation_ratio']:.3f} +- {s4['validation_ratio_poisson_err']:.3f})")
    P(f"Pair statistics: mean relative speed {s4['pair_mean_vrel_kms']:.1f} km/s, median {s4['pair_median_vrel_kms']:.1f} km/s; "
      f"share of pair flux below 30 km/s {s4['pair_flux_frac_dv30']:.3f}; kinetic speed {s4['kinetic_speed_dv30_kms']:.2f} km/s (dv<30), "
      f"{s4['kinetic_speed_any_kms']:.1f} km/s (any)")

    txt = "\n".join(lines)
    with open(os.path.join(DATA, "report_tables.md"), "w") as f:
        f.write(txt + "\n")
    print(txt)
    save_json("results.json", dict(s1=s1, s2=s2, s3=s3, s4=s4, checks=chk))


if __name__ == "__main__":
    main()