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 from scipy.ndimage import gaussian_filter1d 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_reference_beta0p5_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 photon_flux_per_norm( kT_keV: float, abundance_solar: float, low_keV: float, high_keV: float, ) -> float: plasma = xsapec.audit_line_anchor plasma.kT = kT_keV plasma.Abundanc = abundance_solar plasma.redshift = 0.0 plasma.norm = 1.0 edges = np.linspace(low_keV, high_keV, 30001) return float(np.sum(plasma(edges[:-1], edges[1:]))) def response_smoothed_photon_flux_per_norm( kT_keV: float, abundance_solar: float, low_keV: float, high_keV: float, *, fwhm_keV: float, ) -> float: plasma = xsapec.audit_line_anchor_smoothed plasma.kT = kT_keV plasma.Abundanc = abundance_solar plasma.redshift = 0.0 plasma.norm = 1.0 step_keV = 1.0e-4 edges = np.arange(0.3, 1.0 + step_keV, step_keV) midpoint = 0.5 * (edges[:-1] + edges[1:]) photon_density = plasma(edges[:-1], edges[1:]) / step_keV sigma_bins = fwhm_keV / 2.354820045 / step_keV smoothed_density = gaussian_filter1d( photon_density, sigma=sigma_bins, mode="constant", ) selected = (midpoint >= low_keV) & (midpoint < high_keV) return float(np.sum(smoothed_density[selected] * step_keV)) 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 ) reference_kT_keV = 0.15 reference_abundance_solar = 0.1 locatelli_o8_band_emissivity_ph_cm3_s = ( 1.649e-15 * reference_abundance_solar ) target_o8_band_raw_emissivity_ph_cm3_s = ( photon_flux_per_norm( reference_kT_keV, reference_abundance_solar, 0.614, 0.694, ) * 1.0e-14 ) target_o8_band_emissivity_ph_cm3_s = ( response_smoothed_photon_flux_per_norm( reference_kT_keV, reference_abundance_solar, 0.614, 0.694, fwhm_keV=0.080, ) * 1.0e-14 ) o8_line_anchor_scale = ( locatelli_o8_band_emissivity_ph_cm3_s / target_o8_band_emissivity_ph_cm3_s ) 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]] = {} reference_direct_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_beta0p5": reference_direct_flux = flux if reference_direct_flux is None: raise RuntimeError("reference combined beta=0.5 model was not evaluated") reference_flux = reference_direct_flux * o8_line_anchor_scale reference_estimators = side_estimators(table, reference_flux, error_name) direct_estimators = side_estimators(table, reference_direct_flux, error_name) matched_z03_direct_flux = predict_locatelli_fluxes( table, beta=0.5, c_norm=0.046, n0=0.032, radial_scale_kpc=6.2, height_scale_kpc=1.1, kT_keV=reference_kT_keV, abundance_solar=0.3, ) matched_z03_o8_emissivity = ( response_smoothed_photon_flux_per_norm( reference_kT_keV, 0.3, 0.614, 0.694, fwhm_keV=0.080, ) * 1.0e-14 ) matched_z03_flux = matched_z03_direct_flux * ( locatelli_o8_band_emissivity_ph_cm3_s / matched_z03_o8_emissivity ) matched_z03_estimators = side_estimators(table, matched_z03_flux, error_name) reference_kwargs = { "beta": 0.5, "c_norm": 0.046, "n0": 0.032, "radial_scale_kpc": 6.2, "height_scale_kpc": 1.1, "kT_keV": reference_kT_keV, "abundance_solar": reference_abundance_solar, } marginal_errors = { "c_norm": 0.001, "n0": 0.004, "radial_scale_kpc": 0.4, "height_scale_kpc": 0.1, } parameter_deltas: dict[str, float] = {} for parameter, sigma in marginal_errors.items(): low_kwargs = reference_kwargs.copy() high_kwargs = reference_kwargs.copy() low_kwargs[parameter] -= sigma high_kwargs[parameter] += sigma low_flux = ( predict_locatelli_fluxes(table, **low_kwargs) * o8_line_anchor_scale ) high_flux = ( predict_locatelli_fluxes(table, **high_kwargs) * o8_line_anchor_scale ) 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())))) ) emission_measures = np.array( [ emission_measure_kpc_cm6( float(row["l_rad"]), float(row["b_rad"]), beta=reference_kwargs["beta"], c_norm=reference_kwargs["c_norm"], n0=reference_kwargs["n0"], radial_scale_kpc=reference_kwargs["radial_scale_kpc"], height_scale_kpc=reference_kwargs["height_scale_kpc"], ) for _, row in table.iterrows() ] ) o8_band_intensity_lu = ( emission_measures * KPC_CM * locatelli_o8_band_emissivity_ph_cm3_s / (4.0 * np.pi) ) 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["reference_emission_measure_n2_kpc_cm-6"] = emission_measures footprint["reference_intrinsic_o8_0p614_0p694_intensity_lu"] = ( o8_band_intensity_lu ) footprint["direct_em_target_apec_absorbed_0p4_1p25_primary_fluxunit"] = ( reference_direct_flux ) footprint["o8_response_matched_emissivity_closure_scale"] = o8_line_anchor_scale footprint["reference_response_matched_absorbed_0p4_1p25_primary_fluxunit"] = ( reference_flux ) footprint[ "line_normalized_z0p3_absorbed_0p4_1p25_primary_fluxunit" ] = matched_z03_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(reference_flux.min()), "field_p16": float(np.quantile(reference_flux, 0.16)), "field_median": float(np.median(reference_flux)), "field_p84": float(np.quantile(reference_flux, 0.84)), "field_high": float(reference_flux.max()), "north": reference_estimators["north_primary_fluxunit"], "south": reference_estimators["south_primary_fluxunit"], "side_balanced": reference_estimators["side_balanced_primary_fluxunit"], "all_inverse_variance": reference_estimators[ "all_inverse_variance_primary_fluxunit" ], "direct_em_mapping_field_low": float(reference_direct_flux.min()), "direct_em_mapping_field_high": float(reference_direct_flux.max()), "direct_em_mapping_side_balanced": direct_estimators[ "side_balanced_primary_fluxunit" ], "line_normalized_z0p3_field_low": float(matched_z03_flux.min()), "line_normalized_z0p3_field_high": float(matched_z03_flux.max()), "line_normalized_z0p3_side_balanced": matched_z03_estimators[ "side_balanced_primary_fluxunit" ], "line_normalized_z0p3_all_inverse_variance": matched_z03_estimators[ "all_inverse_variance_primary_fluxunit" ], "table1_independent_marginal_sigma": parameter_sigma, "table1_independent_marginal_low": reference_estimators[ "side_balanced_primary_fluxunit" ] - parameter_sigma, "table1_independent_marginal_high": reference_estimators[ "side_balanced_primary_fluxunit" ] + parameter_sigma, } ledger = pd.read_csv(LEDGER) expected = ledger.loc[ ledger["prior_id"].eq( "locatelli2024_reference_beta0p5_m31_footprint" ) ].iloc[0] for value, ledger_key in ( (summary["field_low"], "primary_low"), (summary["side_balanced"], "primary_central"), (summary["field_high"], "primary_high"), (summary["all_inverse_variance"], "all_field_inverse_variance"), (summary["table1_independent_marginal_low"], "parameter_error_low"), (summary["table1_independent_marginal_high"], "parameter_error_high"), (summary["north"], "north_primary"), (summary["south"], "south_primary"), (summary["direct_em_mapping_field_low"], "conversion_check_low"), (summary["direct_em_mapping_field_high"], "conversion_check_high"), ): 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-v3", "input_measurements": str(MEASUREMENTS.relative_to(ROOT)), "reference_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_reference_model": { "name": "combined_beta0p5_reference", "geometry": "spherical beta plus exponential disk", "parameters": reference_kwargs, "published_table1_marginal_errors": marginal_errors, "o8_line_anchor": { "locatelli_intrinsic_0p614_0p694_emissivity_ph_cm3_s": ( locatelli_o8_band_emissivity_ph_cm3_s ), "target_apec_raw_0p614_0p694_emissivity_ph_cm3_s": ( target_o8_band_raw_emissivity_ph_cm3_s ), "target_apec_fwhm80ev_smoothed_0p614_0p694_emissivity_ph_cm3_s": ( target_o8_band_emissivity_ph_cm3_s ), "target_apec_z0p3_fwhm80ev_smoothed_emissivity_ph_cm3_s": ( matched_z03_o8_emissivity ), "multiplicative_scale": o8_line_anchor_scale, "purpose": ( "close the target APEC emissivity against Locatelli's " "FWHM~80 eV-smoothed O VIII-band emissivity before converting " "the reference density model to absorbed 0.4--1.25 keV" ), }, "parameter_effects_on_side_balanced_fluxunit": parameter_deltas, "summary": summary, "figure3_interval_definition": ( "reference-model response-matched O VIII-band 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" ), "absolute_conversion_caveat": ( "Locatelli fitted an FWHM~80 eV-smoothed 0.614--0.694 keV " "O VIII-band photon intensity, not an absorbed broadband energy " "flux. The direct-EM target-APEC mapping before the response-matched " "emissivity closure is retained as a conversion check." ), }, "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", "zhang2026_baryon_budget_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()