"""Step 11: ALFALFA HI sources without an optical counterpart ("dark" sources) as the f -> 1 limit. A disc whose starlight is ~100% intercepted would still be an HI source with a broad (rotation-dominated) profile but no optical counterpart. We count such sources in a.100 and characterise their environment. Revision: a closed accounting of the HI-selected parent sample. HI selection does not depend on starlight, so every broad-line HI source inside the SDSS footprint is either (a) usable for the Tully-Fisher test, (b) has good SDSS photometry but no unWISE stellar mass (i-band residual computed here), (c) has an optical counterpart with bad or no SDSS photometry (Durbala flags 2, 3), or (d) has no optical counterpart. Bins (b)-(d) are where a disc dimmed out of the photometric catalogue would have to hide; their size bounds that leakage without any model of the optical selection function. Hand-entered input (not regenerated): data/no_photometry_notes.csv, the by-eye class of the 35 W50 >= 200 km/s sources of bin (c). Outputs: data/dark_sources.csv, data/dark_summary.csv, data/no_photometry_sources.csv, data/hi_parent_accounting.csv ; JPG cutouts of the broad-line dark sources """ import os import numpy as np import pandas as pd from scipy.spatial import cKDTree from common import RAW, DATA, C_KMS, radec_to_xyz, chord_for_arcsec, poisson_ul, http_get VIEWER = "https://www.legacysurvey.org/viewer/" def main(): d = pd.read_csv(os.path.join(DATA, "alfalfa_all.csv"), low_memory=False) code1 = d[d.HI == 1] dark = code1[code1.RAO.isna()].copy() print(f" a.100 code-1 sources: {len(code1)}; without optical counterpart: {len(dark)} ({len(dark) / len(code1):.4f})") # environment: nearest a.100 source WITH an optical counterpart and overlapping velocity (tidal-debris test) oc = d[d.RAO.notna()] tree = cKDTree(radec_to_xyz(oc.RAO, oc.DEO)) xyz = radec_to_xyz(dark.RAJ2000, dark.DEJ2000) nb = tree.query_ball_point(xyz, chord_for_arcsec(1800.0)) sep_min, proj_kpc, mhi_nb = [], [], [] for i, lst in enumerate(nb): best = (np.nan, np.nan, np.nan) if lst: lst = np.asarray(lst) dv = np.abs(oc.Vhel.values[lst] - dark.Vhel.values[i]) ok = dv < 300 if ok.any(): sep = np.degrees(2 * np.arcsin(np.linalg.norm(tree.data[lst[ok]] - xyz[i], axis=1) / 2)) * 60 j = np.argmin(sep) best = (sep[j], sep[j] / 60 * np.pi / 180 * dark.Dist.values[i] * 1000, oc.logMHI.values[lst[ok]][j]) sep_min.append(best[0]); proj_kpc.append(best[1]); mhi_nb.append(best[2]) dark["nb_sep_arcmin"], dark["nb_proj_kpc"], dark["nb_logMHI"] = sep_min, proj_kpc, mhi_nb dark["has_neighbour_300kpc"] = dark.nb_proj_kpc < 300 cols = ["AGC", "RAJ2000", "DEJ2000", "Vhel", "W50", "e_W50", "SNR", "Dist", "logMHI", "Flag", "nb_sep_arcmin", "nb_proj_kpc", "nb_logMHI", "has_neighbour_300kpc"] dark[cols].to_csv(os.path.join(DATA, "dark_sources.csv"), index=False) rows = [] for snr in (6.5, 10.0): for wmin in (100, 150, 200, 300): par = code1[(code1.SNR >= snr) & (code1.Dist >= 40) & (code1.W50 >= wmin)] dk = dark[(dark.SNR >= snr) & (dark.Dist >= 40) & (dark.W50 >= wmin)] iso = dk[~dk.has_neighbour_300kpc] rows.append(dict(SNR_min=snr, W50_min=wmin, N_HI=len(par), N_dark=len(dk), frac_dark=len(dk) / len(par), UL95_frac_dark=poisson_ul(len(dk)) / len(par), N_dark_isolated=len(iso), UL95_frac_dark_isolated=poisson_ul(len(iso)) / len(par))) s = pd.DataFrame(rows) s.to_csv(os.path.join(DATA, "dark_summary.csv"), index=False) print(s.round(5).to_string(index=False)) print(f" dark sources with an a.100 optical galaxy within 300 kpc and 300 km/s: " f"{dark.has_neighbour_300kpc.sum()} of {len(dark)} ({dark.has_neighbour_300kpc.mean():.2f})") # cutouts for the broad-line, high-S/N dark sources (wide field, centred on the HI centroid) top = dark[(dark.SNR >= 10) & (dark.W50 >= 150) & (dark.Dist >= 40)] for r in top.itertuples(): for layer, tag in (("ls-dr9", "dark"), ("unwise-neo7", "darkW")): p = os.path.join(RAW, "cutouts", f"{tag}_{r.AGC}.jpg") if not (os.path.exists(p) and os.path.getsize(p) > 500): try: blob = http_get(f"{VIEWER}cutout.jpg?ra={r.RAJ2000:.5f}&dec={r.DEJ2000:.5f}&layer={layer}" f"&pixscale=1.0&size=300", timeout=120) open(p, "wb").write(blob) except Exception as e: print(f" WARN cutout AGC {r.AGC}: {e}") # ---- broad-line sources that HAVE an optical counterpart but no usable SDSS photometry (Durbala flag 2/3): # these never enter the Tully-Fisher sample, so they are a potential hiding place for very dim discs b = (d.HI == 1) & (d.SNR >= 10) & (d.Dist >= 40) & (d.W50 >= 200) hid = d[b & d.RAO.notna() & d.Flag.isin([2, 3])] n_in = int((b & d.Flag.isin([1, 2, 3])).sum()) print(f" broad-line HI sources inside the SDSS footprint: {n_in}; good photometry {int((b & (d.Flag == 1)).sum())}; " f"counterpart but bad/no SDSS photometry {len(hid)}; no counterpart {int((b & d.RAO.isna()).sum())}") hid = hid[["AGC", "Name", "RAO", "DEO", "Vhel", "W50", "SNR", "Dist", "logMHI", "Flag", "imag"]] notes_path = os.path.join(DATA, "no_photometry_notes.csv") notes = pd.read_csv(notes_path) if os.path.exists(notes_path) else pd.DataFrame(columns=["AGC", "vis_class", "note"]) hid = hid.merge(notes, on="AGC", how="left") hid.to_csv(os.path.join(DATA, "no_photometry_sources.csv"), index=False) print(f" by-eye class of these {len(hid)} sources: {hid.vis_class.fillna('not inspected').value_counts().to_dict()}") # ---- closed accounting of the HI-selected parent (inside the SDSS footprint: Durbala flags 1-3) fits = pd.read_csv(os.path.join(DATA, "alfalfa_tf_fits.csv")).set_index("relation").loc["i"] r_i = d.logLi - (fits.c0 + fits.c1 * (d.logV - fits.x0) + fits.c2 * (d.logV - fits.x0) ** 2) rows = [] for wmin in (120, 200): par = (d.HI == 1) & (d.SNR >= 10) & (d.Dist >= 40) & (d.W50 >= wmin) & d.Flag.isin([1, 2, 3]) usable = par & (d.Flag == 1) & d.logMsM.notna() & d.logLi.notna() & d.ba.notna() nomass = par & (d.Flag == 1) & ~usable testable = nomass & (d.ba <= 0.70) & (d.V >= 80) & r_i.notna() # same geometry cuts as the TF sample badphot = par & d.Flag.isin([2, 3]) & d.RAO.notna() nooc = par & d.RAO.isna() n_par = int(par.sum()) # (b): only objects that are untestable or lie >= 1.5 mag below the i-band relation can hide a dimmed disc nb_low = int((testable & (r_i < np.log10(0.25))).sum()) nb_untest = int((nomass & ~testable).sum()) n_leak = nb_low + int(badphot.sum()) + int(nooc.sum()) row = dict(W50_min=wmin, N_parent=n_par, N_usable=int(usable.sum()), N_good_phot_no_W1_mass=int(nomass.sum()), of_which_inclined_and_testable=int(testable.sum()), of_which_below_1p5mag_in_i=nb_low, of_which_face_on_or_untestable=nb_untest, N_counterpart_bad_or_no_phot=int(badphot.sum()), N_no_counterpart=int(nooc.sum()), N_possible_hidden=n_leak, frac_possible_hidden=n_leak / n_par, UL95_frac_possible_hidden=poisson_ul(n_leak) / n_par) if wmin == 200 and len(notes): unexpl = int((hid.vis_class != "ordinary").sum()) n2 = nb_low + unexpl + int(nooc.sum()) row.update(N_bad_phot_not_ordinary_by_eye=unexpl, N_possible_hidden_after_inspection=n2, UL95_frac_after_inspection=poisson_ul(n2) / n_par) rows.append(row) print(f" HI-selected parent (W50>={wmin}, S/N>=10, D>=40, in SDSS footprint): {n_par}; usable {int(usable.sum())}; " f"good photometry but no W1 mass {int(nomass.sum())} (testable {int(testable.sum())}, >=1.5 mag low in i: {nb_low}); " f"bad/no photometry {int(badphot.sum())}; no counterpart {int(nooc.sum())}; " f"possible hidden {n_leak} = {n_leak / n_par:.4f} (95% UL {poisson_ul(n_leak) / n_par:.4f})") pd.DataFrame(rows).to_csv(os.path.join(DATA, "hi_parent_accounting.csv"), index=False) # Does faintness push the SDSS axis ratio above the b/a <= 0.70 cut (which would hide a dimmed disc from the # TF test without putting it in any of the bins above)? Fraction passing the cut vs apparent magnitude. par = d[(d.HI == 1) & (d.SNR >= 10) & (d.Dist >= 40) & (d.W50 >= 120) & (d.Flag == 1) & d.logMsM.notna()] edges = [10, 14, 15, 16, 16.5, 17, 17.5, 18, 21] brow = [] for lo, hi in zip(edges[:-1], edges[1:]): b_ = par[(par.imag > lo) & (par.imag <= hi)] brow.append(dict(imag_lo=lo, imag_hi=hi, N=len(b_), frac_ba_le_0p7=(b_.ba <= 0.70).mean(), median_ba=b_.ba.median())) pd.DataFrame(brow).round(3).to_csv(os.path.join(DATA, "ba_vs_magnitude.csv"), index=False) print(" fraction of usable HI-selected sources passing b/a <= 0.70, by SDSS i magnitude:") print(pd.DataFrame(brow).round(3).to_string(index=False)) for r in hid.itertuples(): p = os.path.join(RAW, "cutouts", f"hid_{r.AGC}.jpg") if not (os.path.exists(p) and os.path.getsize(p) > 500): try: blob = http_get(f"{VIEWER}cutout.jpg?ra={r.RAO:.5f}&dec={r.DEO:.5f}&layer=ls-dr9&pixscale=1.0&size=200", timeout=120) open(p, "wb").write(blob) except Exception as e: print(f" WARN cutout AGC {r.AGC}: {e}") print(f" broad-line (W50>=150), SNR>=10, D>=40 dark sources: {len(top)}") print(top[cols].round(2).to_string(index=False)) if __name__ == "__main__": main()