#!/usr/bin/env python3 """Head-start summary table, the Sun's rank, and the "Fermi clock". Inputs (written by the analysis scripts): data/xr_age_distribution.csv, data/xr_results.csv, data/gaia_sg_age_distribution.csv and the published SFHs transcribed in published_sfh.py. Outputs: data/headstart_summary.csv, data/fermi_clock.csv, data/sun_rank.csv, data/fermi_assumptions.csv, and (revision) data/headline_estimates.csv (every estimate of the fraction older than the Sun, and the adopted value with its systematic range), data/fermi_clock_range.csv (N_eff and the limit across that range), data/fermi_marginal_limit.csv (limit when the uncertainty of N_eff is marginalised), data/galaxy_wide_bracket.csv. ASSUMPTIONS (not measured here; all stated in the report) --------------------------------------------------------- N_STARS number of stars in the Milky Way, 2e11 (commonly quoted 1-4e11). IMF Kroupa (2001) broken power law, 0.08-120 Msun -> number fraction of stars born with 0.8-1.2 Msun (computed below). ETA_EARTH Bryson et al. (2021, AJ 161, 36): occurrence of 0.5-1.5 R_earth planets in the conservative habitable zone of 4800-6300 K stars, 0.37 (+0.48/-0.21) [model 1 lower bound] to 0.60 (+0.90/-0.36) [upper bound]. D_GAL 1e5 light-years: distance a settlement front must cover. Emergence "Earth clock": a planet that produces a settling civilisation does so at age 4.567 Gyr. "hard steps n": cumulative probability proportional to (planet age)^n up to the end of the habitable lifetime T_HAB (5.5 Gyr), normalised to p at Earth's present age. """ import os import numpy as np import pandas as pd import common as C import published_sfh as P T_SUN = C.T_SUN D_GAL_LY = 1.0e5 SPEEDS = [1e-5, 3e-5, 1e-4, 3e-4, 1e-3, 1e-2, 1e-1] # units of c N_STARS = (2.0e11, 1.0e11, 4.0e11) # value, low, high ETA = {"Bryson+21 conservative HZ, lower bound": (0.37, 0.16, 0.85), "Bryson+21 conservative HZ, upper bound": (0.60, 0.24, 1.50)} T_HAB = 5.5 POISSON95 = -np.log(0.05) # 2.996 NSUB = 50 # sub-bin sampling def kroupa_number_fraction(m1=0.8, m2=1.2, lo=0.08, hi=120.0): """Number fraction of stars born with m1 < M < m2 (Kroupa 2001: slope 1.3 for 0.08-0.5 Msun, 2.3 above).""" def integ(a, b): tot = 0.0 for (x1, x2, s, k) in [(0.08, 0.5, 1.3, 1.0), (0.5, 1e3, 2.3, 0.5 ** (2.3 - 1.3))]: u1, u2 = max(a, x1), min(b, x2) if u2 > u1: tot += k * (u1 ** (1 - s) - u2 ** (1 - s)) / (s - 1) return tot return integ(m1, m2) / integ(lo, hi) def fine_ages(phi, nsub=NSUB): """Sample the piecewise-uniform age distribution on a fine grid.""" lo, hi = C.AGE_EDGES[:-1], C.AGE_EDGES[1:] f = (np.arange(nsub) + 0.5) / nsub ages = (lo[:, None] + (hi - lo)[:, None] * f[None, :]).ravel() w = np.repeat(phi / phi.sum() / nsub, nsub) return ages, w def emergence_weight(planet_age, model): """c(t): cumulative probability (in units of p = value at Earth's age) that a planet of age t has launched a settling civilisation.""" t = np.asarray(planet_age, float) if model == "Earth clock": return (t >= T_SUN).astype(float) n = int(model.split("=")[1]) return np.where(t > 0, (np.clip(t, 0, T_HAB) / T_SUN) ** n, 0.0) MODELS = ["Earth clock", "hard steps n=1", "hard steps n=3", "hard steps n=5"] def clock_factor(phi, t_settle, model): """Sum over planets of c(age - t_settle): the number of 'Earth-equivalent chances' per planet whose settlement wave would already be complete.""" ages, w = fine_ages(phi) return float(np.sum(w * emergence_weight(ages - t_settle, model))) XR = "LAMOST subgiants (X&R22) | " NAT = "Nataf+24 subgiants | " GSG = "Gaia 200 pc subgiants (this work) | " FID = XR + "formed: all [Fe/H]" def _read(name): path = C.DATA + "/" + name return pd.read_csv(path) if os.path.exists(path) else None def distributions(): """name -> (kind, phi, frac_pass). frac_pass: fraction of all Sun-like stars that pass the metallicity floor of that variant. kind: "tracer" (column, stars ever formed), "tracer-sub" (a living or mid-plane or sub-population variant, not used for the clock), "column" / "midplane" (published SFHs).""" out = {} for prefix, dist_file, res_file in ( (XR, "xr_age_distribution.csv", "xr_results.csv"), (NAT, "nataf_age_distribution.csv", "nataf_results.csv"), (GSG, "gaia_sg_age_distribution.csv", "gaia_sg_results.csv")): d, r = _read(dist_file), _read(res_file) if d is None or r is None: continue r = r.set_index("variant") for k in d.columns: if k in ("age_lo", "age_hi") or k == "no error deconvolution": continue sub = ("living" in k or "mid-plane" in k or "low-alpha" in k) out[prefix + k] = ("tracer-sub" if sub else "tracer", d[k].values, float(r.loc[k, "frac_pass"])) for name, (kind, phi) in P.all_published().items(): out["published SFH | " + name] = (kind, phi, 1.0) return out def summary_table(dists): rows = [] xr = pd.read_csv(C.DATA + "/xr_results.csv").set_index("variant") gs = pd.read_csv(C.DATA + "/gaia_sg_results.csv").set_index("variant") nat = _read("nataf_results.csv") nat = nat.set_index("variant") if nat is not None else None for name, (kind, phi, fp) in dists.items(): st = C.headstart_stats(phi) row = dict(distribution=name, kind=kind, frac_pass_floor=fp, **st) row["young_lt2"] = P.young_fraction(phi) src, _, var = name.partition(" | ") tab = (xr if src.startswith("LAMOST") else gs if src.startswith("Gaia") else nat if src.startswith("Nataf") else None) if tab is not None and var in tab.index and "f_gt_0_lo" in tab.columns: for k in ["f_gt_0", "f_gt_0.1", "f_gt_1", "f_gt_2", "f_gt_3", "f_gt_5", "mean_headstart_older", "median_headstart_older", "mean_offset_all"]: if k + "_lo" in tab.columns: row[k + "_lo"] = tab.loc[var, k + "_lo"] row[k + "_hi"] = tab.loc[var, k + "_hi"] rows.append(row) return pd.DataFrame(rows) def clock_table(dists, n_sun): """N_eff and the 95% upper limit on p for every distribution, emergence model and expansion speed, for both Bryson+21 eta_Earth bounds.""" rows = [] for name, (kind, phi, fp) in dists.items(): if kind in ("tracer-sub", "midplane"): continue # the clock needs stars ever formed, per column for model in MODELS: for v in SPEEDS: ts = D_GAL_LY / v / 1e9 # Gyr cf = clock_factor(phi, ts, model) for eta_name, (eta, eta_lo, eta_hi) in ETA.items(): n_pl = n_sun[0] * eta * fp neff = n_pl * cf rows.append(dict( distribution=name, model=model, v_over_c=v, t_settle_gyr=ts, clock_factor=cf, eta_earth=eta_name, n_planets=n_pl, n_eff=neff, p95=POISSON95 / neff if neff > 0 else np.inf, # extremes of N_sunlike and eta_Earth (1-sigma) together p95_low=(POISSON95 / (n_sun[2] * eta_hi * fp * cf) if cf > 0 else np.inf), p95_high=(POISSON95 / (n_sun[1] * eta_lo * fp * cf) if cf > 0 else np.inf))) return pd.DataFrame(rows) def sun_rank(dists): """The reverse question: where does the Sun sit?""" rows = [] for name, (kind, phi, fp) in dists.items(): ages, w = fine_ages(phi) f_old = float(w[ages > T_SUN].sum()) y2 = P.young_fraction(phi) # formed in the last 2 Gyr rate = y2 / 2.0 # per Gyr, recent row = dict(distribution=name, f_older_to_date=f_old, recent_rate_per_gyr=rate) # future star formation (ASSUMPTION): recent rate continues for T_f for tf in (5.0, 10.0): row[f"f_older_incl_future_const_{tf:.0f}Gyr"] = f_old / (1 + rate * tf) # (an exponentially declining future rate with e-folding time T_f # forms the same number of stars as a constant rate lasting T_f) # Observer-centred statistic. If observers appear when their star is # 4.567 Gyr old, an observer on a star of present age a saw, at that # moment, a fraction r(a) = P(age > a) / P(age > a - 4.567) of the # then-existing Sun-like stars older than its own. r(Sun) = f_old. order = np.argsort(ages) a_s, w_s = ages[order], w[order] surv = 1.0 - np.cumsum(w_s) + 0.5 * w_s # P(age > a) r = surv / np.interp(a_s - T_SUN, a_s, surv, left=1.0) past = a_s >= T_SUN # already observers ww = w_s[past] / w_s[past].sum() o = np.argsort(r[past]) cdf = np.cumsum(ww[o]) q = np.interp([0.16, 0.5, 0.84], cdf, r[past][o]) row.update(obs_r_p16=q[0], obs_r_median=q[1], obs_r_p84=q[2]) # slowest expansion speeds that the clock can still constrain hs = np.sort(ages[ages > T_SUN]) - T_SUN wh = w[np.argsort(ages)][np.sort(ages) > T_SUN] ch = np.cumsum(wh) / wh.sum() h_med = float(np.interp(0.5, ch, hs)) h_max = float(hs[np.searchsorted(ch, 0.999)]) row.update(v_half_over_c=D_GAL_LY / (h_med * 1e9), v_min_over_c=D_GAL_LY / (h_max * 1e9)) rows.append(row) return pd.DataFrame(rows) # Galaxy-wide brackets (ASSUMPTIONS, illustrative only). # (1) The solar cylinder is taken as representative of the disc; a bulge/bar # holding F_BULGE of the Galaxy's stars (Bland-Hawthorn & Gerhard 2016: # ~1.4-1.7e10 of ~5e10 Msun) is assumed to have a fraction F_BULGE_OLD of # its Sun-like stars older than the Sun. # (2) Xiang+ 2025 (arXiv:2410.09705) measure, with a modelled selection # function, a present-day stellar mass of 2e10 Msun for the whole # high-alpha disc, formed mostly ~11 Gyr ago. With a total stellar mass # of 4-6e10 Msun (Bland-Hawthorn & Gerhard 2016: 5 +/- 1e10; Licquia & # Newman 2015: 6.1 +/- 1.1e10) the high-alpha stars are 33-50% of the # Galaxy's stars; all of them are older than the Sun. The remaining # (low-alpha) stars are given the age mix of the local low-alpha stars. # Present-day mass shares understate the share of stars ever formed for # old populations, so this is conservative in that respect; radial age # gradients of the low-alpha disc are ignored. F_BULGE = (0.25, 0.15, 0.35) F_BULGE_OLD = (0.90, 0.80, 1.00) F_HIGH_ALPHA_GALAXY = (0.40, 0.33, 0.50) def galaxy_bracket(f_local): out = [] for fb in F_BULGE: for fo in F_BULGE_OLD: out.append((1 - fb) * f_local + fb * fo) return out[0], min(out), max(out) def galaxy_bracket_alpha(f_low_alpha_local): out = [F + (1 - F) * f_low_alpha_local for F in F_HIGH_ALPHA_GALAXY] return out[0], min(out), max(out) # ------------------------------------------------ marginalised upper limit def p95_marginal_factor(sigma_ln, n_grid=4001): """Factor k by which the 95% limit p < 3.0 / N0 loosens when N_eff is not known but log-normally distributed around N0 (median) with dispersion sigma_ln, for a flat prior on p and zero events: posterior(p) ~ E_N[exp(-p N)], CDF(p) = E[(1 - exp(-p N)) / N] / E[1 / N]. Returns k = p95 * N0 / 3.0.""" from scipy.optimize import brentq if sigma_ln <= 0: return 1.0 x = np.linspace(-8, 8, n_grid) wq = np.exp(-0.5 * x * x) wq /= wq.sum() N = np.exp(sigma_ln * x) # in units of N0 den = np.sum(wq / N) def cdf(p): return np.sum(wq * (1 - np.exp(-p * N)) / N) / den - 0.95 return brentq(cdf, 1e-3, 1e6) / POISSON95 # sigma_ln(N_eff): eta_Earth 0.37 (+0.48/-0.21) is a 68% interval, i.e. # sigma_ln = 0.83; the number of Sun-like stars is "a factor 2 either way", # read either as 2 sigma (0.35) or as 1 sigma (0.69); the fraction with # enough head start adds 0.07. SIGMA_LN_CASES = [ (0.80, "referee's lower case"), (float(np.sqrt(0.83 ** 2 + 0.35 ** 2 + 0.07 ** 2)), "eta_Earth 68% range; star count factor 2 = 2 sigma"), (1.00, "referee's upper case"), (float(np.sqrt(0.83 ** 2 + 0.69 ** 2 + 0.07 ** 2)), "eta_Earth 68% range; star count factor 2 = 1 sigma"), ] # --------------------------------------------------- headline estimates def headline_estimates(): """Every estimate of the column fraction of Sun-like stars older than the Sun made in this study, and the adopted value. tier = "core": the two catalogue samples with published individual ages and a selection box that can be evaluated on isochrones (LAMOST / Xiang & Rix; Nataf+24 primary sample), each on the four published age scales, plus the star-by-star age swap and the Nataf sample restricted to tracers older than 3 Gyr. "bracket": samples that are weaker but not refuted: this study's Gaia-only 200 pc sample rescaled onto published scales (low side) and the Nataf+24 sample including its metal-poor annex (high side). "excluded": known to be biased or disfavoured by a test; listed with the reason. The adopted value is the midpoint of the core envelope; the adopted systematic range is the envelope of core + bracket.""" xr = _read("xr_results.csv").set_index("variant") gs = _read("gaia_sg_results.csv").set_index("variant") nat = _read("nataf_results.csv") nat = nat.set_index("variant") if nat is not None else None swap = _read("age_swap_test.csv") rows = [] def add(label, tab, var, tier, why="", group=""): if tab is None or var not in tab.index: return r = tab.loc[var] rows.append(dict( estimate=label, group=group, tier=tier, f_older=r.f_gt_0, stat_lo=r.get("f_gt_0_lo", np.nan), stat_hi=r.get("f_gt_0_hi", np.nan), f_gt_1=r.f_gt_1, f_gt_3=r.f_gt_3, f_gt_5=r.f_gt_5, mean_headstart_older=r.mean_headstart_older, median_headstart_older=r.median_headstart_older, mean_offset_all=r.mean_offset_all, n=int(r.n), in_range=tier != "excluded", comment=why)) g = "LAMOST subgiants (Xiang & Rix 2022)" add("LAMOST, X&R ages, dynamical column weights (primary)", xr, "formed: all [Fe/H]", "core", group=g) for v in ("ages on the Gaia FLAME scale", "ages on the Nataf+24 scale", "ages on the StarHorse scale"): add("LAMOST, " + v, xr, v, "core", "constant factor measured for stars in common", group=g) if swap is not None: t = swap[swap.subset == "in common, Nataf quality flags"].set_index("ages") d = float(t.loc["Nataf+24 ages", "f_older_column"] - t.loc["X&R ages", "f_older_column"]) base = xr.loc["formed: all [Fe/H]"] rows.append(dict( estimate="LAMOST, Nataf+24 ages substituted star by star " "(paired shift applied to the primary)", group=g, tier="core", f_older=base.f_gt_0 + d, stat_lo=np.nan, stat_hi=np.nan, f_gt_1=np.nan, f_gt_3=np.nan, f_gt_5=np.nan, mean_headstart_older=np.nan, median_headstart_older=np.nan, mean_offset_all=np.nan, n=int(t.n.iloc[0]), in_range=True, comment=f"paired difference {d:+.3f} on the " f"{int(t.n.iloc[0])} stars in common")) add("LAMOST, star-count scale heights (Mazzi+24 exp)", xr, "star-count heights (Mazzi+24 exp)", "excluded", "disfavoured: under-predicts the old fraction LAMOST observes at " "200-800 pc (layer test)", group=g) add("LAMOST, star-count scale heights (Mazzi+24 sech2)", xr, "star-count heights (Mazzi+24 sech2)", "excluded", "disfavoured by the layer test", group=g) add("LAMOST, star-count heights, thick disc 600/900 pc", xr, "star-count heights (exp, thick disc 600/900 pc)", "excluded", "disfavoured by the layer test, less strongly", group=g) add("LAMOST, ages x0.87", xr, "ages x0.87 (referee: X&R / Gaia-only ages)", "excluded", "no published catalogue is this much younger than X&R", group=g) add("LAMOST, ages x1.146", xr, "ages x1.146 (referee: Gaia-only scale)", "excluded", "scale of this study's Gaia-only ages, which are biased " "old", group=g) g = "Nataf+ 2024 subgiants" add("Nataf+24 primary sample, own ages", nat, "formed: all [Fe/H]", "core", group=g) add("Nataf+24 primary sample, ages on the X&R scale", nat, "ages on the X&R22 scale", "core", group=g) add("Nataf+24 primary sample, ages on the Gaia FLAME scale", nat, "ages on the Gaia FLAME scale", "core", group=g) add("Nataf+24 primary sample, ages on the StarHorse scale", nat, "ages on the StarHorse scale", "core", group=g) add("Nataf+24 primary sample, tracers older than 3 Gyr", nat, "tracer ages > 3 Gyr (young prior for 3 Gyr)", "core", "young fraction (< 3 Gyr) from published SFHs", group=g) add("Nataf+24 primary + metal-poor annex, own ages", nat, "primary + annex", "bracket", "annex selection only approximately modelled; annex stars are " "kinematically cold for their ages (possible contamination)", group=g) add("Nataf+24 primary + metal-poor annex, ages on the X&R scale", nat, "primary + annex, ages on the X&R22 scale", "bracket", "as above", group=g) g = "Gaia 200 pc subgiants (this study's own ages)" add("Gaia 200 pc, ages as computed", gs, "formed: all [M/H]", "excluded", "ages more than 10% older than three published sets for the same stars", group=g) add("Gaia 200 pc, ages on the X&R scale", gs, "ages on the X&R22 scale: formed, all [M/H]", "bracket", "constant rescaling of home-made ages; high binary fraction", group=g) add("Gaia 200 pc, ages on the Nataf+24 scale", gs, "ages on the Nataf+24 scale: formed, all [M/H]", "bracket", "as above", group=g) add("Gaia 200 pc, ages on the StarHorse scale", gs, "ages on the StarHorse scale: formed, all [M/H]", "bracket", "as above", group=g) df = pd.DataFrame(rows) core = df[df.tier == "core"] env = df[df.tier != "excluded"] lo, hi = float(env.f_older.min()), float(env.f_older.max()) clo, chi = float(core.f_older.min()), float(core.f_older.max()) prim = df.iloc[0] adopted = dict( f_adopted=round(0.5 * (clo + chi), 2), f_core_lo=clo, f_core_hi=chi, f_range_lo=lo, f_range_hi=hi, f_primary=float(prim.f_older), primary_stat_lo=float(prim.stat_lo), primary_stat_hi=float(prim.stat_hi), n_core=int(len(core)), n_in_range=int(len(env)), core_lowest=core.loc[core.f_older.idxmin(), "estimate"], core_highest=core.loc[core.f_older.idxmax(), "estimate"], lowest=env.loc[env.f_older.idxmin(), "estimate"], highest=env.loc[env.f_older.idxmax(), "estimate"]) for key, col in (("f1", "f_gt_1"), ("f3", "f_gt_3"), ("f5", "f_gt_5"), ("mean_hs", "mean_headstart_older"), ("median_hs", "median_headstart_older"), ("offset_all", "mean_offset_all")): nm = key + ("_range" if key.startswith("f") else "") adopted[nm + "_lo"] = float(env[col].min()) adopted[nm + "_hi"] = float(env[col].max()) adopted[key + "_core_lo"] = float(core[col].min()) adopted[key + "_core_hi"] = float(core[col].max()) return df, adopted def main(): frac = kroupa_number_fraction() n_sun = tuple(n * frac for n in N_STARS) dists = distributions() summ = summary_table(dists) summ.to_csv(C.DATA + "/headstart_summary.csv", index=False) clock = clock_table(dists, n_sun) clock.to_csv(C.DATA + "/fermi_clock.csv", index=False) rank = sun_rank(dists) rank.to_csv(C.DATA + "/sun_rank.csv", index=False) pd.DataFrame([ dict(quantity="N_stars Milky Way (assumed)", value=N_STARS[0], low=N_STARS[1], high=N_STARS[2]), dict(quantity="Kroupa IMF number fraction 0.8-1.2 Msun", value=frac, low=np.nan, high=np.nan), dict(quantity="N Sun-like stars ever formed", value=n_sun[0], low=n_sun[1], high=n_sun[2]), dict(quantity="eta_Earth (Bryson+21 cons. HZ lower bound)", value=0.37, low=0.16, high=0.85), dict(quantity="eta_Earth (Bryson+21 cons. HZ upper bound)", value=0.60, low=0.24, high=1.50), dict(quantity="settling distance [ly]", value=D_GAL_LY, low=np.nan, high=np.nan), dict(quantity="T_hab for hard-step models [Gyr]", value=T_HAB, low=np.nan, high=np.nan), ]).to_csv(C.DATA + "/fermi_assumptions.csv", index=False) pd.set_option("display.width", 250, "display.max_columns", 40, "display.max_colwidth", 70) # ---- every estimate of the headline number, and the adopted value est, adopted = headline_estimates() est.to_csv(C.DATA + "/headline_estimates.csv", index=False) pd.DataFrame([adopted]).to_csv(C.DATA + "/headline_adopted.csv", index=False) print("[clock] estimates of the fraction older than the Sun:") print(est[["estimate", "tier", "f_older", "stat_lo", "stat_hi", "f_gt_5", "mean_headstart_older"]].round(3).to_string(index=False)) print("[clock] adopted: f = %.2f; core range %.3f-%.3f (%s; %s); full " "range %.3f-%.3f (lowest: %s; highest: %s)" % (adopted["f_adopted"], adopted["f_core_lo"], adopted["f_core_hi"], adopted["core_lowest"], adopted["core_highest"], adopted["f_range_lo"], adopted["f_range_hi"], adopted["lowest"], adopted["highest"])) # ---- Fermi clock across the range of age distributions c = clock[clock.eta_earth.str.contains("lower") & (clock.model == "Earth clock")] piv = c.pivot(index="distribution", columns="v_over_c", values="n_eff") piv_p = c.pivot(index="distribution", columns="v_over_c", values="p95") # all-metallicity column distributions on published age scales: the # tracer variants that enter the adopted range plus the published SFHs def allmet(d): if ">-0." in d: return False # metallicity floors if d.startswith(GSG) and "scale" not in d: return False # Gaia-only ages as computed if "star-count" in d: return False # disfavoured vertical model return True keep = [d for d in piv.index if allmet(d)] floors = [d for d in piv.index if ">-0." in d and "star-count" not in d and (not d.startswith(GSG) or "scale" in d)] rng = pd.DataFrame({ "v_over_c": piv.columns, "n_eff_fiducial": piv.loc[FID].values, "n_eff_min": piv.loc[keep].min().values, "n_eff_max": piv.loc[keep].max().values, "p95_fiducial": piv_p.loc[FID].values, "p95_min": piv_p.loc[keep].min().values, "p95_max": piv_p.loc[keep].max().values, "n_distributions": len(keep), "n_eff_min_with_floor": piv.loc[floors].min().values, "n_eff_max_with_floor": piv.loc[floors].max().values, "p95_min_with_floor": piv_p.loc[floors].min().values, "p95_max_with_floor": piv_p.loc[floors].max().values}) rng.to_csv(C.DATA + "/fermi_clock_range.csv", index=False) pd.DataFrame({"distribution": keep}).assign(set="all metallicities") \ .to_csv(C.DATA + "/fermi_clock_range_members.csv", index=False) print("[clock] N_eff and conditional p95 across %d all-metallicity " "column distributions (tracer samples on published age scales " "and published SFHs); '_with_floor': %d metallicity-floor " "variants" % (len(keep), len(floors))) print(rng.to_string(index=False, float_format=lambda v: f"{v:.3g}")) # ---- limit with the uncertainty of N_eff marginalised base = clock[(clock.distribution == FID) & (clock.model == "Earth clock") & clock.eta_earth.str.contains("lower")].set_index("v_over_c") mrows = [] for sig, why in SIGMA_LN_CASES: k = p95_marginal_factor(sig) mrows.append({"sigma_ln": sig, "meaning": why, "factor": k, "p95_at_1e-4c": k * base.loc[1e-4, "p95"], "p95_at_0.1c": k * base.loc[1e-1, "p95"]}) marg = pd.DataFrame(mrows) marg.to_csv(C.DATA + "/fermi_marginal_limit.csv", index=False) print("[clock] 95% limit with N_eff uncertain (lognormal, flat prior on p):") print(marg.to_string(index=False, float_format=lambda v: f"{v:.3g}")) # ---- Galaxy-wide brackets (assumptions) s_ = summ.set_index("distribution") f_loc = float(s_.loc[FID, "f_gt_0"]) xr_res = _read("xr_results.csv").set_index("variant") f_low = float(xr_res.loc["low-alpha stars only ([alpha/Fe] <= 0.16)", "f_gt_0"]) grows = [] for lab, f in (("primary (LAMOST)", f_loc), ("adopted", adopted["f_adopted"]), ("lower edge of range", adopted["f_range_lo"])): gb = galaxy_bracket(f) grows.append(dict(method="bulge assumption", f_local_used=lab, f_local=f, f_galaxy=gb[0], f_galaxy_min=gb[1], f_galaxy_max=gb[2])) ga = galaxy_bracket_alpha(f_low) grows.append(dict(method="Xiang+25 high-alpha disc mass + local " "low-alpha age mix", f_local_used="local low-alpha stars", f_local=f_low, f_galaxy=ga[0], f_galaxy_min=ga[1], f_galaxy_max=ga[2])) gal = pd.DataFrame(grows) gal.to_csv(C.DATA + "/galaxy_wide_bracket.csv", index=False) print("[clock] Galaxy-wide brackets (assumptions, illustrative):") print(gal.round(3).to_string(index=False)) print(f"[clock] Kroupa number fraction 0.8-1.2 Msun = {frac:.4f}; " f"N_sunlike = {n_sun[0]:.2e} ({n_sun[1]:.1e}-{n_sun[2]:.1e})") cols = ["distribution", "f_gt_0", "f_gt_0.1", "f_gt_1", "f_gt_2", "f_gt_3", "f_gt_5", "mean_headstart_older", "median_headstart_older", "mean_offset_all", "median_age", "frac_pass_floor"] print(summ[cols].round(3).to_string(index=False)) c = clock[(clock.distribution == FID) & clock.eta_earth.str.contains("lower")] print(c.pivot(index="v_over_c", columns="model", values="clock_factor") .round(3).to_string()) print(c.pivot(index="v_over_c", columns="model", values="p95") .to_string(float_format=lambda v: f"{v:.2e}")) print(rank.round(3).to_string(index=False)) return summ, clock, rank if __name__ == "__main__": main()