"""Write report.md from the pipeline outputs, so that every number in the text is the one in the tables.

The prose template lives in this file; every quantity marked <<name>> is filled from data/key_numbers.json,
data/hi_line_scenarios.csv, data/dyson_limits.csv, data/grimaldi_comparison.csv, data/surveys_compilation.csv
and data/radio_model_components.csv.  A missing or unused token is an error.
"""
import json
import os
import re

import numpy as np
import pandas as pd

from common import DATA, STUDY, L_KII


def num(x, sf=2):
    """Counts: 2 significant figures; thousands separators below 1e7, 'a.b x 10^n' above."""
    if x is None or (isinstance(x, float) and np.isnan(x)):
        return "n/a"
    if np.isinf(x):
        return "unbounded"
    if x == 0:
        return "0"
    e = int(np.floor(np.log10(abs(x))))
    r = round(x, -e + sf - 1)
    if abs(r) < 10:
        return f"{r:.{sf - 1}f}"
    if abs(r) < 1e7:
        return f"{int(round(r)):,}"
    m = r / 10 ** e
    if round(m, sf - 1) >= 10:
        m, e = m / 10, e + 1
    return f"{m:.{sf - 1}f} × 10^{e}"


def sci(x, sf=2):
    """Small probabilities / fractions."""
    if x == 0:
        return "0"
    if x >= 0.01:
        return f"{x:.{sf}g}"
    e = int(np.floor(np.log10(x)))
    m = round(x / 10 ** e, sf - 1)
    if m >= 10:
        m, e = m / 10, e + 1
    return f"{m:.{sf - 1}f} × 10^{e}"


def pct(x, sf=2):
    v = 100 * x
    return (f"{v:.0f}%" if v >= 10 else f"{v:.{sf}g}%")


