"""Part 1: heat-producing-element (HPE) budgets of Sun-like stars versus formation epoch. MEASURED quantities (differential photospheric abundances relative to the Sun): Th : Botelho et al. 2019 (solar twins, Th II 4019 A) -> direct 232Th Eu : Spina et al. 2018 (solar twins), Delgado Mena et al. 2017 (HARPS-GTO), Battistini & Bensby 2016 (F/G dwarfs) -> r-process proxy for Th and U K : GALAH DR3 and DR4 (K I 7699 A, non-LTE) -> stable K (39K) Mg, Si (rock-forming reference), Fe, isochrone ages from the same sources. MEASURED versus MODEL in the age trend (revision): the europium-per-rock ratio [Eu/rock] is data. The conversion of Eu to INITIAL Th adds a one-zone free-decay factor (`ism_factor`), which is a model and by itself contributes about +0.010 dex/Gyr. The only empirical test of that factor, [Th/Eu] at birth versus age in the solar twins, gives about 0.4 of the predicted gradient. Both the raw fits (`fits_eu_raw`) and the fits with the full factor (`fits_eu`) are tabulated; the population model uses the raw fits plus a strength parameter for the decay factors (population.py). INFERRED quantities: U : never measured in these stars; assumed co-produced with Th (and Eu) in the r-process. 40K/K, 235U/238U and the free-decay corrections: a one-zone, uniform-production model of the interstellar medium (see `ism_factor`), not data. planet/star ratio: Th/(Mg,Si) and U/(Mg,Si) of the planet's silicate part are assumed equal to the star's (refractory lithophiles); K is moderately volatile and its planetary depletion is a prior, not a measurement. Definition. For each star, log10 f_X = [X/R] + corrections, where R is the mean of Mg and Si (the two dominant mantle cations) and f_X is the INITIAL concentration of X in the silicate part of a rocky planet relative to the Earth's initial concentration. By construction the Sun/Earth has f = 1. """ import os import numpy as np import pandas as pd from scipy.optimize import minimize_scalar from common import RAW_DIR, DATA_DIR, HALF_LIFE, AGE_SUN, T_UNIVERSE, LN2, SEED T_NUC_START = 0.5 # Gyr of cosmic time before nucleosynthesis starts (assumption) DEX_PER_GYR_TH = np.log10(2.0) / HALF_LIFE["Th232"] # 0.0215 dex/Gyr # --------------------------------------------------------------------------------------------- # one-zone interstellar-medium decay model (assumption, used only for nuclides not measured) # --------------------------------------------------------------------------------------------- def ism_factor(nuclide, age_gyr, age_ref=AGE_SUN): """log10 of [N_rad/N_stable](birth of star of age `age_gyr`) / same(birth of the Sun), for a radionuclide produced at a constant rate per unit stable co-produced yield since T_NUC_START. N_rad/N_stable at Galactic nucleosynthesis age T is (1 - exp(-lam T)) / (lam T).""" lam = LN2 / HALF_LIFE[nuclide] def g(age): T = np.clip(T_UNIVERSE - T_NUC_START - np.asarray(age, dtype=float), 0.05, None) return (1.0 - np.exp(-lam * T)) / (lam * T) return np.log10(g(age_gyr) / g(age_ref)) # --------------------------------------------------------------------------------------------- # sample construction # --------------------------------------------------------------------------------------------- def _hip(name): digits = "".join(ch for ch in str(name) if ch.isdigit()) return int(digits) if digits else -1 def load_solar_twins(): """Spina+2018 / Bedell+2018 solar twins (79) with Eu, Mg, Si, ages; Th from Botelho+2019.""" s2 = pd.read_csv(os.path.join(RAW_DIR, "spina2018_table2.csv")) s3 = pd.read_csv(os.path.join(RAW_DIR, "spina2018_table3.csv")) s4 = pd.read_csv(os.path.join(RAW_DIR, "spina2018_table4.csv")) be = pd.read_csv(os.path.join(RAW_DIR, "bedell2018_table2.csv")) be["HIP"] = be["Name"].map(_hip) b1 = pd.read_csv(os.path.join(RAW_DIR, "botelho2019_table1.csv")) b3 = pd.read_csv(os.path.join(RAW_DIR, "botelho2019_table3.csv")) b1["HIP"] = b1["Name"].map(_hip) b3["HIP"] = b3["Name"].map(_hip) t = s2[["HIP", "Teff", "logg", "[Fe/H]", "e_[Fe/H]"]].merge( s3[["HIP", "[EuII/H]", "e_[EuII/H]"]], on="HIP").merge( s4[["HIP", "Agemean", "Agestd", "AgeLL1s", "AgeUL1s"]], on="HIP").merge( be[["HIP", "[MgI/H]", "e_[MgI/H]", "[SiI/H]", "e_[SiI/H]"]], on="HIP", how="left") t = t.rename(columns={"[Fe/H]": "feh", "e_[Fe/H]": "e_feh", "[EuII/H]": "eu_h", "e_[EuII/H]": "e_eu_h", "Agemean": "age", "[MgI/H]": "mg_h", "e_[MgI/H]": "e_mg_h", "[SiI/H]": "si_h", "e_[SiI/H]": "e_si_h"}) t["e_age"] = np.maximum(0.5 * (t["AgeUL1s"] - t["AgeLL1s"]), t["Agestd"]) b3["th_published"] = b3["Prev"].astype(str).str.strip() == "*" # the 58 stars of the printed table 3 th = b3[b3["HIP"] > 0][["HIP", "[Th/H]", "e_[Th/H]", "SNR", "th_published"]].rename( columns={"[Th/H]": "th_h", "e_[Th/H]": "e_th_h", "SNR": "snr_th"}) t = t.merge(th, on="HIP", how="left") # Botelho's own ages for the Th stars (same isochrone method) are used as a cross-check only t = t.merge(b1[b1["HIP"] > 0][["HIP", "Age"]].rename(columns={"Age": "age_botelho"}), on="HIP", how="left") # Exclusions of Botelho et al. (their sect. 4): four alpha-rich old twins and one s-rich star excl = {28066, 30476, 73241, 74432, 64150} t["excluded_botelho"] = t["HIP"].isin(excl) t["sample"] = "solar_twins" return t def load_delgado_mena(): d1 = pd.read_csv(os.path.join(RAW_DIR, "delgadomena2017_table1.csv")) d3 = pd.read_csv(os.path.join(RAW_DIR, "delgadomena2017_table3.csv")) d2 = pd.read_csv(os.path.join(RAW_DIR, "delgadomena2019_table2.csv")) for d in (d1, d2, d3): d["Star"] = d["Star"].str.strip() t = d1[["Star", "Teff", "logg", "[Fe/H]", "e_[Fe/H]"]].merge( d3[["Star", "S/N", "pop", "[EuII/Fe]", "e_[EuII/Fe]", "[MgI/Fe]", "e_[MgI/Fe]", "[SiI/Fe]", "e_[SiI/Fe]"]], on="Star").merge(d2[["Star", "Mass", "Age", "e_Age"]], on="Star") t = t.rename(columns={"[Fe/H]": "feh", "e_[Fe/H]": "e_feh", "Age": "age", "e_Age": "e_age"}) t["eu_h"] = t["[EuII/Fe]"] + t["feh"] t["e_eu_h"] = t["e_[EuII/Fe]"] t["mg_h"] = t["[MgI/Fe]"] + t["feh"] t["si_h"] = t["[SiI/Fe]"] + t["feh"] t["e_mg_h"] = t["e_[MgI/Fe]"] t["e_si_h"] = t["e_[SiI/Fe]"] t["pop"] = t["pop"].str.strip() t["sample"] = "harps_gto" return t def load_battistini(): b = pd.read_csv(os.path.join(RAW_DIR, "battistini2016_abundances.csv")) c = pd.read_csv(os.path.join(RAW_DIR, "bensby2014_tablec3.csv")) t = b[["HIP", "[Eu/H]"]].merge( c[["HIP", "Teff", "logg", "Age", "b_Age", "B_Age", "Fe/H", "e_Fe/H", "Mg/Fe", "e_Mg/Fe", "Si/Fe", "e_Si/Fe", "td/d"]], on="HIP") t = t.rename(columns={"[Eu/H]": "eu_h", "Fe/H": "feh", "e_Fe/H": "e_feh", "Age": "age"}) t["e_age"] = 0.5 * (t["B_Age"] - t["b_Age"]) t["mg_h"] = t["Mg/Fe"] + t["feh"] t["si_h"] = t["Si/Fe"] + t["feh"] t["e_mg_h"] = t["e_Mg/Fe"] t["e_si_h"] = t["e_Si/Fe"] t["e_eu_h"] = 0.08 # Battistini & Bensby quote typical random errors of 0.05-0.1 dex for Eu t["sample"] = "battistini_bensby" return t def load_galah(release): if release == 3: g = pd.read_csv(os.path.join(RAW_DIR, "galah_dr3_sunlike.csv")) g = g.rename(columns={"snrc2iraf": "snr", "Age": "age", "e_Age": "e_age"}) else: g = pd.read_csv(os.path.join(RAW_DIR, "galah_dr4_sunlike.csv")) g = g.rename(columns={"SNRpxCCD2": "snr", "Age": "age"}) g["e_age"] = np.nan g = g.rename(columns={"[Fe/H]": "feh", "e_[Fe/H]": "e_feh"}) for el in ("Mg", "Si", "K", "Eu"): g[f"{el.lower()}_h"] = g[f"[{el}/Fe]"] + g["feh"] g[f"e_{el.lower()}_h"] = g[f"e_[{el}/Fe]"] g[f"flag_{el.lower()}"] = g[f"f_[{el}/Fe]"] g["sample"] = f"galah_dr{release}" return g def add_rock_ratios(t, ref="MgSi"): """[X/R] with R = mean of Mg and Si (default), Mg only, or Si only.""" if ref == "MgSi": r = 0.5 * (t["mg_h"] + t["si_h"]) er = 0.5 * np.hypot(t["e_mg_h"], t["e_si_h"]) elif ref == "Mg": r, er = t["mg_h"], t["e_mg_h"] elif ref == "Si": r, er = t["si_h"], t["e_si_h"] else: raise ValueError(ref) out = t.copy() out["rock_h"] = r out["e_rock_h"] = er for x in ("th", "eu", "k"): if f"{x}_h" in out: out[f"{x}_rock"] = out[f"{x}_h"] - r out[f"e_{x}_rock"] = np.hypot(out[f"e_{x}_h"], er) return out # --------------------------------------------------------------------------------------------- # linear model with intrinsic scatter: y_i ~ N(X_i . coef, s^2 + e_i^2) # Maximum likelihood by profiling: for a given s the coefficients are the weighted least-squares # solution, so only a one-dimensional bounded search over ln s is needed. Errors: bootstrap. # --------------------------------------------------------------------------------------------- def _ml_linear(X, y, e): def solve(ls): w = 1.0 / (np.exp(2.0 * ls) + e ** 2) XtW = X.T * w coef = np.linalg.solve(XtW @ X, XtW @ y) r = y - X @ coef return 0.5 * np.sum(r ** 2 * w - np.log(w)), coef res = minimize_scalar(lambda ls: solve(ls)[0], bounds=(np.log(1e-4), np.log(1.0)), method="bounded", options={"xatol": 1e-5}) return solve(res.x)[1], float(np.exp(res.x)) def _fit_with_bootstrap(X, y, e, n_boot, seed): coef, s = _ml_linear(X, y, e) rng = np.random.default_rng(seed) boots = [] for _ in range(n_boot): i = rng.integers(0, len(y), len(y)) c, sb = _ml_linear(X[i], y[i], e[i]) boots.append([*c, sb]) boots = np.array(boots) return coef, s, boots.std(axis=0), boots def fit_line_scatter(x, y, e, x0=0.0, n_boot=300, seed=SEED): x, y, e = (np.asarray(v, dtype=float) for v in (x, y, e)) ok = np.isfinite(x) & np.isfinite(y) & np.isfinite(e) x, y, e = x[ok] - x0, y[ok], e[ok] X = np.column_stack([np.ones_like(x), x]) coef, s, err, _ = _fit_with_bootstrap(X, y, e, n_boot, seed) res = y - X @ coef return {"n": int(len(x)), "a": float(coef[0]), "b": float(coef[1]), "s_int": s, "e_a": float(err[0]), "e_b": float(err[1]), "e_s": float(err[2]), "rms": float(np.sqrt(np.mean(res ** 2))), "med_err": float(np.median(e)), "x0": x0} def fit_plane_scatter(x1, x2, y, e, x10=0.0, x20=0.0, n_boot=300, seed=SEED): """y_i ~ N(a + b (x1 - x10) + c (x2 - x20), s^2 + e_i^2); bootstrap errors.""" x1, x2, y, e = (np.asarray(v, dtype=float) for v in (x1, x2, y, e)) ok = np.isfinite(x1) & np.isfinite(x2) & np.isfinite(y) & np.isfinite(e) x1, x2, y, e = x1[ok] - x10, x2[ok] - x20, y[ok], e[ok] X = np.column_stack([np.ones_like(x1), x1, x2]) coef, s, err, _ = _fit_with_bootstrap(X, y, e, n_boot, seed) return {"n": int(len(y)), "a": float(coef[0]), "b": float(coef[1]), "c": float(coef[2]), "s_int": s, "e_a": float(err[0]), "e_b": float(err[1]), "e_c": float(err[2]), "e_s": float(err[3]), "x10": x10, "x20": x20} def robust_trend_scatter(x, y, x0=0.0, clip=4.0): """For the large GALAH samples: ordinary least squares with sigma clipping; returns slope, intercept and the robust (MAD-based) residual scatter. Measurement errors are subtracted in quadrature by the caller.""" x = np.asarray(x, dtype=float) - x0 y = np.asarray(y, dtype=float) ok = np.isfinite(x) & np.isfinite(y) x, y = x[ok], y[ok] keep = np.ones(len(x), bool) for _ in range(5): b, a = np.polyfit(x[keep], y[keep], 1) res = y - a - b * x mad = 1.4826 * np.median(np.abs(res[keep] - np.median(res[keep]))) keep = np.abs(res) < clip * mad n = keep.sum() sxx = np.sum((x[keep] - x[keep].mean()) ** 2) return {"n": int(n), "a": a, "b": b, "scatter": mad, "e_b": mad / np.sqrt(sxx), "e_a": mad / np.sqrt(n), "x0": x0} # --------------------------------------------------------------------------------------------- # main analysis # --------------------------------------------------------------------------------------------- def build_samples(ref="MgSi"): tw = add_rock_ratios(load_solar_twins(), ref) dm = add_rock_ratios(load_delgado_mena(), ref) bb = add_rock_ratios(load_battistini(), ref) g3 = add_rock_ratios(load_galah(3), ref) g4 = add_rock_ratios(load_galah(4), ref) # --- solar twins: initial Th (decay-corrected), all relative to the Sun ------------------- tw["logf_th_direct"] = tw["th_rock"] + DEX_PER_GYR_TH * (tw["age"] - AGE_SUN) tw["e_logf_th_direct"] = np.hypot(tw["e_th_rock"], DEX_PER_GYR_TH * tw["e_age"]) tw["th_eu_zams"] = tw["th_h"] - tw["eu_h"] + DEX_PER_GYR_TH * (tw["age"] - AGE_SUN) tw["e_th_eu_zams"] = np.sqrt(tw["e_th_h"] ** 2 + tw["e_eu_h"] ** 2 + (DEX_PER_GYR_TH * tw["e_age"]) ** 2) for t in (tw, dm, bb): # DATA ONLY: europium per unit rock, relative to the Sun t["logf_eu_rock"] = t["eu_rock"] t["e_logf_eu_rock"] = t["e_eu_rock"] # DATA + MODEL: Eu per rock plus the one-zone free-decay factor for Th/Eu at the star's birth # (full strength). The factor alone contributes about +0.010 dex/Gyr to the age trend. t["logf_th_from_eu"] = t["eu_rock"] + ism_factor("Th232", t["age"]) t["e_logf_th_from_eu"] = t["e_eu_rock"] # raw ratios to each reference element separately (the trend depends on which is used) t["eu_mg"] = t["eu_h"] - t["mg_h"] t["eu_si"] = t["eu_h"] - t["si_h"] # --- selection cuts (stated in the report) ------------------------------------------------ tw["use_th"] = tw["th_h"].notna() & (tw["th_published"] == True) & ~tw["excluded_botelho"] # noqa: E712 tw["use_eu"] = tw["eu_h"].notna() & ~tw["excluded_botelho"] dm["use_eu"] = (dm["eu_h"].notna() & dm["e_eu_h"].notna() & (dm["e_eu_h"] <= 0.15) & dm["Teff"].between(5300, 6400) & (dm["logg"] > 3.8) & dm["age"].notna() & (dm["e_age"] <= 2.0) & (dm["feh"] > -1.0)) bb["use_eu"] = (bb["eu_h"].notna() & bb["Teff"].between(5300, 6400) & (bb["logg"] > 3.8) & bb["age"].notna() & (bb["e_age"] <= 2.0) & (bb["feh"] > -1.0)) for g in (g3, g4): g["use_k"] = ((g["flag_k"] == 0) & (g["flag_mg"] == 0) & (g["flag_si"] == 0) & g["age"].notna() & (g["snr"] >= 60) & g["Teff"].between(5500, 6100) & g["logg"].between(3.8, 4.6) & (g["feh"] > -1.0)) g3["use_k"] &= (g3["e_age"] / g3["age"].clip(lower=0.1) <= 0.3) g4["use_k"] &= g4["age"].between(0.5, 13.5) for g in (g3, g4): g["logf_k_star"] = g["k_rock"] # stellar [K/R] only; ISM 40K factor and volatility added later return {"solar_twins": tw, "harps_gto": dm, "battistini_bensby": bb, "galah_dr3": g3, "galah_dr4": g4} def analyse(ref="MgSi", verbose=True, n_boot=300): S = build_samples(ref) tw, dm, bb, g3, g4 = (S[k] for k in ("solar_twins", "harps_gto", "battistini_bensby", "galah_dr3", "galah_dr4")) rows = [] def add(label, quantity, xname, fit, extra=""): rows.append({"dataset": label, "quantity": quantity, "x": xname, "n": fit["n"], "value_at_x0": fit["a"], "e_value": fit["e_a"], "slope": fit["b"], "e_slope": fit["e_b"], "intrinsic_scatter_dex": fit.get("s_int", np.nan), "e_scatter": fit.get("e_s", np.nan), "rms_or_robust_scatter": fit.get("rms", fit.get("scatter", np.nan)), "median_meas_err": fit.get("med_err", np.nan), "x0": fit["x0"], "note": extra}) # ---- Th, direct (solar twins) ---- m = tw["use_th"] f = fit_line_scatter(tw.loc[m, "age"], tw.loc[m, "logf_th_direct"], tw.loc[m, "e_logf_th_direct"], x0=AGE_SUN, n_boot=n_boot) add("Botelho+2019 solar twins", "log f_Th (direct, decay-corrected)", "age [Gyr]", f) fit_th_age = f f = fit_line_scatter(tw.loc[m, "feh"], tw.loc[m, "logf_th_direct"], tw.loc[m, "e_logf_th_direct"], x0=0.0, n_boot=n_boot) add("Botelho+2019 solar twins", "log f_Th (direct, decay-corrected)", "[Fe/H]", f) # [Th/Fe] and [Th/H] at ZAMS versus age: the "more r-process per iron in older stars" check y = tw["th_h"] - tw["feh"] + DEX_PER_GYR_TH * (tw["age"] - AGE_SUN) f = fit_line_scatter(tw.loc[m, "age"], y[m], np.hypot(tw.loc[m, "e_th_h"], DEX_PER_GYR_TH * tw.loc[m, "e_age"]), x0=AGE_SUN, n_boot=n_boot) add("Botelho+2019 solar twins", "[Th/Fe] at birth", "age [Gyr]", f, "Botelho+ report no trend, mean +0.086") # Th/Eu at birth versus age: test of the uniform-production ISM model f = fit_line_scatter(tw.loc[m, "age"], tw.loc[m, "th_eu_zams"], tw.loc[m, "e_th_eu_zams"], x0=AGE_SUN, n_boot=n_boot) pred = float((ism_factor("Th232", 8.0) - ism_factor("Th232", 1.0)) / 7.0) add("Botelho+2019 x Spina+2018", "[Th/Eu] at birth", "age [Gyr]", f, f"uniform-production model predicts slope {pred:+.4f} dex/Gyr") fit_theu = f theu_model_slope = pred # ---- Eu per rock: data only (raw) and with the full-strength model Th/Eu factor ---- fits_eu, fits_eu_raw = {}, {} for label, t, key in (("Spina+2018 solar twins", tw, "solar_twins"), ("Delgado Mena+2017/2019 HARPS-GTO", dm, "harps_gto"), ("Battistini & Bensby 2016", bb, "battistini_bensby")): m = t["use_eu"] f0 = fit_line_scatter(t.loc[m, "age"], t.loc[m, "logf_eu_rock"], t.loc[m, "e_logf_eu_rock"], x0=AGE_SUN, n_boot=n_boot) add(label, "[Eu/rock] (observed, data only)", "age [Gyr]", f0, "no model factor") fits_eu_raw[key] = f0 f0z = fit_line_scatter(t.loc[m, "feh"], t.loc[m, "logf_eu_rock"], t.loc[m, "e_logf_eu_rock"], x0=0.0, n_boot=n_boot) add(label, "[Eu/rock] (observed, data only)", "[Fe/H]", f0z, "no model factor") f = fit_line_scatter(t.loc[m, "age"], t.loc[m, "logf_th_from_eu"], t.loc[m, "e_logf_th_from_eu"], x0=AGE_SUN, n_boot=n_boot) add(label, "log f_Th (from Eu + full one-zone decay factor)", "age [Gyr]", f, "data + model") fits_eu[key] = f f2 = fit_line_scatter(t.loc[m, "feh"], t.loc[m, "logf_th_from_eu"], t.loc[m, "e_logf_th_from_eu"], x0=0.0, n_boot=n_boot) add(label, "log f_Th (from Eu + full one-zone decay factor)", "[Fe/H]", f2, "data + model") y = t["eu_h"] - t["feh"] f3 = fit_line_scatter(t.loc[m, "age"], y[m], t.loc[m, "e_eu_h"], x0=AGE_SUN, n_boot=n_boot) add(label, "[Eu/Fe] (observed)", "age [Gyr]", f3, "r-process per iron versus age") # ---- K (GALAH) ---- fits_k, galah_eu = {}, {} for label, g, key in (("GALAH DR3", g3, "galah_dr3"), ("GALAH DR4", g4, "galah_dr4")): m = g["use_k"] f = robust_trend_scatter(g.loc[m, "age"], g.loc[m, "k_rock"], x0=AGE_SUN) med_err = float(np.median(g.loc[m, "e_k_rock"])) f["s_int"] = float(np.sqrt(max(f["scatter"] ** 2 - med_err ** 2, 0.0))) f["med_err"] = med_err add(label, "[K/R] (observed, stable K)", "age [Gyr]", f, "intrinsic scatter = robust scatter minus median quoted error in quadrature") fits_k[key] = f fz = robust_trend_scatter(g.loc[m, "feh"], g.loc[m, "k_rock"], x0=0.0) fz["med_err"] = med_err add(label, "[K/R] (observed, stable K)", "[Fe/H]", fz) f["n_pass_cuts"] = int(m.sum()) # GALAH europium in Sun-like dwarfs: detection statistics (reported as unusable). Denominator: # ALL Sun-like dwarfs passing the star-level cuts (count query), not the pre-filtered download. det = (g["flag_eu"] == 0) n_all = int(pd.read_csv(os.path.join(RAW_DIR, f"{key}_sunlike_count.csv")).iloc[0, 0]) galah_eu[key] = {"n_detected": int(det.sum()), "n_sunlike_all": n_all, "n_download_subset": int(len(g)), "fraction_detected": float(det.sum() / n_all), "median_eu_fe_of_detections": float(np.median(g.loc[det, "[Eu/Fe]"]))} rows.append({"dataset": label, "quantity": "[Eu/Fe] detections among Sun-like dwarfs (unusable: detection-limited)", "x": "-", "n": int(det.sum()), "value_at_x0": float(np.median(g.loc[det, "[Eu/Fe]"])), "note": f"detected in {det.sum() / n_all * 100:.2f}% of {n_all} Sun-like dwarfs passing the star-level cuts; median [Eu/Fe] of detections"}) table = pd.DataFrame(rows) table.to_csv(os.path.join(DATA_DIR, f"abundance_fits_{ref}.csv"), index=False) # ---- per-star table of the small samples ---- keep = [] for key, t in (("solar_twins", tw), ("harps_gto", dm), ("battistini_bensby", bb)): cols = [c for c in ("HIP", "Star", "sample", "Teff", "logg", "feh", "age", "e_age", "mg_h", "si_h", "eu_h", "e_eu_h", "th_h", "e_th_h", "logf_th_direct", "e_logf_th_direct", "logf_eu_rock", "e_logf_eu_rock", "logf_th_from_eu", "e_logf_th_from_eu", "pop", "use_th", "use_eu") if c in t] keep.append(t[cols]) stars = pd.concat(keep, ignore_index=True) stars.to_csv(os.path.join(DATA_DIR, f"hpe_stars_{ref}.csv"), index=False) # ---- binned GALAH K table (small file for the figure) ---- bins = np.arange(0.0, 14.0, 1.0) gb = [] for key, g in (("galah_dr3", g3), ("galah_dr4", g4)): m = g["use_k"] idx = np.digitize(g.loc[m, "age"], bins) for i in range(1, len(bins)): v = g.loc[m, "k_rock"].values[idx == i] if len(v) >= 30: gb.append({"sample": key, "age_lo": bins[i - 1], "age_hi": bins[i - 1] + 1.0, "n": len(v), "median": np.median(v), "p16": np.percentile(v, 16), "p84": np.percentile(v, 84)}) pd.DataFrame(gb).to_csv(os.path.join(DATA_DIR, f"galah_k_binned_{ref}.csv"), index=False) summary = {"fit_th_age": fit_th_age, "fit_theu": fit_theu, "theu_model_slope": theu_model_slope, "fits_eu": fits_eu, "fits_eu_raw": fits_eu_raw, "fits_k": fits_k, "galah_eu": galah_eu, "n": {k: int(v[[c for c in v.columns if c.startswith("use_")][0]].sum()) for k, v in S.items()}, "n_th": int(tw["use_th"].sum()), "n_eu_twins": int(tw["use_eu"].sum())} if verbose: with pd.option_context("display.width", 250, "display.max_columns", 30, "display.max_colwidth", 60): print(table[["dataset", "quantity", "x", "n", "value_at_x0", "e_value", "slope", "e_slope", "intrinsic_scatter_dex", "e_scatter", "rms_or_robust_scatter", "median_meas_err"]].round(4).to_string()) return S, table, summary if __name__ == "__main__": analyse()