from __future__ import annotations import json from pathlib import Path import astropy from astropy.coordinates import SkyCoord import astropy.units as u import numpy as np import pandas as pd import scipy from scipy.integrate import trapezoid import sherpa from sherpa.astro.ui import set_xsabund, set_xsxsect, xsapec, xsphabs HERE = Path(__file__).resolve().parent ROOT = HERE.parent MEASUREMENTS = ( ROOT / "independent_cgm_results_v19_mos14to20keV_excluded_hi4pi_nh_per_pointing" / "independent_cgm_v19_all22_measurements.csv" ) LEDGER = HERE / "m31_cgmsum_conditional_prior_ledger.csv" OUTPUT = HERE / "m31_cgmsum_conditional_prior_audit.json" PREFERRED_FOOTPRINT_OUTPUT = ( HERE / "m31_cgmsum_locatelli2024_preferred_m31_footprint_predictions.csv" ) KPC_CM = 3.0856775814913673e21 ARCMIN2_PER_SR = (180.0 * 60.0 / np.pi) ** 2 def energy_flux_per_norm( kT_keV: float, abundance_solar: float, nh_1e22: float | None, low_keV: float, high_keV: float, ) -> float: plasma = xsapec.audit_plasma plasma.kT = kT_keV plasma.Abundanc = abundance_solar plasma.Redshift = 0.0 plasma.norm = 1.0 model = plasma if nh_1e22 is not None: absorption = xsphabs.audit_absorption absorption.nH = nh_1e22 model = absorption * plasma edges = np.linspace(low_keV, high_keV, 30001) midpoint = 0.5 * (edges[:-1] + edges[1:]) photons = model(edges[:-1], edges[1:]) return float(np.sum(photons * midpoint * 1.602176634e-9)) def emission_measure_kpc_cm6( longitude_rad: float, latitude_rad: float, *, beta: float, c_norm: float, n0: float, radial_scale_kpc: float, height_scale_kpc: float, ) -> float: distance_kpc = np.linspace(0.001, 350.0, 35001) solar_radius_kpc = 8.2 cylindrical_radius = np.sqrt( solar_radius_kpc**2 + (distance_kpc * np.cos(latitude_rad)) ** 2 - 2.0 * solar_radius_kpc * distance_kpc * np.cos(latitude_rad) * np.cos(longitude_rad) ) height = distance_kpc * np.sin(latitude_rad) spherical_radius = np.sqrt(cylindrical_radius**2 + height**2) spherical_density = ( c_norm * spherical_radius ** (-3.0 * beta) if c_norm > 0.0 else np.zeros_like(distance_kpc) ) disk_density = n0 * np.exp(-cylindrical_radius / radial_scale_kpc) * np.exp( -np.abs(height) / height_scale_kpc ) return float( trapezoid(spherical_density**2 + disk_density**2, distance_kpc) ) def weighted_mean(values: np.ndarray, errors: np.ndarray) -> float: weights = 1.0 / errors**2 return float(np.sum(weights * values) / np.sum(weights)) def predict_locatelli_fluxes( table: pd.DataFrame, *, beta: float, c_norm: float, n0: float, radial_scale_kpc: float, height_scale_kpc: float, kT_keV: float, abundance_solar: float, ) -> np.ndarray: fluxes = [] for _, row in table.iterrows(): emission_measure = emission_measure_kpc_cm6( float(row["l_rad"]), float(row["b_rad"]), beta=beta, c_norm=c_norm, n0=n0, radial_scale_kpc=radial_scale_kpc, height_scale_kpc=height_scale_kpc, ) emissivity = ( energy_flux_per_norm( kT_keV, abundance_solar, float(row["nh_hi4pi_1e22_cm-2"]), 0.4, 1.25, ) * 1.0e-14 ) fluxes.append( emissivity * emission_measure * KPC_CM / (4.0 * np.pi) / ARCMIN2_PER_SR / 1.0e-15 ) return np.asarray(fluxes) def side_estimators( table: pd.DataFrame, flux_primary_units: np.ndarray, error_name: str, ) -> dict[str, float]: means = {} for side in ("North/NW", "South/SE"): mask = table["side"].eq(side).to_numpy() means[side] = weighted_mean( flux_primary_units[mask], table.loc[mask, error_name].to_numpy() ) weights = 1.0 / table[error_name].to_numpy() ** 2 return { "north_primary_fluxunit": means["North/NW"], "south_primary_fluxunit": means["South/SE"], "side_balanced_primary_fluxunit": 0.5 * (means["North/NW"] + means["South/SE"]), "all_inverse_variance_primary_fluxunit": float( np.sum(weights * flux_primary_units) / np.sum(weights) ), } def main() -> None: set_xsabund("angr") set_xsxsect("vern") band_ratio = energy_flux_per_norm(0.175306, 0.3, 0.056, 0.4, 1.25) / energy_flux_per_norm( 0.175306, 0.3, None, 0.5, 2.0 ) variants = ( ("disk_like", 0.5, 0.0, 0.017, 8.0, 3.3, 0.15, 0.1), ("combined_beta0p3", 0.3, 0.0009, 0.017, 8.0, 3.3, 0.15, 0.1), ("combined_beta0p5", 0.5, 0.046, 0.032, 6.2, 1.1, 0.15, 0.1), ("combined_beta0p7", 0.7, 0.176, 0.016, 9.9, 1.6, 0.15, 0.1), ("combined_swcx", 0.5, 0.032, 0.063, 3.9, 0.9, 0.178, 0.1), ("combined_20pct_instr", 0.5, 0.043, 0.055, 4.6, 0.9, 0.15, 0.1), ("combined_high_emissivity", 0.5, 0.0136, 0.0094, 6.2, 1.1, 0.225, 0.3), ) table = pd.read_csv(MEASUREMENTS, dtype={"obsid": str}) table = table.loc[table["quality_primary"]].copy() galactic = SkyCoord( table["ra_deg"].to_numpy() * u.deg, table["dec_deg"].to_numpy() * u.deg, ).galactic table["l_rad"] = galactic.l.rad table["b_rad"] = galactic.b.rad error_name = "absorbed_flux_soft_0p4_1p25_staterr_erg_cm-2_s-1_arcmin-2" model_results: dict[str, dict[str, float]] = {} preferred_flux: np.ndarray | None = None for name, beta, c_norm, n0, radial_scale, height_scale, kT, abundance in variants: flux = predict_locatelli_fluxes( table, beta=beta, c_norm=c_norm, n0=n0, radial_scale_kpc=radial_scale, height_scale_kpc=height_scale, kT_keV=kT, abundance_solar=abundance, ) estimators = side_estimators(table, flux, error_name) characteristic = estimators["side_balanced_primary_fluxunit"] model_results[name] = { "north_primary_fluxunit": estimators["north_primary_fluxunit"], "south_primary_fluxunit": estimators["south_primary_fluxunit"], "characteristic_primary_fluxunit": characteristic, "all_inverse_variance_primary_fluxunit": estimators[ "all_inverse_variance_primary_fluxunit" ], "contrast_slope_D_over_T": ( estimators["north_primary_fluxunit"] - estimators["south_primary_fluxunit"] ) / characteristic, } if name == "combined_swcx": preferred_flux = flux if preferred_flux is None: raise RuntimeError("preferred combined_swcx model was not evaluated") preferred_estimators = side_estimators(table, preferred_flux, error_name) preferred_kwargs = { "beta": 0.5, "c_norm": 0.032, "n0": 0.063, "radial_scale_kpc": 3.9, "height_scale_kpc": 0.9, "kT_keV": 0.178, "abundance_solar": 0.1, } marginal_errors = { "c_norm": 0.001, "n0": 0.008, "radial_scale_kpc": 0.2, "height_scale_kpc": 0.1, } parameter_deltas: dict[str, float] = {} for parameter, sigma in marginal_errors.items(): low_kwargs = preferred_kwargs.copy() high_kwargs = preferred_kwargs.copy() low_kwargs[parameter] -= sigma high_kwargs[parameter] += sigma low_flux = predict_locatelli_fluxes(table, **low_kwargs) high_flux = predict_locatelli_fluxes(table, **high_kwargs) low_estimator = side_estimators(table, low_flux, error_name)[ "side_balanced_primary_fluxunit" ] high_estimator = side_estimators(table, high_flux, error_name)[ "side_balanced_primary_fluxunit" ] parameter_deltas[parameter] = 0.5 * (high_estimator - low_estimator) parameter_sigma = float( np.sqrt(np.sum(np.square(list(parameter_deltas.values())))) ) footprint = table[ [ "obsid", "side", "ra_deg", "dec_deg", "nh_hi4pi_1e22_cm-2", "l_rad", "b_rad", ] ].copy() footprint["galactic_l_deg"] = np.rad2deg(footprint.pop("l_rad")) footprint["galactic_b_deg"] = np.rad2deg(footprint.pop("b_rad")) footprint["locatelli_preferred_absorbed_0p4_1p25_primary_fluxunit"] = ( preferred_flux ) footprint.to_csv(PREFERRED_FOOTPRINT_OUTPUT, index=False) variant_values = np.array( [item["characteristic_primary_fluxunit"] for item in model_results.values()] ) summary = { "field_low": float(preferred_flux.min()), "field_p16": float(np.quantile(preferred_flux, 0.16)), "field_median": float(np.median(preferred_flux)), "field_p84": float(np.quantile(preferred_flux, 0.84)), "field_high": float(preferred_flux.max()), "north": preferred_estimators["north_primary_fluxunit"], "south": preferred_estimators["south_primary_fluxunit"], "side_balanced": preferred_estimators["side_balanced_primary_fluxunit"], "all_inverse_variance": preferred_estimators[ "all_inverse_variance_primary_fluxunit" ], "table1_independent_marginal_sigma": parameter_sigma, "table1_independent_marginal_low": preferred_estimators[ "side_balanced_primary_fluxunit" ] - parameter_sigma, "table1_independent_marginal_high": preferred_estimators[ "side_balanced_primary_fluxunit" ] + parameter_sigma, } ledger = pd.read_csv(LEDGER) expected = ledger.loc[ ledger["prior_id"].eq("locatelli2024_m31_footprint_model_family") ].iloc[0] for value, ledger_key in ( (summary["field_low"], "primary_low"), (summary["side_balanced"], "primary_central"), (summary["field_high"], "primary_high"), (summary["table1_independent_marginal_low"], "parameter_error_low"), (summary["table1_independent_marginal_high"], "parameter_error_high"), (summary["north"], "north_primary"), (summary["south"], "south_primary"), ): if not np.isclose(value, expected[ledger_key], rtol=2.0e-4): raise RuntimeError( f"Locatelli audit mismatch for {ledger_key}: " f"{value} != {expected[ledger_key]}" ) variant_envelope = { "low": float(variant_values.min()), "high": float(variant_values.max()), "status": "sensitivity context only; not used as the Figure 3 MW interval", } payload = { "schema": "m31-cgmsum-conditional-prior-audit-v2", "input_measurements": str(MEASUREMENTS.relative_to(ROOT)), "preferred_footprint_output": str(PREFERRED_FOOTPRINT_OUTPUT.relative_to(HERE)), "primary_field_count": int(len(table)), "spectral_convention": { "abundance_table": "angr", "cross_sections": "vern", "conversion_template_kT_keV": 0.175306, "conversion_template_abundance_solar": 0.3, "conversion_template_nh_1e22": 0.056, "absorbed_0p4_1p25_over_intrinsic_0p5_2p0": band_ratio, }, "locatelli2024_preferred_model": { "name": "combined+swcx", "geometry": "spherical beta plus exponential disk", "parameters": preferred_kwargs, "published_table1_marginal_errors": marginal_errors, "parameter_effects_on_side_balanced_fluxunit": parameter_deltas, "summary": summary, "figure3_interval_definition": ( "central-model min--max across the 14 actual M31 fields" ), "parameter_error_caveat": ( "first-order independent propagation of published marginal errors; " "the posterior covariance was not published and the pale envelope is " "not a formal joint credible interval" ), }, "locatelli2024_all_variant_envelope": variant_envelope, "locatelli2024_variants": model_results, "software": { "python_sherpa": sherpa.__version__, "astropy": astropy.__version__, "numpy": np.__version__, "pandas": pd.__version__, "scipy": scipy.__version__, }, "sources": { "locatelli2024_doi": "10.1051/0004-6361/202347061", "schellenberger2026_doi": "10.1051/0004-6361/202556835", "pan2024_xleap1_doi": "10.3847/1538-4365/ad2ea0", "qu2024_xleap2_doi": "10.3847/1538-4357/ad31a0", "zhang2024_doi": "10.1051/0004-6361/202449412", "grayson2025_doi": "10.3847/1538-4357/ae100f", "henley_shelton2013_doi": "10.1088/0004-637X/773/2/92", }, } OUTPUT.write_text(json.dumps(payload, indent=2, sort_keys=True) + "\n") print(json.dumps(summary, sort_keys=True)) print(OUTPUT) if __name__ == "__main__": main()