TEMPLATE = r"""# What 65 years of null technosignature searches actually exclude: a unified limit

*Revised after referee review (2026-10-12). This file is written by `code/make_report.py`; every number is read from the tables produced by `code/run_all.py`.*

**Headline.** I merged <<n_comp_base>> published radio survey components (from a <<n_rows>>-row compilation, 1960-2026) into one coverage function: the number of Milky Way stars (model total <<NMW>>) at which a continuous narrowband transmitter of given power would have been detected. All limits are 95% upper limits on the number N of such transmitters operating now, at detection efficiency 0.5; bracketed ranges span efficiency 0.25-1 and every robustness variant. (1) Always-on isotropic 10^17 W beacons somewhere in 1-10 GHz: **N < <<n17_110>>** (<<n17_110_lo>> to <<n17_110_hi>>). (2) The same on the 21-cm line: **N < <<n17_hi>>**, but only if the 1980s-90s sky surveys META and META II are trusted at nominal thresholds and the carrier arrives free of drift; the defensible range is **<<n17_hi_lo>> to <<n17_hi_hi>>** (<<n17_hi_sah>> if SETI@home's unresolved candidates are all rejected; <<n17_hi_surf>> for a transmitter on a rotating planet if they are not). (3) Kardashev-II (3.8 × 10^26 W) beacons: **N < <<kii_hi>>** on the 21-cm line, which is only the efficiency floor (about <<kii_hi_c>> if just SETI@home and post-2015 searches are trusted; <<kii_hi_j>> to <<kii_hi_gc>> if neither legacy surveys nor SETI@home count), and **N < <<kii_110>>** (<<kii_110_lo>> to <<kii_110_hi>>) anywhere in 1-10 GHz, because dedicated narrowband searches have covered only <<cov_kii_110>> of the stars × log-frequency plane, <<beta_share>> of that through one survey (BETA) at a nominal threshold; post-2015 searches alone cover <<cov_kii_110_mod>>. (4) Arecibo-class radars and Earth-level leakage are unconstrained: the limits exceed the number of stars. (5) Dyson swarms: the only published limits are < 4.8%, < 0.048% and < 0.0086% of stars within 300 pc for covering fractions above 10%, 50% and 90% (Suazo et al. 2022). **Main caveats:** transmitters are assumed distributed like stars, always on and narrower than a few Hz; the 21-cm and K-II coverage rests mostly on legacy surveys whose thresholds I could check only in the Technosearch database; SETI@home is left out of the baseline because its final candidates have no published disposition; in-beam star counts are tied to Mason et al. (2025); wide-band continuum surveys, which would probably have seen a K-II beacon, are not included.

## 1. Question

What does the whole null record, rather than one survey at a time, exclude as a function of transmitter power, frequency, duty cycle and beaming? How many civilizations could the Galaxy contain that behave as (i) always-on isotropic beacons, (ii) Arecibo-like beamed radars, or (iii) Earth in 2025 (leakage only), without contradicting any published null?

## 2. Data

- **Compilation** (`data/surveys_compilation.csv`): <<n_rows>> rows (<<n_radio>> radio, <<n_opt>> optical, <<n_ir>> infrared, <<n_art>> artefact), typed by hand; every row states a result. <<n_full>> rows were checked against full text, <<n_abs>> against the arXiv abstract only, <<n_ts>> only against the Technosearch database, <<n_oth>> against a secondary table or listing. All <<n_arxiv>> arXiv identifiers resolve (`data/arxiv_verification.csv`). Quantitative columns are sparse: dwell time is filled in <<fill_dwell>> of rows, flux threshold in <<fill_smin>>, drift range in <<fill_drift>>.
- **Model inputs** (`data/radio_model_components.csv`): <<n_comp>> components: <<n_tierA>> tier A (numbers from the papers, mostly post-2015), <<n_tierB>> tier B (legacy, nominal thresholds), <<n_tierX>> extragalactic (never counted as Milky Way stars).
- **Downloads** (`code/fetch_data.py`): Technosearch lists (<<n_tsr>> radio, <<n_tso>> optical, <<n_tsa>> archival entries); VizieR tables of Wlodarczyk-Sroka et al. (2020), Harp et al. (2016), Isaacson et al. (2017), Baumgardt & Hilker (2018) and the Gaia Catalogue of Nearby Stars.
- **Attrition.** The Breakthrough Listen target lists recovered from VizieR hold <<bl_L_uniq>> unique stars at GBT L band (<<bl_L_rows>> rows), <<bl_S_uniq>> at S band (<<bl_S_rows>>) and <<bl_P_uniq>> at Parkes (<<bl_P_rows>>), against 882, 1005 and 189 in Price et al. (2020); the rest are not modelled.

## 3. Method

1. **Coverage.** For each component (flux threshold S_min, band, footprint), N_cov(L) is the number of stars in the footprint closer than sqrt(L / 4π S_min). Stars are named targets with distances, or come from an analytic Galaxy model (two discs, halo, bulge) normalised to the Gaia 100 pc sample: <<NMW>> objects, about <<Mstar>> solar masses, some 25% below the usual 5-6 × 10^10. Limits scale at most linearly with this.
2. **Union and priors.** At each frequency, components are merged zone by zone (maximum where footprints overlap, sum where they do not). N_eff(L) averages N_cov over a frequency prior: the 21-cm line as received in our frame (1420.2-1420.6 MHz, META's window); the 21-cm line in the star's rest frame (sky surveys then count only stars whose velocity falls in their window: <<vf_meta>> for META, <<vf_sah>> for SETI@home); the water hole; log-uniform 1-10 GHz and 0.1-100 GHz.
3. **Likelihood.** One transmitter is detected with probability p = η × duty × beam fraction × N_eff / N_MW. The all-null record has probability (1 - p)^N exactly, so N_max = ln 0.05 / ln(1 - p); `code/seti_likelihood.py` returns N ln(1 - p).
4. **Drift.** Each survey's threshold is raised by max(1, drift × integration time / channel width) for three transmitter classes: drift-free (in free fall far from a star, or Doppler-corrected towards us), Earth-like orbit, Earth-like rotating surface. META and META II (0.05 Hz, 20 s) lose factors of <<meta_orbit>> and <<meta_surf>>; Breakthrough Listen <<bl_L_surf>> to <<bl_X_surf>>. Beacons default to drift-free, radars and leakage to the surface class.
5. **In-beam stars** of each Breakthrough Listen pointing are counted with a Gaussian beam taper and rescaled, in five power shells, to the Besançon counts of Mason et al. (2025); unscaled, my model gives <<raw_gbt_lo>> to <<raw_gbt_hi>> times their GBT counts.
6. **SETI@home** (Anderson et al. 2025) is excluded from the baseline: about 200 final candidates, 92 re-observed with FAST, had no published disposition by October 2026. A switch counts it as a completed null over 30% of the sky.
7. **Dyson swarms.** I count stars around which a blackbody swarm would be a > 5σ WISE W3+W4 source: a sensitivity count, not a limit.

## 4. Results

**Table 1. Isotropic always-on beacons, baseline record (no SETI@home, drift-free, η = 0.5).**

| Frequency prior | N_eff at 10^17 W | N_max at 10^17 W (range) | N_max at K-II (range) | Stars covered at K-II |
|---|---|---|---|---|
<<table1>>

At 10^17 W, META plus META II supply <<meta_share>> of the 21-cm N_eff. In 1-10 GHz, <<inbeam_share>> of N_eff is in-beam stars of Breakthrough Listen pointings, and these are concentrated: ten pointings carry <<top10_L>> (L band) and <<top10_S>> (S band) of the count, <<blt10_L>> lies within 10° of the Galactic plane, and the median star is <<med_L>> kpc away. In 0.1-100 GHz, <<ser_share>> of N_eff is one survey, SERENDIP III at 424-436 MHz, at a nominal threshold. Above about 10^21 W every curve is flat (Figs. 2 and 3): unsearched frequencies, not sensitivity, set the limit.

**Table 2. What the limits depend on (N_max, η = 0.5).**

| Scenario | 10^17 W, 21 cm our frame | 10^17 W, 21 cm star frame | 10^17 W, 1-10 GHz | K-II, 21 cm our frame | K-II, 1-10 GHz |
|---|---|---|---|---|---|
<<table2>>

The 10^17 W 21-cm limit moves by a factor of <<hi_span>> across these rows, the 1-10 GHz limit by <<span_110>>. The K-II 21-cm limit sits at the floor only while the legacy sky surveys count: a given K-II 21-cm beacon had a <<p_kii_a>> chance of being caught if they do, <<p_kii_c>> if only SETI@home and modern searches do, and <<p_kii_j>> if neither does. With the Galactic-centre star counts of Gajjar et al. (2021) instead of mine, the K-II entries of row J become <<kii_hi_gc>> and <<kii_110_gc>>.

**Table 3. The three behaviours (1-10 GHz prior).**

| Behaviour | EIRP (W) | Duty × beam fraction | p per civilization | N_max (range) | Detections expected if every star hosted one |
|---|---|---|---|---|---|
<<table3>>

Even aimed at the Sun, an Arecibo-class radar population is limited only to <<n_radar_aimed>>, within a factor of a few of the number of stars. Leakage is given two ways. Literal: the 4 GW peak of Saide et al. (2023) is summed over all LTE bands; spread over 20 MHz (my assumption, not theirs) it leaves <<lte_chan>> W per 2.79 Hz channel, visible to FAST out to <<lte_au>> AU, so no star is covered. Generous: all 4 GW in one carrier would reach <<lte_pc>> pc, covering <<lte_nstar>> FAST targets.

**Posterior.** With a log-uniform prior on N from 1 to N_MW (prior 95th percentile <<prior_q95>>), the posterior 95th percentile is <<post17_hi>> for 10^17 W 21-cm beacons, <<post17_110>> for 10^17 W beacons in 1-10 GHz, <<postkii_110>> for K-II beacons in 1-10 GHz, and essentially the prior for radars and leakage (Fig. 4).

**Table 4. Dyson swarms at 300 K: fraction of stars.**

| Covering fraction above | Published limit, < 100 pc | Published limit, < 300 pc | Published limit, < 5 kpc | Sensitivity floor, < 300 pc (not a limit) |
|---|---|---|---|---|
<<table4>>

The published limits are the photometric-compatibility fractions of Suazo et al. (2022). The floor is what a fully complete search of the <<dy_sens_lo>> to <<dy_sens_hi>> million stars within 300 pc around which WISE could see such a swarm would give. No published pipeline has measured its completeness, and two of the seven candidates of Suazo et al. (2024) with fitted covering fraction above 0.1 are not yet refuted, so the 10% floor allows for two survivors.

Figures (in `figures/`): `fig1_transmitter_fraction_limit.png` (limit on the fraction of stars with transmitters against power, with published points), `fig2_coverage_map.png` (frequency × power coverage), `fig3_max_civilizations.png` (N_max against power; scenario ladder), `fig4_posteriors.png`, `fig5_dyson_limits.png`.

## 5. Robustness

N_max for 10^17 W beacons under each variant (`data/robustness_variants.csv`):

- **De-duplication:** 1-10 GHz <<v_lo_110>> (plain sum) to <<v_hi_110>> (maximum only).
- **In-beam counts:** unscaled model <<v_raw_110>>; calibration without the two densest inner-Galaxy pointings <<v_ex2_110>>; all in-beam counts halved <<v_half_110>>.
- **Efficiency:** Breakthrough Listen at the 12.7% injection recovery of Margot et al. (2023): <<v_eff_110>>. All thresholds doubled: <<v_x2_110>>.
- **Legacy thresholds:** tier B × 10 moves the 21-cm limit to <<v_tb_hi>> and the 0.1-100 GHz limit from <<n17_wide>> to <<v_tb_wide>>; META × 10, 30, 60 gives <<v_m10>>, <<v_m30>>, <<v_m60>>.
- **Drift:** the surface class cuts 21-cm coverage to <<dr_hi>> of the drift-free value and 1-10 GHz coverage to <<dr_110>>.
- **SETI@home sky fraction** 25% instead of 30%: <<v_sah25>> against <<n17_hi_sah>>. Scintillation weighting: <<v_scint_110>> in 1-10 GHz.
- **Galactic-centre fields:** my model holds <<gc_ratio>> times the stars quoted by Gajjar et al. (2021); this matters only where legacy surveys are dropped (Table 2).

## 6. Limitations

- Continuous narrowband transmitters only; pulsed, swept or broadband signals enter only through the duty factor.
- Transmitters follow the stars. A population concentrated away from the searched sight lines is less constrained.
- Legacy surveys enter at nominal thresholds, without drift search and, for META, with a frame assumption; I could not read their original papers. They carry the 21-cm and K-II results.
- **Radio astronomy at large is not included.** A K-II narrowband beacon 8-25 kpc away would be a <<nvss_far>> to <<nvss_near>> Jy source in NVSS and <<gb6_far>> to <<gb6_near>> Jy in GB6/PMN, so continuum surveys have probably excluded it at 1.4 and 4.85 GHz already. I have not quantified this (interferometric flagging could remove such a line), so the true K-II constraint is probably tighter than Table 1. Lacki (2024, 2025) treats non-SETI bounds; his collective bound of 10^6 (L / 10^13 W)^-1 L-band broadcasts per star per GHz is, by my arithmetic, about <<lacki_kii>> K-II beacons per GHz, weaker than the limits here.
- The Parkes in-beam calibration needs factors from <<pk_fac_hi>> (near) to <<pk_fac_lo>> (far), which I cannot explain; Parkes carries under 10% of the in-beam count.
- Optical, artefact and most infrared searches are compiled, not modelled.

## 7. Closest prior work and what is new

Per-survey transmitter fractions: Enriquez et al. (2017), Price et al. (2020), Margot et al. (2023). In-beam stars: Wlodarczyk-Sroka et al. (2020) and Mason et al. (2025); my in-beam counting re-derives theirs and is calibrated to it, so it is not new. Search-volume fractions: Tarter et al. (2010), Wright et al. (2018). Population inference: Grimaldi & Marcy (2018), Grimaldi (2023), Lacki (2024, 2025). Dyson limits: Suazo et al. (2022, 2024).

What I believe is new: (1) one cross-survey coverage function in units of Milky Way stars, de-duplicated across surveys, with an exact likelihood module; (2) the conditional structure of the 21-cm limit (SETI@home status × transmitter drift × legacy thresholds, Table 2); (3) for dedicated SETI searches only, the statement that above about 10^21 W frequency coverage rather than sensitivity sets the limit.

**Comparison with Grimaldi & Marcy (2018)**, uniform 1-2 GHz prior: the sampled fraction of the Galaxy in their idealised prior search, 2.6 × 10^-8 (L / 2 × 10^13 W)^1.5, against my coverage fraction.

| EIRP (W) | Grimaldi & Marcy | This work | Ratio |
|---|---|---|---|
<<table5>>

Mine is lower because theirs assumes the whole sky searched to 10^-23 W/m² at every frequency in the band. Grimaldi (2023) bounds the emission rate at under 1-5 per century on the premise that no detectable emission has crossed Earth in 60 years. My coverage supports that premise only for 21-cm beacons above about 10^20 W, and only if the legacy surveys are trusted (<<cov1e20_hi>> of stars covered, against <<cov1e20_110>> in 1-10 GHz). His quantity is a birth rate that depends on signal longevity, so I give no like-for-like number.

## 8. Implications for the Fermi Paradox

The radio "Great Silence" is a narrow observational fact. It excludes a Galaxy holding more than a handful of deliberate K-II 21-cm beacons, more than <<kii_110>> (legacy surveys trusted) to <<kii_110_mod>>-<<kii_110_gc>> (post-2015 searches only) K-II beacons elsewhere in 1-10 GHz, or more than <<n17_hi_lo>> to <<n17_hi_hi>> planet-scale (10^17 W) 21-cm beacons. It does not exclude a Galaxy in which every star hosts a civilization like ours: the whole record would then expect <<mu_leak>> detections under the generous leakage bound, none under the literal one, and <<mu_radar>> from randomly pointed planetary radars. The silence therefore constrains behaviour, not abundance: if a fraction f of civilizations run 10^17 W isotropic beacons somewhere in 1-10 GHz, then N × f < <<n17_110>>. An integrator should call `seti_likelihood.log_likelihood` with an explicit behaviour mixture and treat the drift class and the SETI@home switch as nuisance choices; for any population without loud, persistent, omnidirectional beacons the radio term is flat. The paradox is then carried by other evidence: the absence of settlement or artefacts, which this study only compiles, and waste heat, where fewer than about 1 in 10^4 nearby stars can host a near-complete swarm.

## 9. Response to referee

1. **No report.** Written; generated by the pipeline so that text and tables cannot diverge.
2. **21-cm limit conditional.** Accepted. SETI@home is excluded by default; Table 2 gives <<n17_hi_sah>> (SETI@home counted), <<n17_hi>> (META at nominal thresholds, drift-free) and about 10^5 otherwise; the quoted range runs to <<n17_hi_hi>>.
3. **Drift as a per-survey factor.** Accepted; computed from each survey's channel width and integration time (`data/drift_factors.csv`), with SETI@home's non-barycentric sensitivities for drifting transmitters.
4. **K-II 21 cm.** Accepted in substance: the limit is stated as the efficiency floor with the exact binomial value <<kii_hi>>, and the 5%-50% detection chance is given. Two of the referee's numbers are not adopted. In the star's rest frame the K-II limit is <<kii_rest>>, not about 16, because BETA's 1400-1720 MHz band contains every Doppler shift; <<n1e20_rest>> is what I get at 10^20 W. And "4 to 60" is not the full range: the limit is <<kii_hi_j>> or more if neither legacy surveys nor SETI@home count.
5. **1-10 GHz at 10^17 W.** Accepted. The baseline is now calibrated shell by shell to Mason et al.; the central value is <<n17_110>> (the referee's rebuild gave 4.8 × 10^5) and the range is kept. Like for like, with the beam taper, my unscaled model is <<raw_gbt_lo>> to 0.9 times their GBT counts inside 3 kpc and <<raw_gbt_mid>> to <<raw_gbt_hi>> times beyond, which confirms the excess at large distance. The pointing-concentration sentence is in Section 4.
6. **K-II 1-10 GHz.** Accepted: coverage statement reworded, limit quoted as <<kii_110>> to <<kii_110_mod>> (<<kii_110_gc>> with Gajjar et al.'s star counts), continuum caveat and Lacki citations added. The continuum constraint is stated, not quantified.
7. **Dyson.** Accepted: the published Suazo et al. (2022) fractions lead; the 10^-6 figures are labelled a sensitivity floor; the seven-candidate column is gone; two survivors are allowed in the 10% bin; `p_detect_dyson` requires an explicit completeness.
8. **Likelihood.** Accepted: `log_likelihood` returns N ln(1 - p); the SETI@home treatment is documented and switchable (default off, stricter than requested).
9. **2 × 10^13 W rows.** Fixed: every value is interpolated at the stated power.
10. **SERENDIP III.** Stated in Section 4 (<<ser_share>> on recomputation; the referee found 85%).
11. **LTE bandwidth.** Labelled as my assumption.
12. **Fig. 1.** Redrawn at η = 0.5 with numbered, non-overlapping labels.
13. **Compilation.** Result column filled for all rows, structured infrared columns added, verification levels stated. The "69 rows without a result" came from pandas reading the string "null" as missing; the rows now say "no detection".
14. **Novelty.** Claim (4) removed, claim (3) restricted to dedicated SETI searches, in-beam counting described as a re-derivation of Mason et al.
15. **Grimaldi.** Numerical comparison with Grimaldi & Marcy (2018) added; for Grimaldi (2023) I explain why none is made.

Minor points: target-list attrition, the star-count normalisation and the SETI@home sky fraction (30%, variant 25%) are now stated.
"""


