#!/usr/bin/env python3 """Generate reproducible MW/M31 profile-explorer tables and figures. The Milky Way calculation reproduces the Locatelli et al. (2024) reference Combined beta=0.5 geometry used by the v19 conditional-prior audit. Halo and disk emission measures are integrated separately and added without a density cross term. Spectral conversion is delegated to the authoritative Sherpa audit; continuous fixed-NH curves reuse one conversion factor. """ from __future__ import annotations import argparse import hashlib import importlib.util import json import math from pathlib import Path from types import ModuleType from typing import Any import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np import pandas as pd import scienceplots # noqa: F401 from astropy.coordinates import SkyCoord import astropy.units as u from scipy.integrate import trapezoid HERE = Path(__file__).resolve().parent ROOT = HERE.parent PAPER_DIR = ROOT / "paper_apj_v19_cgmsum_draft" AUDIT_SCRIPT = PAPER_DIR / "audit_apj_v19_cgmsum_conditional_priors.py" PUBLIC_CSV = PAPER_DIR / "m31_cgmsum_v19_primary_measurements_public.csv" LEDGER_CSV = PAPER_DIR / "m31_cgmsum_conditional_prior_ledger.csv" FROZEN_REFERENCE_CSV = ( PAPER_DIR / "m31_cgmsum_locatelli2024_reference_beta0p5_m31_footprint_predictions.csv" ) AUDIT_JSON = PAPER_DIR / "m31_cgmsum_conditional_prior_audit.json" M31_RA_DEG = 10.68470833 M31_DEC_DEG = 41.26875 M31_DISTANCE_KPC = 780.0 NOMINAL_APERTURE_RADIUS_ARCMIN = 15.0 NOMINAL_APERTURE_RADIUS_KPC = M31_DISTANCE_KPC * math.tan( math.radians(NOMINAL_APERTURE_RADIUS_ARCMIN / 60.0) ) SOLAR_RADIUS_KPC = 8.2 FIXED_NH_1E22 = 0.056 KPC_CM = 3.0856775814913673e21 ARCMIN2_PER_SR = (180.0 * 60.0 / np.pi) ** 2 PRIMARY_FLUX_UNIT = 1.0e-15 REFERENCE_BETA = 0.5 REFERENCE_C_NORM = 0.046 REFERENCE_N0_CM3 = 0.032 REFERENCE_RADIAL_SCALE_KPC = 6.2 REFERENCE_HEIGHT_SCALE_KPC = 1.1 REFERENCE_KT_KEV = 0.15 REFERENCE_ABUNDANCE_SOLAR = 0.1 LOCATELLI_O8_EMISSIVITY_PH_CM3_S = 1.649e-16 O8_RESPONSE_CLOSURE = 0.9498760469603238 ZHANG_S20_PRIMARY = 0.9418332370377434 ZHANG_RC_KPC = 7.39 ZHANG_BETA = 0.37 PROFILE_COLUMNS = ( "em_halo_kpc_cm-6", "em_disk_kpc_cm-6", "em_total_kpc_cm-6", "o8_halo_lu", "o8_disk_lu", "o8_total_lu", "direct_halo_absorbed_0p4_1p25_primary_fluxunit", "direct_disk_absorbed_0p4_1p25_primary_fluxunit", "direct_total_absorbed_0p4_1p25_primary_fluxunit", "response_matched_halo_absorbed_0p4_1p25_primary_fluxunit", "response_matched_disk_absorbed_0p4_1p25_primary_fluxunit", "response_matched_total_absorbed_0p4_1p25_primary_fluxunit", ) FIGURE_STEMS = ( "mw_latitude_profile", "mw_signed_axis_profile", "cgmsum_m31_radius_profile", "figure3_field_conversion", ) DATA_FILENAMES = ( "mw_latitude_profile.csv", "mw_m31_gc_signed_axis_profile.csv", "m31_field_profile_decomposition.csv", "conditional_m31_templates.csv", "explorer_summary.json", ) _AUDIT_MODULE: ModuleType | None = None def load_authoritative_audit() -> ModuleType: """Load the frozen audit module without modifying its source directory.""" global _AUDIT_MODULE if _AUDIT_MODULE is None: spec = importlib.util.spec_from_file_location( "m31_conditional_prior_audit", AUDIT_SCRIPT ) if spec is None or spec.loader is None: raise RuntimeError(f"Could not import authoritative audit: {AUDIT_SCRIPT}") module = importlib.util.module_from_spec(spec) spec.loader.exec_module(module) _AUDIT_MODULE = module return _AUDIT_MODULE def sha256(path: Path) -> str: digest = hashlib.sha256() with path.open("rb") as handle: for chunk in iter(lambda: handle.read(1024 * 1024), b""): digest.update(chunk) return digest.hexdigest() def m31_coordinate() -> SkyCoord: return SkyCoord(M31_RA_DEG * u.deg, M31_DEC_DEG * u.deg, frame="icrs") def galactic_center_coordinate() -> SkyCoord: return SkyCoord(l=0.0 * u.deg, b=0.0 * u.deg, frame="galactic").icrs def m31_to_gc_position_angle() -> u.Quantity: """Astronomical PA (north through east) from M31 toward Galactic center.""" return m31_coordinate().position_angle(galactic_center_coordinate()) def emission_measure_components_kpc_cm6( longitude_rad: float, latitude_rad: float, ) -> tuple[float, float, float]: """Return halo, disk, and sum of squared-density LOS integrals. The total is integral(n_h**2 + n_d**2) ds. There is deliberately no 2*n_h*n_d cross term because these are additive emission components. """ distance_kpc = np.linspace(0.001, 350.0, 35001) cos_b = np.cos(latitude_rad) cylindrical_radius = np.sqrt( SOLAR_RADIUS_KPC**2 + (distance_kpc * cos_b) ** 2 - 2.0 * SOLAR_RADIUS_KPC * distance_kpc * cos_b * np.cos(longitude_rad) ) height = distance_kpc * np.sin(latitude_rad) spherical_radius = np.sqrt(cylindrical_radius**2 + height**2) halo_density = REFERENCE_C_NORM * spherical_radius ** (-3.0 * REFERENCE_BETA) disk_density = ( REFERENCE_N0_CM3 * np.exp(-cylindrical_radius / REFERENCE_RADIAL_SCALE_KPC) * np.exp(-np.abs(height) / REFERENCE_HEIGHT_SCALE_KPC) ) halo_em = float(trapezoid(halo_density**2, distance_kpc)) disk_em = float(trapezoid(disk_density**2, distance_kpc)) return halo_em, disk_em, halo_em + disk_em def spectral_constants() -> dict[str, float]: """Evaluate fixed-NH Sherpa conversion constants once per process.""" audit = load_authoritative_audit() audit.set_xsabund("angr") audit.set_xsxsect("vern") energy_flux_per_norm = audit.energy_flux_per_norm( REFERENCE_KT_KEV, REFERENCE_ABUNDANCE_SOLAR, FIXED_NH_1E22, 0.4, 1.25, ) direct_per_em = ( energy_flux_per_norm * 1.0e-14 * KPC_CM / (4.0 * np.pi) / ARCMIN2_PER_SR / PRIMARY_FLUX_UNIT ) o8_lu_per_em = KPC_CM * LOCATELLI_O8_EMISSIVITY_PH_CM3_S / (4.0 * np.pi) return { "absorbed_energy_flux_per_apec_norm_erg_cm-2_s-1": float(energy_flux_per_norm), "direct_primary_fluxunit_per_em_kpc_cm-6": float(direct_per_em), "response_matched_primary_fluxunit_per_em_kpc_cm-6": float( direct_per_em * O8_RESPONSE_CLOSURE ), "o8_lu_per_em_kpc_cm-6": float(o8_lu_per_em), } def add_profile_quantities( table: pd.DataFrame, constants: dict[str, float], ) -> pd.DataFrame: output = table.copy() for component in ("halo", "disk", "total"): emission_measure = output[f"em_{component}_kpc_cm-6"].to_numpy() output[f"o8_{component}_lu"] = ( emission_measure * constants["o8_lu_per_em_kpc_cm-6"] ) output[f"direct_{component}_absorbed_0p4_1p25_primary_fluxunit"] = ( emission_measure * constants["direct_primary_fluxunit_per_em_kpc_cm-6"] ) output[f"response_matched_{component}_absorbed_0p4_1p25_primary_fluxunit"] = ( emission_measure * constants["response_matched_primary_fluxunit_per_em_kpc_cm-6"] ) return output def build_latitude_profile(constants: dict[str, float]) -> pd.DataFrame: m31_galactic = m31_coordinate().galactic rows: list[dict[str, float]] = [] for abs_b_deg in np.linspace(5.0, 90.0, 171): b_deg = -float(abs_b_deg) em_halo, em_disk, em_total = emission_measure_components_kpc_cm6( m31_galactic.l.rad, np.deg2rad(b_deg), ) icrs = SkyCoord( l=m31_galactic.l, b=b_deg * u.deg, frame="galactic", ).icrs rows.append( { "galactic_l_deg": float(m31_galactic.l.deg), "galactic_b_deg": b_deg, "abs_galactic_b_deg": float(abs_b_deg), "ra_deg": float(icrs.ra.deg), "dec_deg": float(icrs.dec.deg), "nh_1e22_cm-2": FIXED_NH_1E22, "em_halo_kpc_cm-6": em_halo, "em_disk_kpc_cm-6": em_disk, "em_total_kpc_cm-6": em_total, } ) return add_profile_quantities(pd.DataFrame(rows), constants) def coordinate_on_m31_gc_axis(signed_radius_kpc: float) -> tuple[SkyCoord, float]: """Return exact great-circle point; positive radius points toward the GC.""" center = m31_coordinate() angular_offset_rad = math.atan(abs(signed_radius_kpc) / M31_DISTANCE_KPC) pa = m31_to_gc_position_angle() if signed_radius_kpc < 0.0: pa = pa + 180.0 * u.deg point = center.directional_offset_by(pa, angular_offset_rad * u.rad) return point, math.copysign(angular_offset_rad, signed_radius_kpc) def build_axis_profile(constants: dict[str, float]) -> pd.DataFrame: rows: list[dict[str, float]] = [] for signed_radius_kpc in np.linspace(-50.0, 50.0, 201): point, signed_theta_rad = coordinate_on_m31_gc_axis(float(signed_radius_kpc)) galactic = point.galactic em_halo, em_disk, em_total = emission_measure_components_kpc_cm6( galactic.l.rad, galactic.b.rad, ) rows.append( { "signed_axis_kpc": float(signed_radius_kpc), "angular_offset_deg": float(np.rad2deg(signed_theta_rad)), "ra_deg": float(point.ra.deg), "dec_deg": float(point.dec.deg), "galactic_l_deg": float(galactic.l.deg), "galactic_b_deg": float(galactic.b.deg), "nh_1e22_cm-2": FIXED_NH_1E22, "axis_positive_direction": "toward Galactic center", "em_halo_kpc_cm-6": em_halo, "em_disk_kpc_cm-6": em_disk, "em_total_kpc_cm-6": em_total, } ) return add_profile_quantities(pd.DataFrame(rows), constants) def gnomonic_projection_kpc(coordinates: SkyCoord) -> tuple[np.ndarray, np.ndarray]: """Exact TAN/gnomonic coordinates along and perpendicular to the GC axis.""" center = m31_coordinate() ra = coordinates.icrs.ra.rad dec = coordinates.icrs.dec.rad ra0 = center.ra.rad dec0 = center.dec.rad delta_ra = ra - ra0 denominator = np.sin(dec0) * np.sin(dec) + np.cos(dec0) * np.cos(dec) * np.cos( delta_ra ) if np.any(denominator <= 0.0): raise ValueError("Gnomonic projection is undefined beyond 90 degrees") east = np.cos(dec) * np.sin(delta_ra) / denominator north = ( np.cos(dec0) * np.sin(dec) - np.sin(dec0) * np.cos(dec) * np.cos(delta_ra) ) / denominator pa = m31_to_gc_position_angle().to_value(u.rad) signed = M31_DISTANCE_KPC * (east * np.sin(pa) + north * np.cos(pa)) perpendicular = M31_DISTANCE_KPC * (east * np.cos(pa) - north * np.sin(pa)) return np.asarray(signed), np.asarray(perpendicular) def _assert_source_aliases_unchanged( source: pd.DataFrame, merged: pd.DataFrame ) -> None: aliases = { "cgmsum_kt_keV": "mwhalo_kt_keV", "cgmsum_kt_err_keV": "mwhalo_kt_err_keV", "cgmsum_norm": "mwhalo_norm", "cgmsum_norm_err": "mwhalo_norm_err", "cgmsum_cov_kt_norm": "mwhalo_cov_kt_norm", } for column in source.columns: if not source[column].equals(merged[column]): raise RuntimeError( f"Observed source column was changed during merge: {column}" ) for alias, historical in aliases.items(): if not np.allclose(source[alias], source[historical], equal_nan=True): raise RuntimeError( f"Public alias no longer matches historical field: {alias}" ) def build_field_table(constants: dict[str, float]) -> pd.DataFrame: public = pd.read_csv(PUBLIC_CSV, dtype={"obsid": str}) frozen = pd.read_csv(FROZEN_REFERENCE_CSV, dtype={"obsid": str}) if set(public["obsid"]) != set(frozen["obsid"]): raise RuntimeError("Public and frozen reference products have different fields") coordinates = SkyCoord( public["ra_deg"].to_numpy() * u.deg, public["dec_deg"].to_numpy() * u.deg, frame="icrs", ) signed, perpendicular = gnomonic_projection_kpc(coordinates) merged = public.copy() merged["gc_axis_signed_gnomonic_kpc"] = signed merged["gc_axis_perpendicular_gnomonic_kpc"] = perpendicular merged["gnomonic_radius_kpc"] = np.hypot(signed, perpendicular) merged["gc_axis_positive_direction"] = "toward Galactic center" # Stable explorer aliases are appended; the frozen public columns above # retain their names, values, order, and historical CGMsum aliases. merged["signed_axis_r_kpc"] = signed merged["signed_perpendicular_kpc"] = perpendicular merged["cgmsum_absorbed_0p4_1p25_primary_fluxunit"] = ( merged["absorbed_flux_soft_0p4_1p25_erg_cm-2_s-1_arcmin-2"] / PRIMARY_FLUX_UNIT ) merged["cgmsum_absorbed_0p4_1p25_staterr_primary_fluxunit"] = ( merged["absorbed_flux_soft_0p4_1p25_staterr_erg_cm-2_s-1_arcmin-2"] / PRIMARY_FLUX_UNIT ) merged["nominal_aperture_radius_arcmin"] = NOMINAL_APERTURE_RADIUS_ARCMIN merged["nominal_aperture_radial_extent_kpc"] = NOMINAL_APERTURE_RADIUS_KPC merged["rproj_aperture_min_kpc"] = np.maximum( 0.0, merged["rproj_kpc"] - NOMINAL_APERTURE_RADIUS_KPC ) merged["rproj_aperture_max_kpc"] = merged["rproj_kpc"] + NOMINAL_APERTURE_RADIUS_KPC frozen_extra = frozen.drop( columns=["side", "ra_deg", "dec_deg", "nh_hi4pi_1e22_cm-2"] ) merged = merged.merge(frozen_extra, on="obsid", how="left", validate="one_to_one") # Recompute the actual-NH branch through the authoritative Sherpa audit and # fail closed if it no longer agrees with the frozen audit CSV. audit = load_authoritative_audit() audit.set_xsabund("angr") audit.set_xsxsect("vern") audit_input = merged.copy() audit_input["l_rad"] = np.deg2rad(audit_input["galactic_l_deg"]) audit_input["b_rad"] = np.deg2rad(audit_input["galactic_b_deg"]) recomputed_direct = audit.predict_locatelli_fluxes( audit_input, beta=REFERENCE_BETA, c_norm=REFERENCE_C_NORM, n0=REFERENCE_N0_CM3, radial_scale_kpc=REFERENCE_RADIAL_SCALE_KPC, height_scale_kpc=REFERENCE_HEIGHT_SCALE_KPC, kT_keV=REFERENCE_KT_KEV, abundance_solar=REFERENCE_ABUNDANCE_SOLAR, ) frozen_direct = merged[ "direct_em_target_apec_absorbed_0p4_1p25_primary_fluxunit" ].to_numpy() frozen_matched = merged[ "reference_response_matched_absorbed_0p4_1p25_primary_fluxunit" ].to_numpy() if not np.allclose(recomputed_direct, frozen_direct, rtol=2.0e-12, atol=0.0): raise RuntimeError("Authoritative actual-NH direct predictions changed") if not np.allclose( recomputed_direct * O8_RESPONSE_CLOSURE, frozen_matched, rtol=2.0e-12, atol=0.0, ): raise RuntimeError("Authoritative response-matched predictions changed") merged["actual_nh_recomputed_direct_primary_fluxunit"] = recomputed_direct merged["actual_nh_recomputed_response_matched_primary_fluxunit"] = ( recomputed_direct * O8_RESPONSE_CLOSURE ) # Geometry-only fixed-NH branch: retain exactly the same field geometry/EM. emission_measure = merged["reference_emission_measure_n2_kpc_cm-6"].to_numpy() merged["fixed_nh_1e22_cm-2"] = FIXED_NH_1E22 merged["fixed_nh_geometry_direct_primary_fluxunit"] = ( emission_measure * constants["direct_primary_fluxunit_per_em_kpc_cm-6"] ) merged["fixed_nh_geometry_response_matched_primary_fluxunit"] = ( emission_measure * constants["response_matched_primary_fluxunit_per_em_kpc_cm-6"] ) merged["fixed_nh_geometry_o8_lu"] = ( emission_measure * constants["o8_lu_per_em_kpc_cm-6"] ) merged["mw_reference_fixed_nh_primary_fluxunit"] = merged[ "fixed_nh_geometry_response_matched_primary_fluxunit" ] merged["mw_reference_actual_nh_primary_fluxunit"] = merged[ "reference_response_matched_absorbed_0p4_1p25_primary_fluxunit" ] _assert_source_aliases_unchanged(public, merged) return merged def zhang_primary_profile(radius_kpc: np.ndarray | float) -> np.ndarray: radius = np.asarray(radius_kpc, dtype=float) exponent = -3.0 * ZHANG_BETA + 0.5 shape = (1.0 + (radius / ZHANG_RC_KPC) ** 2) ** exponent shape_at_20 = (1.0 + (20.0 / ZHANG_RC_KPC) ** 2) ** exponent return ZHANG_S20_PRIMARY * shape / shape_at_20 def build_conditional_templates() -> pd.DataFrame: ledger = pd.read_csv(LEDGER_CSV) rows: list[dict[str, Any]] = [] for radius, value in zip( np.linspace(10.0, 30.0, 81), zhang_primary_profile(np.linspace(10.0, 30.0, 81)), strict=True, ): rows.append( { "template_id": "zhang2024_m31_mass_beta_profile", "label": "Zhang et al. 2024 M31-mass stack", "representation": "sampled_beta_profile_10_to_30_kpc", "radius_kpc": float(radius), "bin_low_kpc": np.nan, "bin_high_kpc": np.nan, "primary_central": float(value), "primary_low": float(value), "primary_high": float(value), } ) grayson = ledger.loc[ledger["prior_id"].str.startswith("grayson2025_")] for _, row in grayson.iterrows(): rows.append( { "template_id": row["prior_id"], "label": row["label"], "representation": "single_10_to_30_kpc_bin_band_no_interpolation", "radius_kpc": np.nan, "bin_low_kpc": 10.0, "bin_high_kpc": 30.0, "primary_central": float(row["primary_central"]), "primary_low": float(row["primary_low"]), "primary_high": float(row["primary_high"]), } ) return pd.DataFrame(rows) def weighted_mean(values: np.ndarray, errors: np.ndarray) -> float: weights = 1.0 / np.square(errors) return float(np.sum(weights * values) / np.sum(weights)) def field_estimators( table: pd.DataFrame, value_column: str, error_column: str, ) -> dict[str, float]: values = table[value_column].to_numpy(dtype=float) errors = table[error_column].to_numpy(dtype=float) side_values: dict[str, float] = {} for side in ("North/NW", "South/SE"): mask = table["side"].eq(side).to_numpy() side_values[side] = weighted_mean(values[mask], errors[mask]) return { "north_inverse_variance": side_values["North/NW"], "south_inverse_variance": side_values["South/SE"], "side_balanced": 0.5 * (side_values["North/NW"] + side_values["South/SE"]), "all_inverse_variance": weighted_mean(values, errors), "unweighted_mean": float(np.mean(values)), "median": float(np.median(values)), "p16": float(np.quantile(values, 0.16)), "p84": float(np.quantile(values, 0.84)), "minimum": float(np.min(values)), "maximum": float(np.max(values)), } def _save_figure(fig: plt.Figure, output_dir: Path, stem: str, column: str) -> None: png_path = output_dir / f"{stem}_{column}.png" pdf_path = output_dir / f"{stem}_{column}.pdf" fig.savefig( png_path, dpi=400, metadata={"Software": "generate_explorer_products.py"}, ) fig.savefig( pdf_path, metadata={ "Creator": "generate_explorer_products.py", "CreationDate": None, "ModDate": None, }, ) plt.close(fig) def _apply_column_style(column: str) -> None: base = 7.0 if column == "1col" else 8.5 plt.rcParams.update( { "font.size": base, "axes.titlesize": base + 1.0, "axes.labelsize": base, "xtick.labelsize": base - 0.5, "ytick.labelsize": base - 0.5, "legend.fontsize": base - 1.0, } ) def plot_latitude_profile(table: pd.DataFrame, output_dir: Path, column: str) -> None: _apply_column_style(column) width = 3.35 if column == "1col" else 7.0 fig, axes = plt.subplots( 2, 1, figsize=(width, width * 0.9), sharex=True, layout="constrained" ) x = table["abs_galactic_b_deg"] axes[0].plot(x, table["em_halo_kpc_cm-6"], label="Halo $n_h^2$", color="#3b6fb6") axes[0].plot(x, table["em_disk_kpc_cm-6"], label="Disk $n_d^2$", color="#d9792b") axes[0].plot(x, table["em_total_kpc_cm-6"], label="Total", color="black") axes[0].set_ylabel(r"EM (kpc cm$^{-6}$)") axes[0].legend(fontsize=6 if column == "1col" else 7, ncol=3) axes[1].plot( x, table["direct_total_absorbed_0p4_1p25_primary_fluxunit"], label="Direct EM mapping", color="#777777", linestyle="--", ) axes[1].plot( x, table["response_matched_total_absorbed_0p4_1p25_primary_fluxunit"], label="O VIII response-matched", color="#6a3d9a", ) m31_abs_b = abs(m31_coordinate().galactic.b.deg) axes[1].axvline( m31_abs_b, color="#2a9d8f", linewidth=0.9, linestyle=":", label="M31" ) axes[1].set_xlabel(r"$|b|$ (deg), fixed $l=l_{\rm M31}$") axes[1].set_ylabel(r"Flux ($10^{-15}$ units)") axes[1].legend(fontsize=6 if column == "1col" else 7) fig.suptitle( "Locatelli reference MW latitude profile", fontsize=8 if column == "1col" else 10, ) _save_figure(fig, output_dir, "mw_latitude_profile", column) def plot_axis_profile( table: pd.DataFrame, fields: pd.DataFrame, output_dir: Path, column: str, ) -> None: _apply_column_style(column) width = 3.35 if column == "1col" else 7.0 fig, axes = plt.subplots( 2, 1, figsize=(width, width * 0.88), sharex=True, layout="constrained" ) x = table["signed_axis_kpc"] axes[0].plot( x, table["response_matched_halo_absorbed_0p4_1p25_primary_fluxunit"], label="Halo", color="#3b6fb6", ) axes[0].plot( x, table["response_matched_disk_absorbed_0p4_1p25_primary_fluxunit"], label="Disk", color="#d9792b", ) axes[0].plot( x, table["response_matched_total_absorbed_0p4_1p25_primary_fluxunit"], label="Total", color="black", ) axes[0].scatter( fields["gc_axis_signed_gnomonic_kpc"], fields["reference_response_matched_absorbed_0p4_1p25_primary_fluxunit"], s=13 if column == "1col" else 20, marker="o", facecolors="none", edgecolors="#6a3d9a", linewidths=0.8, label="Field centers: actual HI4PI NH", zorder=4, ) axes[0].scatter( fields["gc_axis_signed_gnomonic_kpc"], fields["fixed_nh_geometry_response_matched_primary_fluxunit"], s=10 if column == "1col" else 16, marker="x", color="#2a9d8f", linewidths=0.8, label="Field centers: fixed NH", zorder=4, ) axes[0].set_ylabel(r"Matched flux ($10^{-15}$)") axes[0].legend(fontsize=4.7 if column == "1col" else 6.5, ncol=2) center = table.loc[ np.isclose(x, 0.0), "response_matched_total_absorbed_0p4_1p25_primary_fluxunit" ].iloc[0] axes[1].plot( x, table["response_matched_total_absorbed_0p4_1p25_primary_fluxunit"] / center, color="#6a3d9a", ) axes[1].axhline(1.0, color="0.5", linewidth=0.8, linestyle="--") axes[1].set_xlabel("Signed tangent radius (kpc)\n(+ toward Galactic center)") axes[1].set_ylabel("MW / center") fig.suptitle( "Great-circle M31--Galactic-center axis", fontsize=8 if column == "1col" else 10 ) _save_figure(fig, output_dir, "mw_signed_axis_profile", column) def plot_conditional_templates( fields: pd.DataFrame, templates: pd.DataFrame, output_dir: Path, column: str, ) -> None: _apply_column_style(column) width = 3.35 if column == "1col" else 7.0 fig, ax = plt.subplots(figsize=(width, width * 0.78), layout="constrained") x = fields["rproj_kpc"] observed = ( fields["absorbed_flux_soft_0p4_1p25_erg_cm-2_s-1_arcmin-2"] / PRIMARY_FLUX_UNIT ) errors = ( fields["absorbed_flux_soft_0p4_1p25_staterr_erg_cm-2_s-1_arcmin-2"] / PRIMARY_FLUX_UNIT ) for side, color, marker in ( ("North/NW", "#377eb8", "o"), ("South/SE", "#e41a1c", "s"), ): selected = fields["side"].eq(side) ax.errorbar( x[selected], observed[selected], xerr=NOMINAL_APERTURE_RADIUS_KPC, yerr=errors[selected], fmt=marker, ms=3, color=color, ecolor=color, elinewidth=0.7, capsize=1.5, linestyle="none", label=f"Observed CGMsum: {side}", zorder=5, ) zhang = templates.loc[ templates["template_id"].eq("zhang2024_m31_mass_beta_profile") ] radius = zhang["radius_kpc"].to_numpy() values = zhang["primary_central"].to_numpy() ax.plot(radius, values, color="#d95f02", label="Zhang 2024 beta profile") colors = ["#1b9e77", "#377eb8", "#4daf4a", "#984ea3", "#e41a1c"] grayson = templates.loc[templates["representation"].str.startswith("single_")] for color, (_, row) in zip(colors, grayson.iterrows(), strict=True): left, right = row["bin_low_kpc"], row["bin_high_kpc"] ax.fill_between( [left, right], row["primary_low"], row["primary_high"], color=color, alpha=0.16, linewidth=0, ) ax.hlines(row["primary_central"], left, right, color=color, linewidth=1.0) ax.plot( [], [], color=color, label=row["label"].replace("Grayson et al. 2025 ", "") ) ax.set_yscale("log") ax.set_xlim(5, 35) ax.set_ylim(0.16, 25.0) ax.set_xlabel("Circular projected radius from M31 (kpc)") ax.set_ylabel(r"Surface brightness ($10^{-15}$)") ax.set_title( "Observed CGMsum and conditional M31 templates", fontsize=8 if column == "1col" else 10, ) ax.legend(fontsize=4.7 if column == "1col" else 6.8, ncol=2) _save_figure(fig, output_dir, "cgmsum_m31_radius_profile", column) def plot_field_decomposition( fields: pd.DataFrame, output_dir: Path, column: str ) -> None: _apply_column_style(column) width = 3.35 if column == "1col" else 7.0 fig, axes = plt.subplots( 2, 1, figsize=(width, width * 0.95), sharex=True, layout="constrained" ) ordered = fields.sort_values("gc_axis_signed_gnomonic_kpc") x = ordered["gc_axis_signed_gnomonic_kpc"] axes[0].plot( x, ordered["direct_em_target_apec_absorbed_0p4_1p25_primary_fluxunit"], "o", ms=3, label="Actual NH: direct EM", color="#777777", ) axes[0].plot( x, ordered["reference_response_matched_absorbed_0p4_1p25_primary_fluxunit"], "o", ms=3, label="Actual NH: response-matched", color="#6a3d9a", ) axes[0].plot( x, ordered["fixed_nh_geometry_response_matched_primary_fluxunit"], "s", ms=2.6, label="Fixed NH: geometry only", color="#2a9d8f", ) axes[0].set_ylabel(r"Predicted flux ($10^{-15}$)") axes[0].legend(fontsize=5.5 if column == "1col" else 7) axes[1].plot( x, ordered["o8_response_matched_emissivity_closure_scale"], "o", ms=3, color="#d95f02", label="O VIII closure", ) transmission_effect = ( ordered["reference_response_matched_absorbed_0p4_1p25_primary_fluxunit"] / ordered["fixed_nh_geometry_response_matched_primary_fluxunit"] ) axes[1].plot( x, transmission_effect, "s", ms=3, color="#377eb8", label="Actual/fixed NH", ) axes[1].axhline(1.0, color="0.5", linewidth=0.8, linestyle=":") axes[1].set_xlabel("Signed gnomonic radius (kpc)\n(+ toward Galactic center)") axes[1].set_ylabel("Factor") axes[1].legend(fontsize=5.5 if column == "1col" else 7) fig.suptitle( "Figure 3 field-conversion decomposition", fontsize=8 if column == "1col" else 10, ) _save_figure(fig, output_dir, "figure3_field_conversion", column) def expected_output_paths(output_dir: Path = HERE) -> list[Path]: product_dir = output_dir / "products" figure_dir = output_dir / "figures" paths = [product_dir / name for name in DATA_FILENAMES] for stem in FIGURE_STEMS: for column in ("1col", "2col"): for suffix in ("png", "pdf"): paths.append(figure_dir / f"{stem}_{column}.{suffix}") return paths def build_summary( constants: dict[str, float], latitude: pd.DataFrame, axis: pd.DataFrame, fields: pd.DataFrame, ) -> dict[str, Any]: error_column = "absorbed_flux_soft_0p4_1p25_staterr_erg_cm-2_s-1_arcmin-2" fields_for_observed = fields.assign( observed_primary_fluxunit=( fields["absorbed_flux_soft_0p4_1p25_erg_cm-2_s-1_arcmin-2"] / PRIMARY_FLUX_UNIT ), observed_primary_error_fluxunit=(fields[error_column] / PRIMARY_FLUX_UNIT), ) estimators = { "observed_cgmsum": field_estimators( fields_for_observed, "observed_primary_fluxunit", "observed_primary_error_fluxunit", ), "actual_nh_direct_em": field_estimators( fields, "direct_em_target_apec_absorbed_0p4_1p25_primary_fluxunit", error_column, ), "actual_nh_response_matched": field_estimators( fields, "reference_response_matched_absorbed_0p4_1p25_primary_fluxunit", error_column, ), "fixed_nh_geometry_only_direct": field_estimators( fields, "fixed_nh_geometry_direct_primary_fluxunit", error_column, ), "fixed_nh_geometry_only_response_matched": field_estimators( fields, "fixed_nh_geometry_response_matched_primary_fluxunit", error_column, ), } center = axis.loc[np.isclose(axis["signed_axis_kpc"], 0.0)].iloc[0] plus_30 = axis.loc[np.isclose(axis["signed_axis_kpc"], 30.0)].iloc[0] return { "schema": "standalone-mw-m31-profile-explorer-v1", "formula": { "density_components": { "halo": "n_h = c * r^(-3 beta)", "disk": "n_d = n0 exp(-R/Rd) exp(-|z|/zd)", "emission_measure": "integral (n_h^2 + n_d^2) ds; no cross term", }, "gnomonic": ( "exact TAN projection about M31; signed coordinate rotated to the " "great-circle direction from M31 toward Galactic center" ), "axis_radius": "R = D tan(theta)", "zhang_surface_brightness": ( "S(R)=S20*[(1+(R/rc)^2)/(1+(20/rc)^2)]^(-3 beta+1/2)" ), }, "constants": { "m31_ra_deg": M31_RA_DEG, "m31_dec_deg": M31_DEC_DEG, "m31_distance_kpc": M31_DISTANCE_KPC, "nominal_aperture_radius_arcmin": NOMINAL_APERTURE_RADIUS_ARCMIN, "nominal_aperture_radial_extent_kpc": NOMINAL_APERTURE_RADIUS_KPC, "m31_galactic_l_deg": float(m31_coordinate().galactic.l.deg), "m31_galactic_b_deg": float(m31_coordinate().galactic.b.deg), "m31_to_gc_position_angle_deg": float(m31_to_gc_position_angle().deg), "solar_radius_kpc": SOLAR_RADIUS_KPC, "fixed_nh_1e22_cm-2": FIXED_NH_1E22, "reference_beta": REFERENCE_BETA, "reference_c_norm": REFERENCE_C_NORM, "reference_n0_cm-3": REFERENCE_N0_CM3, "reference_radial_scale_kpc": REFERENCE_RADIAL_SCALE_KPC, "reference_height_scale_kpc": REFERENCE_HEIGHT_SCALE_KPC, "reference_kT_keV": REFERENCE_KT_KEV, "reference_abundance_solar": REFERENCE_ABUNDANCE_SOLAR, "o8_emissivity_ph_cm3_s": LOCATELLI_O8_EMISSIVITY_PH_CM3_S, "o8_response_matched_closure": O8_RESPONSE_CLOSURE, "zhang_s20_primary_fluxunit": ZHANG_S20_PRIMARY, "zhang_rc_kpc": ZHANG_RC_KPC, "zhang_beta": ZHANG_BETA, **constants, }, "profiles": { "latitude_rows": int(len(latitude)), "axis_rows": int(len(axis)), "axis_positive_direction": "toward Galactic center", "axis_center_response_matched_primary_fluxunit": float( center["response_matched_total_absorbed_0p4_1p25_primary_fluxunit"] ), "axis_plus30_over_center": float( plus_30["response_matched_total_absorbed_0p4_1p25_primary_fluxunit"] / center["response_matched_total_absorbed_0p4_1p25_primary_fluxunit"] ), }, "field_count": int(len(fields)), "estimators_primary_fluxunit": estimators, "provenance": { path.name: {"absolute_path": str(path), "sha256": sha256(path)} for path in ( AUDIT_SCRIPT, PUBLIC_CSV, LEDGER_CSV, FROZEN_REFERENCE_CSV, AUDIT_JSON, ) }, "expected_products": [ *(f"products/{name}" for name in DATA_FILENAMES), *( f"figures/{stem}_{column}.{suffix}" for stem in FIGURE_STEMS for column in ("1col", "2col") for suffix in ("png", "pdf") ), ], } def generate(output_dir: Path = HERE) -> list[Path]: product_dir = output_dir / "products" figure_dir = output_dir / "figures" product_dir.mkdir(parents=True, exist_ok=True) figure_dir.mkdir(parents=True, exist_ok=True) plt.style.use(["science", "no-latex"]) plt.rcParams["savefig.bbox"] = None constants = spectral_constants() latitude = build_latitude_profile(constants) axis = build_axis_profile(constants) fields = build_field_table(constants) templates = build_conditional_templates() latitude.to_csv(product_dir / "mw_latitude_profile.csv", index=False) axis.to_csv(product_dir / "mw_m31_gc_signed_axis_profile.csv", index=False) fields.to_csv(product_dir / "m31_field_profile_decomposition.csv", index=False) templates.to_csv(product_dir / "conditional_m31_templates.csv", index=False) for column in ("1col", "2col"): plot_latitude_profile(latitude, figure_dir, column) plot_axis_profile(axis, fields, figure_dir, column) plot_conditional_templates(fields, templates, figure_dir, column) plot_field_decomposition(fields, figure_dir, column) summary = build_summary(constants, latitude, axis, fields) (product_dir / "explorer_summary.json").write_text( json.dumps(summary, indent=2, sort_keys=True) + "\n", encoding="utf-8", ) outputs = expected_output_paths(output_dir) missing = [ str(path) for path in outputs if not path.is_file() or path.stat().st_size == 0 ] if missing: raise RuntimeError(f"Missing or empty explorer products: {missing}") return outputs def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument( "--output-dir", type=Path, default=HERE, help="Build root containing products/ and figures/ (default: script directory)", ) return parser.parse_args() def main() -> None: args = parse_args() outputs = generate(args.output_dir.resolve()) print(f"Generated {len(outputs)} products") for path in outputs: print(path.resolve()) if __name__ == "__main__": main()