"""Step 4: injection-recovery of partial Dyson swarms -> completeness function C(f, T | class). For every parent star and every (f, T) the swarm model is added to the star's *actual* AllWISE/2MASS/Gaia measurements (so each star's real noise realisation, depth, photospheric level and quality flags are used), the photometric pipeline of step 3 is re-run on the modified photometry, and the star counts as recovered if it (i) would still be in the parent sample and (ii) exceeds the same detection threshold in a searchable band. Outputs data/completeness.csv : class, T, f, N_parent, N_recovered, completeness data/fmin_quantiles.csv : quantiles of the minimum detectable covering fraction RAW/derived/injection.npz : per-star recovery matrix and minimum detectable f """ import json import os import numpy as np import pandas as pd import photosphere as ph from common import RAW, DATA, T_GRID, F_GRID, jy2mag from swarm import swarm_response import importlib bs = importlib.import_module("02_build_sample") ex = importlib.import_module("03_excess") DERIVED = os.path.join(RAW, "derived") G_LIM = 20.0 # a dimmed host must stay brighter than this to be in the GCNS-based sample EKS_LIM = 0.25 # ... and keep a usable 2MASS Ks measurement (scaled uncertainty, mag) F_FINE = np.unique(np.round(np.concatenate([np.logspace(-5, np.log10(0.89), 40), F_GRID]), 9)) def load_cal(): c = json.load(open(os.path.join(DATA, "calibration.json"))) return {k: (np.array(v[0]), np.array(v[1])) for k, v in c.items()} def load_stars(): d = pd.read_csv(os.path.join(DERIVED, "stars.csv.gz"), dtype={"cc_flags": str, "moon_lev": str, "ph_qual": str, "designation": str}, low_memory=False) for c in ("has_aw", "q_base", "q3", "q4", "det3", "det4", "det"): d[c] = d[c].astype(bool) return d def inject(cal, s, f, resp, thr=ex.THR): """s: dict of numpy arrays for stars with usable anchors. Returns recovered mask and chi's.""" dim = 1.0 - f dmag = -2.5 * np.log10(dim) FKs = dim * s["FKs"] + f * s["fbol"] * resp["Ks"] FW1 = dim * s["FW1"] + f * s["fbol"] * resp["W1"] # W3/W4: remove the fraction f of the (predicted) photosphere, keep the star's own noise # realisation and any real excess, add the swarm F3 = s["FW3"] - f * s["P3K"] + f * s["fbol"] * resp["W3"] F4 = s["FW4"] - f * s["P4K"] + f * s["fbol"] * resp["W4"] # 2MASS uncertainty of the dimmed host: the larger of the star's own and the sample median # uncertainty at the new Ks magnitude eKs = np.maximum(s["e_Ks"], np.interp(jy2mag(FKs, "Ks"), *cal["eKs_curve"])) pred = ph.predict(cal, s["M_G"] + dmag, s["g_rp"], FKs, FW1) sig = ph.significance(cal, eKs, FW1, s["eW1"], F3, s["eW3"], F4, s["eW4"], pred) in_sample = ((s["dM"] + dmag < bs.DM_BELOW) & (s["M_G"] + dmag <= bs.MG_MAX) & (s["G"] + dmag <= G_LIM) & (eKs < EKS_LIM)) d3 = s["q3"] & (sig["chi3"] >= thr) d4 = s["q4"] & (sig["chi4"] >= thr) return in_sample & (d3 | d4), in_sample & d3, in_sample & d4 def star_arrays(d): use = (d.has_aw & np.isfinite(d.e_Ksmag) & np.isfinite(d.FW1) & (d.FW1 > 0) & np.isfinite(d.eW1)).values m = d[use] s = dict(FKs=m.FKs.values, FW1=m.FW1.values, eW1=m.eW1.values, FW3=m.FW3.values, eW3=m.eW3.values, FW4=m.FW4.values, eW4=m.eW4.values, P3K=m.P3K.values, P4K=m.P4K.values, e_Ks=m.e_Ksmag.values, M_G=m.M_G.values, g_rp=m.g_rp.values, dM=m.dM.values, G=m.Gmag.values, fbol=m.fbol.values, q3=m.q3.values, q4=m.q4.values) return use, s def main(): cal = load_cal() d = load_stars() use, s = star_arrays(d) n = len(d) rec = np.zeros((n, len(T_GRID), len(F_FINE)), dtype=bool) rec3 = np.zeros_like(rec) rec4 = np.zeros_like(rec) for it, T in enumerate(T_GRID): resp = swarm_response(T) for jf, f in enumerate(F_FINE): r, r3, r4 = inject(cal, s, f, resp) rec[use, it, jf] = r rec3[use, it, jf] = r3 rec4[use, it, jf] = r4 # minimum detectable covering fraction per star and temperature (first grid f recovered) fmin = np.full((n, len(T_GRID)), np.nan) for it in range(len(T_GRID)): anyrec = rec[:, it, :].any(axis=1) first = rec[:, it, :].argmax(axis=1) fmin[anyrec, it] = F_FINE[first[anyrec]] np.savez_compressed(os.path.join(DERIVED, "injection.npz"), sid=d.sid.values, rec=rec, rec3=rec3, rec4=rec4, fmin=fmin, f_fine=F_FINE, T=np.array(T_GRID)) rows = [] cls = d.cls.values searchable = (d.q3 | d.q4).values for c in ("A", "FGK", "M", "all"): sel = (cls == c) if c != "all" else np.ones(n, bool) for it, T in enumerate(T_GRID): for jf, f in enumerate(F_FINE): if not np.any(np.isclose(f, F_GRID, rtol=1e-6, atol=0)): continue nrec = int(rec[sel, it, jf].sum()) rows.append(dict(cls=c, T_K=T, f=f, N_parent=int(sel.sum()), N_searchable=int((sel & searchable).sum()), N_recovered=nrec, N_rec_W3=int(rec3[sel, it, jf].sum()), N_rec_W4=int(rec4[sel, it, jf].sum()), completeness=nrec / sel.sum(), completeness_searchable=nrec / max(1, (sel & searchable).sum()))) comp = pd.DataFrame(rows) comp.to_csv(os.path.join(DATA, "completeness.csv"), index=False, float_format="%.6g") qrows = [] for c in ("A", "FGK", "M", "all"): sel = ((cls == c) if c != "all" else np.ones(n, bool)) & searchable for it, T in enumerate(T_GRID): v = np.where(np.isfinite(fmin[sel, it]), fmin[sel, it], 2.0) # never recovered -> 2 (>1) q = np.quantile(v, [0.1, 0.25, 0.5, 0.75, 0.9]) qrows.append(dict(cls=c, T_K=T, n_searchable=int(sel.sum()), q10=q[0], q25=q[1], q50=q[2], q75=q[3], q90=q[4])) pd.DataFrame(qrows).to_csv(os.path.join(DATA, "fmin_quantiles.csv"), index=False, float_format="%.4g") piv = comp[comp.T_K == 300].pivot(index="f", columns="cls", values="completeness") print("Completeness (fraction of ALL parent stars) at T = 300 K:\n", piv.round(4).to_string()) for T in T_GRID: piv = comp[(comp.T_K == T) & comp.f.isin([1e-3, 1e-2, 0.1, 0.5, 0.9])].pivot(index="f", columns="cls", values="completeness") print(f"T={T} K\n", piv.round(4).to_string()) print(pd.DataFrame(qrows).to_string()) if __name__ == "__main__": main()