#!/usr/bin/env python3 """Download every raw input used by the head-start study. All raw files go to RAW (default fermi_raw/headstart/, override with the HEADSTART_RAW environment variable). Files that already exist are kept, so the script is idempotent. Only the Python standard library is used, except for the Nataf+24 step, which reads a FITS table with astropy. Sources ------- * VizieR TAP : Casagrande+ 2011 GCS re-analysis (J/A+A/530/A138/catalog) Xiang & Rix 2022 subgiants (J/other/Nat/603.599/table1) and the smaller catalogues listed in VIZIER_TABLES * Gaia archive: DR3 gaia_source x astrophysical_parameters (FLAME), d < 200 pc DR3 FLAME for Hipparcos stars (cross-validation against GCS) * MIST v1.2 : non-rotating UBVRIplus isochrones (with 2MASS Ks) (subgiant selection weights) * Nataf+ 2024 : ages for ~400,000 subgiants (ApJ 976, 87). The catalogue is not on VizieR; the authors' repository (github.com/DavidMoiseNataf/Subgiants, README) links the file TotalSample.fits on Google Drive. The catalogue has no parallaxes, so Gaia DR3 astrometry and radial velocities are fetched by source_id for the stars near the solar annulus. """ import os import sys import time import urllib.parse import urllib.request RAW = os.environ.get("HEADSTART_RAW", "fermi_raw/headstart/") VIZIER_SYNC = "https://tapvizier.cds.unistra.fr/TAPVizieR/tap/sync" GAIA_TAP = "https://gea.esac.esa.int/tap-server/tap" MIST_URL = ("https://mist.science/data/tarballs_v1.2/" "MIST_v1.2_vvcrit0.0_UBVRIplus.txz") VIZIER_TABLES = { # local file name : (VizieR table, MAXREC) "gcs_casagrande2011.csv": ("J/A+A/530/A138/catalog", 20000), "xiangrix2022.csv": ("J/other/Nat/603.599/table1", 400000), } GAIA_200PC_QUERY = """ SELECT g.source_id, g.l, g.b, g.parallax, g.parallax_error, g.pmra, g.pmdec, g.ra, g.dec, g.radial_velocity, g.radial_velocity_error, g.phot_g_mean_mag, g.bp_rp, g.ruwe, a.age_flame, a.age_flame_lower, a.age_flame_upper, a.mass_flame, a.mass_flame_lower, a.mass_flame_upper, a.evolstage_flame, a.flags_flame, a.lum_flame, a.radius_flame, a.teff_gspphot, a.logg_gspphot, a.mh_gspphot, a.ag_gspphot, a.mh_gspspec, a.alphafe_gspspec, a.teff_gspspec, a.logg_gspspec, a.flags_gspspec FROM gaiadr3.gaia_source AS g JOIN gaiadr3.astrophysical_parameters AS a ON g.source_id = a.source_id WHERE g.parallax > 5 AND a.mass_flame BETWEEN 0.7 AND 1.3 """ GAIA_HIP_QUERY = """ SELECT h.original_ext_source_id AS hip, h.source_id, g.parallax, a.age_flame, a.age_flame_lower, a.age_flame_upper, a.mass_flame, a.evolstage_flame, a.flags_flame, a.mh_gspphot, a.teff_gspphot, a.logg_gspphot FROM gaiadr3.hipparcos2_best_neighbour AS h JOIN gaiadr3.gaia_source AS g ON g.source_id = h.source_id JOIN gaiadr3.astrophysical_parameters AS a ON a.source_id = h.source_id WHERE a.age_flame IS NOT NULL """ # Parent sample for the volume-complete subgiant census: every Gaia DR3 source # within 200 pc that is brighter than M_G = 4.8 with 0.55 < BP-RP < 1.25 # (turn-off stars, subgiants, lower giant branch), whether or not it has # spectroscopic parameters. GAIA_PARENT_QUERY = """ SELECT g.source_id, g.l, g.b, g.ra, g.dec, g.parallax, g.parallax_error, g.pmra, g.pmdec, g.radial_velocity, g.radial_velocity_error, g.phot_g_mean_mag, g.phot_bp_mean_mag, g.phot_rp_mean_mag, g.bp_rp, g.ruwe, a.teff_gspphot, a.ag_gspphot, a.ebpminrp_gspphot, a.teff_gspspec, a.logg_gspspec, a.mh_gspspec, a.alphafe_gspspec, a.flags_gspspec, a.age_flame, a.age_flame_lower, a.age_flame_upper, a.mass_flame, a.lum_flame, a.evolstage_flame FROM gaiadr3.gaia_source AS g LEFT OUTER JOIN gaiadr3.astrophysical_parameters AS a ON g.source_id = a.source_id WHERE g.parallax > 5 AND g.phot_g_mean_mag + 5 * LOG10(g.parallax / 100.0) < 4.8 AND g.bp_rp BETWEEN 0.55 AND 1.25 """ # StarHorse (Queiroz et al. 2023) results for Gaia DR3 RVS stars within 210 pc STARHORSE_QUERY = ( 'SELECT GaiaEDR3, GLON, GLAT, Mass16, Mass50, Mass84, Teff50, logg50, ' 'Met16, Met50, Met84, Dist50, AV50, Age05, Age16, Age50, Age84, Age95, ' 'AgeFlag, OutputFlag FROM "J/A+A/673/A155/gaiarvs" WHERE Dist50 < 0.21') SMALL_VIZIER = { "spina2018_solar_twins.csv": 'SELECT * FROM "J/MNRAS/474/2580/table4"', "nissen2020_solar_type.csv": 'SELECT * FROM "J/A+A/640/A81/table1"', "taniguchi2026_solar_twins.csv": ( 'SELECT GaiaDR3, Teff, logg, "[M/H]", BJ21rgeo, "Age(logg)", ' '"b_Age(logg)", "B_Age(logg)", "Age(GMAG)", "b_Age(GMAG)", ' '"B_Age(GMAG)", "Age(KsMAG)", "b_Age(KsMAG)", "B_Age(KsMAG)", ' '"Mini(GMAG)", Zmax FROM "J/A+A/707/A260/tablee1"'), } def log(msg): print(f"[fetch] {msg}", flush=True) def _get(url, data=None, timeout=900): req = urllib.request.Request(url, data=data) with urllib.request.urlopen(req, timeout=timeout) as r: return r.read() def vizier_table(fname, table, maxrec): path = os.path.join(RAW, fname) if os.path.exists(path) and os.path.getsize(path) > 1000: log(f"have {fname}") return path q = f'SELECT * FROM "{table}"' params = urllib.parse.urlencode({"REQUEST": "doQuery", "LANG": "ADQL", "FORMAT": "csv", "MAXREC": str(maxrec), "QUERY": q}) log(f"VizieR {table} -> {fname}") body = _get(VIZIER_SYNC + "?" + params) with open(path, "wb") as f: f.write(body) return path def vizier_query(fname, adql, maxrec=500000): """Arbitrary ADQL against VizieR TAP (sync).""" path = os.path.join(RAW, fname) if os.path.exists(path) and os.path.getsize(path) > 200: log(f"have {fname}") return path params = urllib.parse.urlencode({"REQUEST": "doQuery", "LANG": "ADQL", "FORMAT": "csv", "MAXREC": str(maxrec), "QUERY": adql}) log(f"VizieR query -> {fname}") body = _get(VIZIER_SYNC + "?" + params) with open(path, "wb") as f: f.write(body) return path def gaia_sync(adql, timeout=600, retries=3): """One synchronous Gaia-archive TAP query; returns CSV text.""" params = urllib.parse.urlencode({"REQUEST": "doQuery", "LANG": "ADQL", "FORMAT": "csv", "QUERY": " ".join(adql.split())}) last = None for attempt in range(retries): try: body = _get(GAIA_TAP + "/sync?" + params, timeout=timeout).decode() if body.lstrip().startswith("<"): raise RuntimeError("Gaia TAP error: " + body[:500]) return body except Exception as exc: # network hiccup: wait and retry last = exc time.sleep(10 * (attempt + 1)) raise RuntimeError(f"Gaia query failed after {retries} tries: {last}") # Parallax slices of equal volume out to 200 pc. The synchronous service # answers each slice in ~20-60 s (one big asynchronous job sat in the queue for # > 30 min when this was written). The service occasionally returns a # silently truncated table (exactly 8193 rows), so every slice is validated # against a COUNT(*) of the same selection and retried on mismatch. N_SLICES = 16 PLX_EDGES = [1.0e9] + [1000.0 / (200.0 * (k / N_SLICES) ** (1.0 / 3.0)) for k in range(1, N_SLICES + 1)] def _count_query(q): i = q.upper().index("FROM GAIADR3.GAIA_SOURCE") return "SELECT COUNT(*) AS n " + q[i:] def gaia_200pc(fname="gaia_flame_200pc.csv", base_query=None): base_query = " ".join((base_query or GAIA_200PC_QUERY).split()) path = os.path.join(RAW, fname) if os.path.exists(path) and os.path.getsize(path) > 1000: log(f"have {fname}") return path chunks = [] for hi, lo in zip(PLX_EDGES[:-1], PLX_EDGES[1:]): q = base_query.replace( "g.parallax > 5", f"g.parallax > {lo:.6f} AND g.parallax <= {hi:.6f}") assert q != base_query n_expected = int(gaia_sync(_count_query(q)).strip().splitlines()[-1]) for attempt in range(6): lines = gaia_sync(q).strip().splitlines() if len(lines) - 1 == n_expected: break log(f" truncated slice ({len(lines) - 1} of {n_expected}); retry") time.sleep(5) else: raise RuntimeError("Gaia slice kept coming back truncated") log(f"Gaia slice {lo:.3f} < plx <= {hi:.3f}: {n_expected} rows") chunks.append(lines if not chunks else lines[1:]) with open(path, "w") as f: f.write("\n".join("\n".join(c) for c in chunks) + "\n") return path def gaia_query(fname, adql): path = os.path.join(RAW, fname) if os.path.exists(path) and os.path.getsize(path) > 1000: log(f"have {fname}") return path body = gaia_sync(adql) with open(path, "w") as f: f.write(body) log(f"wrote {fname} ({len(body.splitlines()) - 1} rows)") return path def mist_isochrones(): d = os.path.join(RAW, "mist") os.makedirs(d, exist_ok=True) path = os.path.join(d, os.path.basename(MIST_URL)) if os.path.exists(path) and os.path.getsize(path) > 1.5e8: log("have MIST tarball") return path log("downloading MIST isochrones (157 MB)") with urllib.request.urlopen(MIST_URL, timeout=1800) as r, \ open(path, "wb") as f: while True: chunk = r.read(1 << 20) if not chunk: break f.write(chunk) return path # ------------------------------------------------------------ Nataf+ 2024 NATAF_DIR = os.path.join(RAW, "nataf2024") NATAF_FITS = os.path.join(NATAF_DIR, "TotalSample.fits") NATAF_GAIA = os.path.join(NATAF_DIR, "gaia_for_nataf.csv") # link given in the README of github.com/DavidMoiseNataf/Subgiants NATAF_URL = ("https://drive.usercontent.google.com/download?" "id=16s8GL_3smDSTqK1Lze4C5tixFSEYZvzq&export=download&confirm=t") NATAF_GAIA_COLS = ("source_id, ra, dec, l, b, parallax, parallax_error, pmra, " "pmdec, radial_velocity, radial_velocity_error, " "phot_g_mean_mag, bp_rp, ruwe") # pre-selection (catalogue-based distances, see nataf_photometric_distance): # the slab and its surroundings out to 1.3 kpc (shell test), plus the layers # above the slab out to |z| = 0.85 kpc inside d < 0.9 kpc (layer test) NATAF_PRE = dict(zmax_kpc=0.33, zmax_layers_kpc=0.85, dmax_layers_kpc=0.90, dR_kpc=0.60, R0=8.178, zsun_kpc=0.0208) def nataf_catalogue(): os.makedirs(NATAF_DIR, exist_ok=True) if os.path.exists(NATAF_FITS) and os.path.getsize(NATAF_FITS) > 7e7: log("have Nataf+24 TotalSample.fits") return NATAF_FITS log("downloading Nataf+24 TotalSample.fits (76 MB)") body = _get(NATAF_URL, timeout=1800) if not body.startswith(b"SIMPLE"): raise RuntimeError("Nataf+24 download did not return a FITS file") with open(NATAF_FITS, "wb") as f: f.write(body) return NATAF_FITS def nataf_photometric_distance(mass, logg, teff, gmag, av): """Distance [kpc] from the catalogue's own posterior medians: L = M g_sun/g (Teff/Teff_sun)^4, M_G = 4.67 - 2.5 log L, A_G = 0.84 A_V. Against Gaia parallaxes this is 2.3% short with 2% scatter (3,000 stars), hence the division by 0.977. Used only to pre-select which stars to look up in the Gaia archive; every cut uses the parallax.""" import numpy as np logL = np.log10(mass) - (logg - 4.438) + 4 * np.log10(teff / 5772.0) return 10 ** ((gmag - (4.67 - 2.5 * logL) - 0.84 * av + 5) / 5) / 1000 / 0.977 def nataf_preselected_ids(): import numpy as np from astropy.table import Table t = Table.read(nataf_catalogue()) d = nataf_photometric_distance(*(np.asarray(t[c], float) for c in ("mass", "logg", "Teff", "phot_g_mean_mag", "AV_model"))) l, b = np.radians(t["l2000"]), np.radians(t["b2000"]) z = np.abs(d * np.sin(b) + NATAF_PRE["zsun_kpc"]) R = np.hypot(NATAF_PRE["R0"] - d * np.cos(b) * np.cos(l), d * np.cos(b) * np.sin(l)) ann = np.abs(R - NATAF_PRE["R0"]) < NATAF_PRE["dR_kpc"] pre = ann & ((z < NATAF_PRE["zmax_kpc"]) | ((z < NATAF_PRE["zmax_layers_kpc"]) & (d < NATAF_PRE["dmax_layers_kpc"]))) return np.asarray(t["dr3_source_id"])[pre], len(t) def nataf_gaia(batch=4000): """Gaia DR3 astrometry + RVs for the pre-selected Nataf+24 stars (about 236,000; ~8 min). Incremental: only source_ids not yet in the local file are requested, so an interrupted download can be resumed.""" ids, n_all = nataf_preselected_ids() have = set() lines = [] if os.path.exists(NATAF_GAIA) and os.path.getsize(NATAF_GAIA) > 1000: with open(NATAF_GAIA) as f: lines = f.read().strip().splitlines() have = {int(ln.split(",")[0]) for ln in lines[1:]} todo = [i for i in ids.tolist() if i not in have] if not todo: log(f"have Gaia astrometry for {len(have)} Nataf+24 stars") return NATAF_GAIA log(f"Nataf+24: {len(ids)} of {n_all} stars pre-selected; " f"{len(todo)} to look up in the Gaia archive") for a in range(0, len(todo), batch): chunk = todo[a:a + batch] q = (f"SELECT {NATAF_GAIA_COLS} FROM gaiadr3.gaia_source " f"WHERE source_id IN ({','.join(map(str, chunk))})") data = urllib.parse.urlencode({"REQUEST": "doQuery", "LANG": "ADQL", "FORMAT": "csv", "QUERY": q}).encode() for attempt in range(5): try: rows = _get(GAIA_TAP + "/sync", data=data, timeout=600).decode().strip().splitlines() if len(rows) - 1 == len(chunk): break log(f" batch returned {len(rows) - 1} of {len(chunk)}; retry") except Exception as exc: log(f" batch failed ({exc}); retry") time.sleep(5 * (attempt + 1)) else: raise RuntimeError("Gaia look-up for Nataf+24 stars failed") lines.extend(rows if not lines else rows[1:]) with open(NATAF_GAIA, "w") as f: # checkpoint after each batch f.write("\n".join(lines) + "\n") log(f" Gaia look-up {min(a + batch, len(todo))}/{len(todo)}") return NATAF_GAIA def main(): os.makedirs(RAW, exist_ok=True) for fname, (table, maxrec) in VIZIER_TABLES.items(): vizier_table(fname, table, maxrec) gaia_200pc() gaia_200pc("gaia_parent_200pc.csv", GAIA_PARENT_QUERY) gaia_query("gaia_flame_hip.csv", GAIA_HIP_QUERY) vizier_query("starhorse_rvs_210pc.csv", STARHORSE_QUERY) for fname, q in SMALL_VIZIER.items(): vizier_query(fname, q) mist_isochrones() nataf_gaia() log("all raw inputs present in " + RAW) if __name__ == "__main__": sys.exit(main())