"""Step 2: build the ALFALFA-SDSS Tully-Fisher parent sample with automated quality flags. Inputs : a100.csv (Haynes+2018), durbala_t1/t2.csv (Durbala+2020), tempel_dr12.csv, SDSS DR16 Petrosian radii of the gold galaxies (VizieR V/154, fetched by objID and cached) Output : data/alfalfa_all.csv (every a.100 source with derived quantities + flags) """ import io import os import urllib.parse import numpy as np import pandas as pd from scipy.spatial import cKDTree from common import RAW, DATA, C_KMS, VIZIER_TAP, http_get, radec_to_xyz, chord_for_arcsec Q0 = 0.20 # intrinsic axial ratio of an edge-on disc TURB = 6.5 # km/s, linear turbulence correction (Springob et al. 2005) MSUN_I_AB = 4.53 # absolute AB magnitude of the Sun in SDSS i BEAM_R_ARCSEC = 120 # neighbour search radius: ALFALFA beam FWHM is ~3.5' def sdss_sizes(objids, chunk=1000): """i-band Petrosian radius (arcsec) from SDSS DR16 for a list of objIDs; cached in RAW. Durbala et al. give DR15 objIDs; a minority are not found under the same ID in DR16 (they come back without a radius and are treated as 'size unknown' downstream).""" path = os.path.join(RAW, "sdss16_sizes.csv") have = pd.read_csv(path) if os.path.exists(path) else pd.DataFrame(columns=["objID", "iPrad", "rPrad"]) todo = np.array(sorted(set(int(x) for x in objids) - set(have.objID.astype("int64"))), dtype="int64") new = [] for k in range(0, len(todo), chunk): ids = todo[k:k + chunk] q = ('SELECT objID, "iPrad", "rPrad" FROM "V/154/sdss16" WHERE objID IN (%s)' % ",".join(map(str, ids))) blob = http_get(VIZIER_TAP, data=urllib.parse.urlencode( {"REQUEST": "doQuery", "LANG": "ADQL", "FORMAT": "csv", "QUERY": q}).encode(), timeout=600) if blob[:200].lstrip().startswith(b"<"): raise RuntimeError(f"SDSS size query failed: {blob[:400]!r}") got = pd.read_csv(io.BytesIO(blob)).drop_duplicates("objID") got = pd.DataFrame({"objID": ids}).merge(got, on="objID", how="left") # keep misses, so they are not re-queried new.append(got) print(f" SDSS DR16 sizes: chunk {k // chunk + 1}/{(len(todo) - 1) // chunk + 1}, " f"{got.iPrad.notna().sum()}/{len(ids)} found", flush=True) if new: have = pd.concat([have] + new, ignore_index=True) have.to_csv(path, index=False) return have def build(q0=Q0, with_sizes=True): a = pd.read_csv(os.path.join(RAW, "a100.csv")) t1 = pd.read_csv(os.path.join(RAW, "durbala_t1.csv"), dtype={"ObjID": "Int64"}) # 19-digit IDs: keep exact t2 = pd.read_csv(os.path.join(RAW, "durbala_t2.csv")) d = a.drop(columns=["recno"]).merge( t1.drop(columns=["recno", "Dist", "e_Dist", "Dprop", "RVel"]).rename( columns={"RAJ2000": "RAd", "DEJ2000": "DEd"}), on="AGC", how="left") d = d.merge(t2.drop(columns=["recno", "logMHI", "e_logMHI"]), on="AGC", how="left") d = d.rename(columns={"b/a": "ba", "g-i": "g_i", "e_b/a": "e_ba", "e_g-i": "e_g_i"}) # --- kinematics z = d.Vhel / C_KMS d["Wc"] = d.W50 / (1 + z) - TURB # W50 is already corrected for instrumental broadening cosi2 = np.clip((d.ba ** 2 - q0 ** 2) / (1 - q0 ** 2), 0, 1) d["sini"] = np.sqrt(1 - cosi2) d["V"] = d.Wc / (2 * d.sini) # rotation velocity, km/s d["logV"] = np.log10(d.V.where(d.V > 0)) d["logVmin"] = np.log10((d.Wc / 2).where(d.Wc > 0)) # edge-on assumption: strict lower limit on V # --- photometry (all on the a.100 distance scale, H0 = 70) d["logLi"] = -0.4 * (d.iMAG - MSUN_I_AB) # i-band luminosity, Galactic+internal extinction corrected d["logLi_noint"] = -0.4 * (d.iMAG + d.Ai - MSUN_I_AB) # same without the internal-extinction correction d["logLW1"] = d.logMsM # W1 luminosity x 0.45 Msun/Lsun (McGaugh & Schombert 2015) d["logMbar"] = np.log10(10 ** d.logMsM + 1.33 * 10 ** d.logMHI) d["logMbarT"] = np.log10(10 ** d.logMsT + 1.33 * 10 ** d.logMHI) d["fgas"] = 1.33 * 10 ** d.logMHI / 10 ** d.logMbar # --- HI centroid to optical counterpart offset has_oc = d.RAO.notna() cosd = np.cos(np.radians(d.DEJ2000)) d["hi_opt_off"] = np.hypot((d.RAJ2000 - d.RAO) * cosd, d.DEJ2000 - d.DEO) * 3600.0 # --- confusion flag 1: SDSS DR12 spectroscopic neighbours inside the beam tp = pd.read_csv(os.path.join(RAW, "tempel_dr12.csv")) ttree = cKDTree(radec_to_xyz(tp.RAJ2000, tp.DEJ2000)) tcz = tp.zobs.values * C_KMS ra = np.where(has_oc, d.RAO, d.RAJ2000) de = np.where(has_oc, d.DEO, d.DEJ2000) xyz = radec_to_xyz(ra, de) nb = ttree.query_ball_point(xyz, chord_for_arcsec(BEAM_R_ARCSEC)) r_self = chord_for_arcsec(20.0) n_conf = np.zeros(len(d), int) nearest = np.full(len(d), np.nan) self_match = np.zeros(len(d), bool) cz = d.Vhel.values w = d.W50.values for i, lst in enumerate(nb): if not lst: continue lst = np.asarray(lst) sep = np.linalg.norm(ttree.data[lst] - xyz[i], axis=1) same = np.abs(tcz[lst] - cz[i]) < (w[i] / 2 + 200.0) self_match[i] = np.any(same & (sep < r_self)) comp = same & (sep >= r_self) n_conf[i] = comp.sum() if comp.any(): nearest[i] = np.degrees(2 * np.arcsin(sep[comp].min() / 2)) * 3600 d["n_spec_nb"] = n_conf # same-redshift SDSS spectroscopic galaxies 20"-120" away d["nb_sep"] = nearest d["oc_has_sdss_spec"] = self_match # crude footprint flag for the spectroscopic sample n1deg = np.array([len(x) for x in ttree.query_ball_point(xyz, chord_for_arcsec(3600.0))]) d["in_spec_footprint"] = n1deg >= 3 # --- confusion flag 2: another a.100 detection within 4' and overlapping in velocity atree = cKDTree(radec_to_xyz(d.RAJ2000, d.DEJ2000)) pairs = atree.query_pairs(chord_for_arcsec(240.0), output_type="ndarray") hi_nb = np.zeros(len(d), int) for i, j in pairs: if abs(cz[i] - cz[j]) < (w[i] + w[j]) / 2 + 100.0: hi_nb[i] += 1 hi_nb[j] += 1 d["n_hi_nb"] = hi_nb # --- sample definitions ------------------------------------------------ base = (d.HI == 1) & (d.Flag == 1) & d.logMsM.notna() & d.logLi.notna() & d.ba.notna() geom = (d.ba <= 0.70) & (d.Dist >= 40.0) & (d.V >= 80.0) d["gold"] = base & geom & (d.SNR >= 10.0) d["silver"] = base & geom & (d.SNR >= 6.5) & (d.SNR < 10.0) # automated "clean" flag: no known same-redshift neighbour in the beam, HI centroid on the optical galaxy d["auto_clean"] = (d.n_spec_nb == 0) & (d.n_hi_nb == 0) & (d.hi_opt_off <= 40.0) # --- mean i-band surface brightness inside the Petrosian ellipse (gold galaxies only; mag/arcsec^2). # Only used to ask whether a dimmed galaxy would still lie inside the range of surface brightness that # the catalogue actually contains (s07); a crude but uniform proxy. d["iPrad"] = np.nan if with_sizes: g = (d.gold & d.ObjID.notna()).values ids = d.ObjID[g].astype("int64").values sz = sdss_sizes(ids).drop_duplicates("objID").set_index("objID").iPrad.astype(float) d.loc[g, "iPrad"] = sz.reindex(ids).values okr = d.iPrad.between(3.0, 200.0) # 2.97" is the pipeline's failure value; >200" is a bright-star halo d["mu_i"] = np.where(okr, d.imag + 2.5 * np.log10(np.pi * d.iPrad ** 2 * d.ba), np.nan) return d def main(): d = build() out = os.path.join(DATA, "alfalfa_all.csv") d.to_csv(out, index=False) print(f" a.100 sources: {len(d)}") steps = [("HI code 1", d.HI == 1), ("+ Durbala Flag 1 (good SDSS photometry)", d.Flag == 1), ("+ unWISE W1 stellar mass available", d.logMsM.notna()), ("+ b/a <= 0.70", d.ba <= 0.70), ("+ Dist >= 40 Mpc", d.Dist >= 40), ("+ V_rot >= 80 km/s", d.V >= 80), ("+ SNR >= 10 [GOLD]", d.SNR >= 10)] m = np.ones(len(d), bool) rows = [] for name, cut in steps: m &= cut.values rows.append((name, int(m.sum()))) print(f" {name:45s} {m.sum():6d}") pd.DataFrame(rows, columns=["cut", "N"]).to_csv(os.path.join(DATA, "alfalfa_cut_flow.csv"), index=False) g = d[d.gold] print(f" gold: {len(g)}; V>=100: {(g.V >= 100).sum()}; auto_clean: {g.auto_clean.sum()} " f"(spec nb {(g.n_spec_nb > 0).sum()}, HI nb {(g.n_hi_nb > 0).sum()}, offset>40\" {(g.hi_opt_off > 40).sum()})") print(f" gold in SDSS spectroscopic footprint: {g.in_spec_footprint.sum()}") print(f" silver (6.5<=SNR<10): {d.silver.sum()}") print(f" gold with a usable SDSS DR16 Petrosian radius (surface brightness known): {g.mu_i.notna().sum()} of {len(g)}; " f"mu_i percentiles 50/98/99.5: {np.nanpercentile(g.mu_i, [50, 98, 99.5]).round(2)}") if __name__ == "__main__": main()