"""Download every raw catalogue used by the study (cached in RAW).
1. Gaia DR3 astrometric two-body orbits + gaia_source columns (Gaia archive TAP, async)
2. gaiadr3.binary_masses (Gaia Collaboration, Arenou et al. 2023)
3. AllWISE photometry of the NSS sources via the Gaia pre-computed cross-match
4. Shahaf et al. 2023 AMRF classification (VizieR J/MNRAS/518/2991 tables 1 and 2) and
Shahaf et al. 2024 Table A1 (Zenodo 14799795: colour excess of the non-class-I systems)
5. Mamajek "mean dwarf" sequence (mass, M_G, luminosity, colours)
6. Parent-population counts of Gaia DR3 stars within 300 pc in (distance, G, M_G, colour) cells
and a 1-in-40 random subsample of the same population
7. (best effort) arXiv abstracts of the papers behind the hand-compiled literature tables
8. WISE relative spectral response curves (IRSA) for band-averaged fluxes of cold blackbodies
9. AllWISE sources that are bright at 22 um and dark at 3.4 um (VizieR II/328/allwise) in 144
non-overlapping 5-degree cones at Galactic latitude |b| > 30 deg (added in the revision)
"""
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from common import RAW, GAIA_TAP, VIZIER_TAP, tap_async, tap_sync, http_get, fetch_to
def run(base, q, timeout=900):
"""Synchronous TAP first (fast on the Gaia archive, no row cap seen); async as fallback."""
try:
blob = tap_sync(base, q, timeout=timeout)
if blob[:200].lstrip().startswith(b"<"):
raise RuntimeError(blob[:500].decode(errors="replace"))
return blob
except Exception as exc: # noqa: BLE001
print(" sync failed (%s); trying async" % str(exc)[:200], flush=True)
return tap_async(base, q, timeout=5400)
ASTRO_TYPES = ("'Orbital','AstroSpectroSB1','OrbitalAlternative','OrbitalAlternativeValidated',"
"'OrbitalTargetedSearch','OrbitalTargetedSearchValidated'")
Q_NSS = f"""
SELECT n.source_id, n.nss_solution_type, n.parallax AS plx_nss, n.parallax_error AS plx_nss_err,
n.a_thiele_innes, n.b_thiele_innes, n.f_thiele_innes, n.g_thiele_innes,
n.a_thiele_innes_error, n.b_thiele_innes_error, n.f_thiele_innes_error, n.g_thiele_innes_error,
n.c_thiele_innes, n.h_thiele_innes,
n.period, n.period_error, n.eccentricity, n.eccentricity_error,
n.significance, n.goodness_of_fit, n.astrometric_n_good_obs_al, n.flags,
g.ra, g.dec, g.l, g.b, g.parallax AS plx_single, g.parallax_error AS plx_single_err,
g.phot_g_mean_mag, g.bp_rp, g.phot_bp_rp_excess_factor, g.ruwe,
g.visibility_periods_used, g.astrometric_excess_noise, g.radial_velocity,
g.ag_gspphot, g.ebpminrp_gspphot, g.teff_gspphot, g.logg_gspphot, g.mh_gspphot
FROM gaiadr3.nss_two_body_orbit AS n
JOIN gaiadr3.gaia_source AS g ON g.source_id = n.source_id
WHERE n.nss_solution_type IN ({ASTRO_TYPES})
"""
Q_BM = """
SELECT source_id, m1, m1_lower, m1_upper, m2, m2_lower, m2_upper,
fluxratio, fluxratio_lower, fluxratio_upper, combination_method, m1_ref, flag
FROM gaiadr3.binary_masses
"""
Q_WISE = f"""
SELECT n.source_id, x.angular_distance, x.number_of_neighbours, x.number_of_mates,
a.designation, a.w1mpro, a.w1mpro_error, a.w2mpro, a.w2mpro_error,
a.w3mpro, a.w3mpro_error, a.w4mpro, a.w4mpro_error, a.cc_flags, a.ext_flag, a.ph_qual
FROM gaiadr3.nss_two_body_orbit AS n
JOIN gaiadr3.allwise_best_neighbour AS x ON x.source_id = n.source_id
JOIN gaiadr1.allwise_original_valid AS a ON a.allwise_oid = x.allwise_oid
WHERE n.nss_solution_type IN ({ASTRO_TYPES})
"""
# Parent population: all Gaia DR3 sources with a 5-parameter-quality parallax within 300 pc,
# counted in cells of distance (25 pc), apparent G (0.5 mag), absolute G (0.25 mag) and
# BP-RP (0.1 mag). No extinction correction (d < 300 pc; discussed in the report).
Q_PARENT = """
SELECT dbin, gbin, mgbin, cbin, COUNT(*) AS n FROM (
SELECT FLOOR((1000.0/parallax)/25.0) AS dbin,
FLOOR(phot_g_mean_mag/0.5) AS gbin,
FLOOR((phot_g_mean_mag + 5.0*LOG10(parallax) - 10.0)/0.25) AS mgbin,
FLOOR(bp_rp/0.1) AS cbin
FROM gaiadr3.gaia_source
WHERE parallax > 3.3333333 AND parallax_over_error > 10
AND phot_g_mean_mag < 17.5 AND bp_rp IS NOT NULL
) AS t
GROUP BY dbin, gbin, mgbin, cbin
"""
# 1-in-40 random subsample of the same parent population (individual stars), used to calibrate
# the retention of the "single-light" cut and the colour-magnitude locus.
Q_PARENT_RANDOM = """
SELECT source_id, parallax, phot_g_mean_mag, bp_rp, ruwe
FROM gaiadr3.gaia_source
WHERE parallax > 3.3333333 AND parallax_over_error > 10
AND phot_g_mean_mag < 17.5 AND bp_rp IS NOT NULL AND MOD(random_index, 40) = 0
"""
# arXiv identifiers of the papers whose abstracts were used to verify the hand-compiled tables in
# 05_oort.py and 06_microlensing.py (abstracts are archived in RAW/arxiv for the referee).
ARXIV_IDS = """1507.08661 1507.08914 1909.04740 0911.2706 2403.12118 2009.02329 2311.00627 2009.08461
2206.03502 1310.3706 1601.01846 2412.03651 1809.02654 1909.04444 2401.17549 1709.00806 2607.18408
2302.01168 2303.08876 2107.06510 1909.11802 1912.00038 2504.06347 2410.09147 2508.18343 1609.06720
2610.07133 2110.07934 2409.02157 1703.06947 1704.01724 1905.08540 2310.13066 2112.08029 1901.07281
2112.01613 2201.13296 2202.01903 2308.03302 2609.28420 1508.06636 astro-ph/9906472 astro-ph/0207228
2608.24009 2512.07924 2106.05300 1906.02210 2410.23553 2005.14668 1509.05334 astro-ph/9812404
astro-ph/0405155 1309.0809 1704.05063 1604.01407 1808.05603 1811.07911 2012.11477 2005.12018
1608.07954 1706.04193 2411.00088 2411.00654 2209.00828 2309.15143 1809.01673 2201.11123 1804.08351
0811.2376 1508.02406 2507.19594 2401.14471 1208.5517 1208.5518 2405.00089 2409.11447""".split()
def fetch_arxiv_abstracts():
"""Best effort: archive the abstracts of the cited papers (not needed by the analysis code)."""
import json
import re
out = {}
for k in range(0, len(ARXIV_IDS), 25):
ids = ",".join(ARXIV_IDS[k:k + 25])
xml = http_get(f"https://export.arxiv.org/api/query?id_list={ids}&max_results=30").decode("utf-8", "replace")
for e in re.findall(r"(.*?)", xml, flags=re.S):
aid = re.search(r"(.*?)", e).group(1).split("/abs/")[-1]
out[aid] = {"title": re.sub(r"\s+", " ", re.search(r"(.*?)", e, flags=re.S).group(1)),
"first_author": re.search(r"(.*?)", e).group(1),
"published": re.search(r"(.*?)", e).group(1),
"abstract": re.sub(r"\s+", " ", re.search(r"(.*?)", e, flags=re.S).group(1))}
return json.dumps(out, indent=1).encode()
RSR_URL = "https://irsa.ipac.caltech.edu/data/WISE/docs/release/Prelim/expsup/figures/RSR-{b}.txt"
# ---- AllWISE sources bright at 22 um and dark at 3.4 um, in sample fields ---------------------------
W4DARK_RADIUS = 5.0 # deg
W4DARK_RINGS = ((36.0, 24), (48.0, 20), (60.0, 15), (72.0, 10), (84.0, 3)) # |b|, cones per ring
def w4dark_fields():
"""Centres of 2 x 72 non-overlapping cones of radius 5 deg, all at |b| >= 31 deg."""
from astropy.coordinates import SkyCoord
import astropy.units as u
rows = []
for sgn in (1, -1):
for b, n in W4DARK_RINGS:
for k in range(n):
l = (360.0 / n) * (k + (0.5 if sgn < 0 else 0.0))
c = SkyCoord(l=l * u.deg, b=sgn * b * u.deg, frame="galactic").icrs
rows.append((len(rows), l, sgn * b, float(c.ra.deg), float(c.dec.deg)))
return rows
def w4dark_query(ra, dec):
return ("SELECT AllWISE, RAJ2000, DEJ2000, W1mag, e_W1mag, snr1, W2mag, e_W2mag, snr2, "
"W3mag, e_W3mag, snr3, W4mag, e_W4mag, snr4, chi2W4, ccf, ex, var, qph, nb, na, covW4, "
"Jmag, pmRA, e_pmRA, pmDE, e_pmDE "
'FROM "II/328/allwise" WHERE 1=CONTAINS(POINT(\'ICRS\',RAJ2000,DEJ2000),'
f"CIRCLE('ICRS',{ra:.5f},{dec:.5f},{W4DARK_RADIUS})) "
"AND W4mag<7.7 AND snr4>=5 AND (W1mag IS NULL OR W1mag-W4mag>8.0)")
def fetch_w4dark():
"""One VizieR cone query per field (cached per field), concatenated with a field column."""
import io
import time
import pandas as pd
d = os.path.join(RAW, "w4dark")
os.makedirs(d, exist_ok=True)
frames, meta = [], []
for fid, l, b, ra, dec in w4dark_fields():
path = os.path.join(d, f"field_{fid:03d}.csv")
if not (os.path.exists(path) and os.path.getsize(path) > 0):
for attempt in range(5):
try:
blob = tap_sync(VIZIER_TAP, w4dark_query(ra, dec), timeout=280)
if blob[:200].lstrip().startswith(b"<"):
raise RuntimeError(blob[:300].decode(errors="replace"))
break
except Exception as exc: # noqa: BLE001
print(f" field {fid}: retry after {str(exc)[:120]}", flush=True)
time.sleep(5 + 10 * attempt)
else:
raise RuntimeError(f"field {fid} failed")
with open(path + ".part", "wb") as f:
f.write(blob)
os.replace(path + ".part", path)
df = pd.read_csv(path)
df.insert(0, "field", fid)
frames.append(df)
meta.append({"field": fid, "l": l, "b": b, "ra": ra, "dec": dec, "radius_deg": W4DARK_RADIUS, "n_rows": len(df)})
if fid % 12 == 0:
print(f" field {fid}: {len(df)} rows", flush=True)
pd.DataFrame(meta).to_csv(os.path.join(RAW, "w4dark_fields.csv"), index=False)
return pd.concat(frames, ignore_index=True).to_csv(index=False).encode()
def counted(q, what="*"):
"""Run a (small) Gaia query and verify the number of rows against a server-side COUNT(*)."""
import io
import pandas as pd
from_part = q[q.upper().index("FROM"):]
for attempt in range(4):
n_expected = int(run(GAIA_TAP, "SELECT COUNT(*) AS n " + from_part).decode().split()[1])
blob = run(GAIA_TAP, q, timeout=1700)
n = len(pd.read_csv(io.BytesIO(blob)))
print(f" {n} rows, server count {n_expected}", flush=True)
if n == n_expected:
return blob
raise RuntimeError("row count mismatch after retries")
def chunked(q, alias="n", nchunk=12):
"""Run a Gaia-archive query in HEALPix level-0 chunks of source_id and verify the row count
of every chunk against a server-side COUNT(*) (one un-chunked sync download silently lost
3 per cent of the rows)."""
import io
import re
import pandas as pd
step = 2 ** 59
frames = []
from_part = q[q.upper().index("FROM"):]
for k in range(nchunk):
cond = f" AND {alias}.source_id >= {k*step} AND {alias}.source_id < {(k+1)*step}"
for attempt in range(4):
n_expected = int(run(GAIA_TAP, "SELECT COUNT(*) AS n " + from_part + cond).decode().split()[1])
df = pd.read_csv(io.BytesIO(run(GAIA_TAP, q + cond)))
print(f" chunk {k}: {len(df)} rows, server count {n_expected}", flush=True)
if len(df) == n_expected:
break
else:
raise RuntimeError(f"chunk {k}: row count mismatch after retries")
frames.append(df)
out = pd.concat(frames, ignore_index=True)
return out.to_csv(index=False).encode()
def main():
print("[download] Gaia NSS astrometric orbits")
fetch_to(os.path.join(RAW, "nss_astrometric.csv"), lambda: chunked(Q_NSS))
print("[download] gaiadr3.binary_masses")
fetch_to(os.path.join(RAW, "binary_masses.csv"), lambda: run(GAIA_TAP, Q_BM))
print("[download] AllWISE photometry for NSS sources")
fetch_to(os.path.join(RAW, "nss_allwise.csv"), lambda: chunked(Q_WISE))
print("[download] Shahaf+2023 tables")
for t in ("table1", "table2"):
q = f'SELECT * FROM "J/MNRAS/518/2991/{t}"'
fetch_to(os.path.join(RAW, f"shahaf2023_{t}.csv"), lambda q=q: run(VIZIER_TAP, q))
print("[download] Shahaf+2024 Table A1 (Zenodo record 14799795)")
fetch_to(os.path.join(RAW, "shahaf2024_tableA1.fits"), lambda: http_get(
"https://zenodo.org/api/records/14799795/files/table_I_20250126.fits/content"))
print("[download] Mamajek dwarf sequence")
fetch_to(os.path.join(RAW, "mamajek_dwarf.txt"), lambda: http_get(
"https://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt"))
print("[download] Gaia parent-population cell counts (d < 300 pc)")
fetch_to(os.path.join(RAW, "parent_cells_300pc.csv"), lambda: run(GAIA_TAP, Q_PARENT, timeout=1700))
print("[download] Gaia parent-population 1-in-40 random subsample (d < 300 pc)")
fetch_to(os.path.join(RAW, "parent_random_1in40.csv"), lambda: counted(Q_PARENT_RANDOM))
print("[download] arXiv abstracts of the cited papers (best effort)")
os.makedirs(os.path.join(RAW, "arxiv"), exist_ok=True)
try:
fetch_to(os.path.join(RAW, "arxiv", "arxiv_abstracts_cited_v2.json"), fetch_arxiv_abstracts)
except Exception as exc: # noqa: BLE001 (the analysis does not depend on these files)
print(" arXiv abstracts not fetched:", str(exc)[:200])
print("[download] WISE relative spectral response curves")
for b in ("W1", "W2", "W3", "W4"):
fetch_to(os.path.join(RAW, f"RSR-{b}.txt"), lambda b=b: http_get(RSR_URL.format(b=b)))
print("[download] AllWISE 22-um-bright, 3.4-um-dark sources in 144 high-latitude cones (VizieR)")
fetch_to(os.path.join(RAW, "w4dark_sources.csv"), fetch_w4dark)
print("[download] done")
if __name__ == "__main__":
main()