def build_tokens():
    K = json.load(open(os.path.join(DATA, "key_numbers.json")))
    comp = pd.read_csv(os.path.join(DATA, "survey_contributions.csv"))
    rawc = pd.read_csv(os.path.join(DATA, "radio_model_components.csv"))
    cmp_ = pd.read_csv(os.path.join(DATA, "surveys_compilation.csv"), keep_default_na=False)
    lad = pd.read_csv(os.path.join(DATA, "hi_line_scenarios.csv"))
    dy = pd.read_csv(os.path.join(DATA, "dyson_limits.csv"))
    gm = pd.read_csv(os.path.join(DATA, "grimaldi_comparison.csv"))
    val = pd.read_csv(os.path.join(DATA, "galaxy_model_validation.csv"))
    ax = pd.read_csv(os.path.join(DATA, "arxiv_verification.csv"))
    B, R, NMW = K["behaviours"], K["beacon_overall_range"], K["N_MW_model"]
    PR = [("HI_line", "21-cm line, our frame"), ("HI_rest", "21-cm line, star's rest frame"), ("water_hole", "Water hole, 1.42-1.67 GHz"),
          ("log_1_10GHz", "1-10 GHz, log-uniform"), ("log_0p1_100GHz", "0.1-100 GHz, log-uniform")]
    t = {}
    b17, bk = B["i_beacon_1e17W"], B["i_beacon_KII"]
    t["NMW"] = num(NMW)
    t["Mstar"] = num(K["stellar_mass_equiv_Msun"])
    t["n_comp"] = str(len(rawc))
    t["n_comp_base"] = str(int((comp.in_baseline_record & (comp.tier != "X")).sum()))
    t["n_tierA"], t["n_tierB"], t["n_tierX"] = (str(int((rawc.tier == x).sum())) for x in "ABX")
    # ---- compilation statistics
    t["n_rows"] = str(len(cmp_))
    for k, d in (("n_radio", "radio"), ("n_opt", "optical"), ("n_ir", "infrared"), ("n_art", "artifact")):
        t[k] = str(int((cmp_.domain == d).sum()))
    assert (cmp_.result != "").all()
    vf = cmp_.verified_from
    full, abs_only, ts_only = vf.str.contains("full text"), vf.isin(["abstract", "abstract (first lines)"]), vf == "Technosearch"
    t["n_full"], t["n_abs"], t["n_ts"] = str(int(full.sum())), str(int(abs_only.sum())), str(int(ts_only.sum()))
    t["n_oth"] = str(int(len(cmp_) - full.sum() - abs_only.sum() - ts_only.sum()))
    t["n_arxiv"] = str(int(ax.arxiv_resolves.sum()))
    assert ax.arxiv_resolves.all()
    for k, c in (("fill_dwell", "dwell_s"), ("fill_smin", "smin_Wm2"), ("fill_drift", "max_drift_Hz_s")):
        t[k] = pct(float((cmp_[c] != "").mean()))
    for k, f in (("n_tsr", "technosearch_radio.csv"), ("n_tso", "technosearch_optical.csv"), ("n_tsa", "technosearch_archival.csv")):
        t[k] = str(len(pd.read_csv(os.path.join(DATA, f))))
    for tag, key in (("L", "bl_gbtL_targets"), ("S", "bl_gbtS_targets"), ("P", "bl_pks10_targets")):
        t[f"bl_{tag}_rows"], t[f"bl_{tag}_uniq"] = (str(int(x)) for x in K["bl_target_rows"][key])
    # ---- method numbers
    t["vf_meta"] = pct(K["velocity_window_fraction"]["meta"])
    t["vf_sah"] = pct(K["velocity_window_fraction"]["setiathome"])
    dtf = K["drift_threshold_factors"]
    t["meta_orbit"], t["meta_surf"] = f"{dtf['orbit']['meta']:.0f}", f"{dtf['surface']['meta']:.0f}"
    t["bl_L_surf"], t["bl_X_surf"] = f"{dtf['surface']['bl_gbtL_targets']:.1f}", f"{dtf['surface']['bl_gbtX_archive']:.1f}"
    rg = K["bycatch_raw_over_mason"]["GBT"]
    t["raw_gbt_lo"], t["raw_gbt_mid"], t["raw_gbt_hi"] = f"{min(rg):.2f}", f"{rg[3]:.1f}", f"{max(rg):.1f}"
    assert max(rg[:3]) < 0.95 and rg[3] > 1
    fp = K["bycatch_shell_factors"]["Parkes"]
    t["pk_fac_hi"], t["pk_fac_lo"] = f"{max(fp):.1f}", f"{min(fp):.2f}"
    # ---- headline
    d17 = {p: b17[p]["inertial"] for p, _ in PR}
    t["n17_110"], t["n17_110_lo"], t["n17_110_hi"] = num(d17["log_1_10GHz"]["Nmax95"]), num(R["i_beacon_1e17W"]["log_1_10GHz"]["lo"]), num(R["i_beacon_1e17W"]["log_1_10GHz"]["hi"])
    t["n17_hi"], t["n17_hi_lo"], t["n17_hi_hi"] = num(d17["HI_line"]["Nmax95"]), num(R["i_beacon_1e17W"]["HI_line"]["lo"]), num(R["i_beacon_1e17W"]["HI_line"]["hi"])
    t["n17_hi_sah"] = num(b17["HI_line"]["inertial_sah"]["Nmax95"])
    t["n17_hi_surf"] = num(b17["HI_line"]["surface"]["Nmax95"])
    t["n17_wide"] = num(d17["log_0p1_100GHz"]["Nmax95"])
    t["kii_hi"] = num(bk["HI_line"]["inertial"]["Nmax95"])
    t["kii_rest"] = num(bk["HI_rest"]["inertial"]["Nmax95"])
    t["n1e20_rest"] = num(B["i_beacon_1e20W"]["HI_rest"]["inertial"]["Nmax95"])
    t["kii_110"], t["kii_110_lo"], t["kii_110_hi"] = num(bk["log_1_10GHz"]["inertial"]["Nmax95"]), num(R["i_beacon_KII"]["log_1_10GHz"]["lo"]), num(R["i_beacon_KII"]["log_1_10GHz"]["hi"])
    vk, v17 = K["variants_KII"]["inertial"], K["variants_1e17"]["inertial"]
    name = lambda frag: [n for n in v17 if frag in n][0]
    MOD, MODS, MODGC = name("modern surveys only (legacy"), name("modern surveys + SETI@home"), name("modern surveys only, Galactic")
    t["kii_hi_c"], t["kii_hi_j"], t["kii_hi_gc"] = num(vk[MODS]["HI_line"]), num(vk[MOD]["HI_line"]), num(vk[MODGC]["HI_line"])
    t["kii_110_mod"], t["kii_110_gc"] = num(vk[MOD]["log_1_10GHz"]), num(vk[MODGC]["log_1_10GHz"])
    cf = K["coverage_fraction"]
    t["cov_kii_110"], t["cov_kii_110_mod"] = pct(cf["base"]["log_1_10GHz"]), pct(cf["modern"]["log_1_10GHz"])
    sh = K["shares"]
    t["beta_share"], t["meta_share"] = pct(sh["beta_of_KII_1_10GHz"]), pct(sh["meta_meta2_of_1e17_HI"])
    t["inbeam_share"], t["ser_share"] = pct(sh["inbeam_of_1e17_1_10GHz"]), pct(sh["serendip3_of_1e17_0p1_100GHz"])
    bp = K["bycatch_pointings"]
    t["top10_L"], t["top10_S"] = pct(bp["bl_gbtL_bycatch"]["top10_share"]), pct(bp["bl_gbtS_bycatch"]["top10_share"])
    t["blt10_L"], t["med_L"] = pct(bp["bl_gbtL_bycatch"]["share_absb_lt10"]), f"{bp['bl_gbtL_bycatch']['median_d_pc'] / 1e3:.0f}"
    # ---- Table 1
    rows = []
    for p, lab in PR:
        r17, rk = R["i_beacon_1e17W"][p], R["i_beacon_KII"][p]
        lo_k = "< 1" if rk["lo"] < 1 else num(rk["lo"])
        rows.append(f"| {lab} | {num(d17[p]['Neff'])} | {num(d17[p]['Nmax95'])} ({num(r17['lo'])} to {num(r17['hi'])}) | "
                    f"{num(bk[p]['inertial']['Nmax95'])} ({lo_k} to {num(rk['hi'])}) | {pct(cf['base'][p])} |")
    t["table1"] = "\n".join(rows)
    # ---- Table 2
    rows, hi_vals, v110 = [], [], []
    for _, r in lad.iterrows():
        hi_vals.append(r["Nmax95_HI_line_1e17W"]); v110.append(r["Nmax95_log_1_10GHz_1e17W"])
        rows.append(f"| {r['scenario']} | {num(r['Nmax95_HI_line_1e17W'])} | {num(r['Nmax95_HI_rest_1e17W'])} | {num(r['Nmax95_log_1_10GHz_1e17W'])} | "
                    f"{num(r['Nmax95_HI_line_KII'])} | {num(r['Nmax95_log_1_10GHz_KII'])} |")
    t["table2"] = "\n".join(rows)
    t["hi_span"], t["span_110"] = f"{max(hi_vals) / min(hi_vals):.0f}", f"{max(v110) / min(v110):.1f}"
    pk = lambda n: 1.0 - 0.05 ** (1.0 / n)
    t["p_kii_a"], t["p_kii_c"], t["p_kii_j"] = pct(pk(vk[name("baseline")]["HI_line"])), pct(pk(vk[MODS]["HI_line"])), pct(pk(vk[MOD]["HI_line"]))
    # ---- Table 3
    spec = [("i_beacon_1e17W", "(i) isotropic beacon, 10^17 W", "1"), ("i_beacon_KII", "(i) isotropic beacon, K-II", "1"),
            ("i_beacon_2e13W", "(i) isotropic beacon, Arecibo-radar power", "1"),
            ("ii_radar_random", "(ii) Arecibo-class beam, random pointing", "0.05 × 5 × 10^-8"),
            ("ii_radar_aimed", "(ii') same, aimed at the Sun", "0.05 × 1"),
            ("iii_leak_narrow", "(iii) Earth-2025 leakage, generous: 4 GW in one carrier", "1"),
            ("iii_leak_literal", "(iii) Earth-2025 leakage, literal: 4 GW over 20 MHz", "1")]
    eirp = {"i_beacon_1e17W": "10^17", "i_beacon_KII": "3.8 × 10^26", "i_beacon_2e13W": "2 × 10^13", "ii_radar_random": "2 × 10^13",
            "ii_radar_aimed": "2 × 10^13", "iii_leak_narrow": "4 × 10^9", "iii_leak_literal": f"{K['lte']['L_per_channel_W']:.0f} per channel"}
    rows = []
    for bid, lab, db in spec:
        x = B[bid]["log_1_10GHz"][K["default_drift"][bid]]
        ds = ("inertial", "surface") if K["default_drift"][bid] == "inertial" else ("surface",)
        lo = min(B[bid]["log_1_10GHz"][d]["Nmax95_lo"] for d in ds)
        hi = max(B[bid]["log_1_10GHz"][d]["Nmax95_hi"] for d in ds)
        rng = "unbounded" if np.isinf(x["Nmax95"]) else f"{num(x['Nmax95'])} ({num(lo)} to {num(hi)})"
        mu = x["mu_all"]
        rows.append(f"| {lab} | {eirp[bid]} | {db} | {sci(x['p'])} | {rng} | {num(mu) if mu >= 10 else sci(mu)} |")
    t["table3"] = "\n".join(rows)
    xr, xa, xl = (B[b]["log_1_10GHz"]["surface"] for b in ("ii_radar_random", "ii_radar_aimed", "iii_leak_narrow"))
    t["n_radar_aimed"] = num(xa["Nmax95"])
    t["mu_leak"], t["mu_radar"] = f"{xl['mu_all']:.1f}", sci(xr["mu_all"])
    lte = K["lte"]
    t["lte_chan"], t["lte_au"] = f"{lte['L_per_channel_W']:.0f}", f"{lte['horizon_AU']:.0f}"
    t["lte_pc"], t["lte_nstar"] = f"{lte['dmax_pc_4GW_narrow_FAST']:.1f}", f"{lte['n_fast_stars_within_reach_4GW_narrow']:.0f}"
    # ---- posteriors
    P = K["posteriors"]
    t["prior_q95"] = num(P["i_beacon_1e17W"]["HI_line"]["base"]["prior_q95"])
    t["post17_hi"] = num(P["i_beacon_1e17W"]["HI_line"]["base"]["post_q95"])
    t["post17_110"] = num(P["i_beacon_1e17W"]["log_1_10GHz"]["base"]["post_q95"])
    t["postkii_110"] = num(P["i_beacon_KII"]["log_1_10GHz"]["base"]["post_q95"])
    # ---- Table 4
    d3 = dy[np.isclose(dy.T_DS_K, 300.0)].sort_values("covering_fraction")
    rows = []
    for _, r in d3.iterrows():
        k2 = r["floor_frac_300pc_if_100pct_complete_k2"]
        floor = f"{sci(k2)} (two survivors)" if np.isfinite(k2) else sci(r["floor_frac_300pc_if_100pct_complete_k0"])
        rows.append(f"| {r['covering_fraction']:.0%} | {sci(r['Suazo2022_published_frac_100pc'])} | {sci(r['Suazo2022_published_frac_300pc'])} | "
                    f"{sci(r['Suazo2022_published_frac_5kpc'])} | {floor} |")
    t["table4"] = "\n".join(rows)
    t["dy_sens_lo"], t["dy_sens_hi"] = f"{d3.N_sensitive_300pc_model.min() / 1e6:.1f}", f"{d3.N_sensitive_300pc_model.max() / 1e6:.1f}"
    # ---- robustness
    V = lambda frag, p: num(v17[name(frag)][p])
    t["v_lo_110"], t["v_hi_110"] = V("union upper bound", "log_1_10GHz"), V("union lower bound", "log_1_10GHz")
    t["v_raw_110"], t["v_ex2_110"], t["v_half_110"] = V("raw analytic", "log_1_10GHz"), V("omitting the two", "log_1_10GHz"), V("in-beam counts x0.5", "log_1_10GHz")
    t["v_eff_110"], t["v_x2_110"] = V("12.7%", "log_1_10GHz"), V("all thresholds x2", "log_1_10GHz")
    t["v_tb_hi"], t["v_tb_wide"] = V("tier B thresholds", "HI_line"), V("tier B thresholds", "log_0p1_100GHz")
    t["v_m10"], t["v_m30"], t["v_m60"] = V("thresholds x10 (drift", "HI_line"), V("thresholds x30 (drift", "HI_line"), V("thresholds x60 (drift", "HI_line")
    t["v_sah25"], t["v_scint_110"] = V("sky fraction 25%", "HI_line"), V("scintillation", "log_1_10GHz")
    dcr = K["drift_coverage_ratio"]["surface"]
    t["dr_hi"], t["dr_110"] = pct(dcr["HI_line"]), pct(dcr["log_1_10GHz"])
    gc = val[val.check.str.startswith("bl_gc_")]["ratio"]
    t["gc_ratio"] = f"{gc.min():.1f} to {gc.max():.1f}"
    # ---- limitations
    c = K["continuum_flux_density_Jy_KII"]
    t["nvss_near"], t["nvss_far"] = f"{c['NVSS_84MHz']['8.2kpc']:.0f}", f"{c['NVSS_84MHz']['25kpc']:.0f}"
    t["gb6_near"], t["gb6_far"] = f"{c['GB6_PMN_600MHz']['8.2kpc']:.0f}", f"{c['GB6_PMN_600MHz']['25kpc']:.0f}"
    t["lacki_kii"] = num(1e6 * 1e13 / L_KII * NMW)
    # ---- Grimaldi
    rows = []
    for _, r in gm[gm.record == "baseline"].iterrows():
        e_ = int(np.floor(np.log10(r["EIRP_W"]) + 1e-9)); m_ = r["EIRP_W"] / 10 ** e_
        rows.append(f"| {'' if abs(m_ - 1) < 1e-6 else f'{m_:.0f} × '}10^{e_} | {sci(r['GM2018_pi_prior'])} | {sci(r['this_work_coverage_fraction'])} | {r['ratio_this_over_GM']:.2g} |")
    t["table5"] = "\n".join(rows)
    c20 = K["coverage_fraction_1e20"]["base"]
    t["cov1e20_hi"], t["cov1e20_110"] = pct(c20["HI_line"]), pct(c20["log_1_10GHz"])
    return t


def main():
    t = build_tokens()
    used = set(re.findall(r"<<(\w+)>>", TEMPLATE))
    missing, unused = used - set(t), set(t) - used
    assert not missing, f"tokens without a value: {sorted(missing)}"
    assert not unused, f"values never used: {sorted(unused)}"
    text = re.sub(r"<<(\w+)>>", lambda m: t[m.group(1)], TEMPLATE)
    out = os.path.join(STUDY, "report.md")
    with open(out, "w", encoding="utf8") as f:
        f.write(text)
    body = text.split("## 9. Response to referee")[0]
    print(f"[make_report] wrote {out}: {len(body.split())} words before the response section, {len(text.split())} in total")
    return out


if __name__ == "__main__":
    main()