from __future__ import annotations

from pathlib import Path

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
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
OUTPUT_DIR = HERE
PRIMARY_CSV = (
    ROOT
    / "independent_cgm_results_v19_mos14to20keV_excluded_hi4pi_nh_per_pointing"
    / "independent_cgm_v19_all22_measurements.csv"
)
PRIOR_LEDGER_CSV = HERE / "m31_cgmsum_conditional_prior_ledger.csv"
PREDICTIONS_CSV = HERE / "m31_cgmsum_total_difference_predictions.csv"

FLUX_UNIT = 1.0e-15
NORTH_FLUX = 1.1062308841694078
NORTH_ERROR = 0.05357907323160643
SOUTH_FLUX = 0.893850973242182
SOUTH_ERROR = 0.2032861967092256
ALL_FLUX = 1.0923775268463574
ALL_ERROR = 0.05186950369566604

# Locatelli et al. (2024) reference Combined (beta=0.5) model, projected
# through the 14-field M31 primary footprint.  The dark interval preserves the
# fitted intrinsic 0.614--0.694 keV O VIII-band intensity before conversion to
# the target absorbed broadband.  The pale interval is a diagonal Table-1
# marginal-error diagnostic; the published posterior covariance is absent.
LOCATELLI_M31_CENTRAL = 1.5596097515217344
LOCATELLI_M31_ALL_FIELD = 1.6017123021968118
LOCATELLI_M31_FIELD_LOW = 1.4312606908446714
LOCATELLI_M31_FIELD_HIGH = 1.6682262195476179
LOCATELLI_M31_NORTH = 1.6080149355676707
LOCATELLI_M31_SOUTH = 1.5112045674757981
LOCATELLI_M31_STAT_SIGMA = 0.3979910769938221
LOCATELLI_DIRECT_EM_LOW = 1.506786801735674
LOCATELLI_DIRECT_EM_HIGH = 1.7562567504318796

# For an APEC plasma with kT=0.175306 keV, Z=0.3 Z_sun, phabs nH=0.056,
# angr abundances, and vern cross-sections, recomputed with Sherpa 4.18/XSPEC:
# F_abs(0.4-1.25 keV) / F_intrinsic(0.5-2.0 keV).
PRIMARY_TO_INTRINSIC_05_2_RATIO = 0.7719131901134512
KPC_CM = 3.0856775814913673e21
ARCMIN2_SR = (np.pi / (180.0 * 60.0)) ** 2
LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY = (
    ARCMIN2_SR / (4.0 * np.pi * KPC_CM**2) / FLUX_UNIT * PRIMARY_TO_INTRINSIC_05_2_RATIO
)


def expected_output_paths() -> list[Path]:
    stems = (
        "m31_cgmsum_mw_m31_constraint_plane",
        "m31_cgmsum_mw_m31_constraint_plane_loglog",
        "m31_cgmsum_total_difference_predictive_plane",
    )
    return [
        Path(OUTPUT_DIR) / f"{stem}_{column}.{suffix}"
        for stem in stems
        for column in ("1col", "2col")
        for suffix in ("png", "pdf")
    ]


def measurement_total_difference() -> tuple[float, float, np.ndarray]:
    total = 0.5 * (NORTH_FLUX + SOUTH_FLUX)
    difference = NORTH_FLUX - SOUTH_FLUX
    covariance = np.array(
        [
            [
                0.25 * (NORTH_ERROR**2 + SOUTH_ERROR**2),
                0.5 * (NORTH_ERROR**2 - SOUTH_ERROR**2),
            ],
            [0.5 * (NORTH_ERROR**2 - SOUTH_ERROR**2), NORTH_ERROR**2 + SOUTH_ERROR**2],
        ]
    )
    return total, difference, covariance


def _primary_fields() -> pd.DataFrame:
    table = pd.read_csv(PRIMARY_CSV, dtype={"obsid": str})
    return table.loc[table["quality_primary"]].copy()


def _weighted_side_mean(table: pd.DataFrame, values: np.ndarray, side: str) -> float:
    error_name = "absorbed_flux_soft_0p4_1p25_staterr_erg_cm-2_s-1_arcmin-2"
    mask = table["side"].eq(side).to_numpy()
    weights = 1.0 / table.loc[mask, error_name].to_numpy() ** 2
    return float(np.sum(weights * values[mask]) / np.sum(weights))


def _contrast_slope(table: pd.DataFrame, values: np.ndarray) -> float:
    north = _weighted_side_mean(table, values, "North/NW")
    south = _weighted_side_mean(table, values, "South/SE")
    return 2.0 * (north - south) / (north + south)


def _zhang_profile(radius_kpc: np.ndarray) -> np.ndarray:
    core_radius_kpc = 7.39
    beta = 0.37
    return (1.0 + (radius_kpc / core_radius_kpc) ** 2) ** (0.5 - 3.0 * beta)


def _locatelli_emission_shape(
    l_rad: np.ndarray,
    b_rad: np.ndarray,
    *,
    beta: float,
    c_norm: float,
    n0: float,
    radial_scale_kpc: float,
    height_scale_kpc: float,
) -> np.ndarray:
    distance_kpc = np.linspace(0.001, 350.0, 35001)
    solar_radius_kpc = 8.2
    output = []
    for longitude, latitude in zip(l_rad, b_rad, strict=True):
        cylindrical_radius = np.sqrt(
            solar_radius_kpc**2
            + (distance_kpc * np.cos(latitude)) ** 2
            - 2.0
            * solar_radius_kpc
            * distance_kpc
            * np.cos(latitude)
            * np.cos(longitude)
        )
        height = distance_kpc * np.sin(latitude)
        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)
        )
        output.append(trapezoid(spherical_density**2 + disk_density**2, distance_kpc))
    return np.asarray(output)


def _apec_primary_band_transmission(nh_1e22: np.ndarray, kT_keV: float) -> np.ndarray:
    nh_grid = np.linspace(0.04, 0.08, 9)
    transmission_grids = {
        0.15: np.array(
            [
                0.729846506,
                0.701996621,
                0.675279178,
                0.649645958,
                0.625050780,
                0.601449548,
                0.578800060,
                0.557062121,
                0.536197223,
            ]
        ),
        0.178: np.array(
            [
                0.744631128,
                0.718073718,
                0.692547309,
                0.668009357,
                0.644419095,
                0.621737548,
                0.599927377,
                0.578952992,
                0.558780250,
            ]
        ),
        0.225: np.array(
            [
                0.772413462,
                0.748317246,
                0.725066336,
                0.702628298,
                0.680972008,
                0.660067629,
                0.639886507,
                0.620401284,
                0.601585631,
            ]
        ),
    }
    return np.interp(nh_1e22, nh_grid, transmission_grids[kT_keV])


def model_contrast_slopes() -> dict[str, float | tuple[float, float]]:
    table = _primary_fields()
    m31_slope = _contrast_slope(table, _zhang_profile(table["rproj_kpc"].to_numpy()))

    # Locatelli et al. (2024), Table 1. The variants bracket the spatial
    # systematics used in that paper. Only relative projected morphology is
    # needed here. The absorbed-band prediction also includes the HI4PI nH
    # variation across the footprint through APEC transmission lookup tables.
    variants = (
        ("disk", 0.5, 0.0, 0.017, 8.0, 3.3, 0.15),
        ("beta03", 0.3, 0.0009, 0.017, 8.0, 3.3, 0.15),
        ("beta05", 0.5, 0.046, 0.032, 6.2, 1.1, 0.15),
        ("beta07", 0.7, 0.176, 0.016, 9.9, 1.6, 0.15),
        ("swcx", 0.5, 0.032, 0.063, 3.9, 0.9, 0.178),
        ("instr", 0.5, 0.043, 0.055, 4.6, 0.9, 0.15),
        ("high_emissivity", 0.5, 0.0136, 0.0094, 6.2, 1.1, 0.225),
    )
    galactic = SkyCoord(
        table["ra_deg"].to_numpy() * u.deg,
        table["dec_deg"].to_numpy() * u.deg,
    ).galactic
    l_rad = galactic.l.rad
    b_rad = galactic.b.rad
    nh_1e22 = table["nh_hi4pi_1e22_cm-2"].to_numpy()
    mw_slopes: dict[str, float] = {}
    for name, beta, c_norm, n0, radial_scale, height_scale, kT_keV in variants:
        values = _locatelli_emission_shape(
            l_rad,
            b_rad,
            beta=beta,
            c_norm=c_norm,
            n0=n0,
            radial_scale_kpc=radial_scale,
            height_scale_kpc=height_scale,
        )
        values *= _apec_primary_band_transmission(nh_1e22, kT_keV)
        mw_slopes[name] = _contrast_slope(table, values)
    family = (min(mw_slopes.values()), max(mw_slopes.values()))
    return {
        "Zhang 2024 symmetric M31": m31_slope,
        "Locatelli 2024 MW reference": mw_slopes["beta05"],
        "Locatelli 2024 MW sensitivity family": family,
        "Uniform component": 0.0,
    }


def build_prior_ledger() -> pd.DataFrame:
    rows: list[dict[str, object]] = []

    hs_conversion = PRIMARY_TO_INTRINSIC_05_2_RATIO / 3600.0 * 1.0e3
    rows.append(
        {
            "prior_id": "henley_shelton2013_high_latitude_population",
            "component": "MW",
            "label": "Henley & Shelton 2013 high-|b| population",
            "reference": "Henley & Shelton (2013)",
            "doi": "10.1088/0004-637X/773/2/92",
            "original_quantity": "intrinsic 0.5-2 keV halo surface brightness",
            "original_central": 1.5,
            "original_low": 0.5,
            "original_high": 7.0,
            "original_units": "1e-12 erg cm-2 s-1 deg-2",
            "conversion_to_primary": hs_conversion,
            "primary_central": 1.5 * hs_conversion,
            "primary_low": 0.5 * hs_conversion,
            "primary_high": 7.0 * hs_conversion,
            "scope": "population context; not a local M31-direction measurement",
            "caveat": "M31 lies outside the paper's |b|>30 degree domain; APEC conversion is conditional.",
            "provenance": "published range and median; converted with the v19 CGMsum APEC template",
        }
    )
    rows.append(
        {
            "prior_id": "locatelli2024_reference_beta0p5_m31_footprint",
            "component": "MW",
            "label": "Locatelli et al. 2024 reference beta=0.5 geometry",
            "reference": "Locatelli et al. (2024)",
            "doi": "10.1051/0004-6361/202347061",
            "original_quantity": "reference O VIII-fitted 3D density model projected to the primary footprint",
            "original_central": np.nan,
            "original_low": np.nan,
            "original_high": np.nan,
            "original_units": "Table 1 density parameters; intrinsic 0.614-0.694 keV photon emissivity",
            "conversion_to_primary": np.nan,
            "primary_central": LOCATELLI_M31_CENTRAL,
            "primary_low": LOCATELLI_M31_FIELD_LOW,
            "primary_high": LOCATELLI_M31_FIELD_HIGH,
            "parameter_error_low": LOCATELLI_M31_CENTRAL - LOCATELLI_M31_STAT_SIGMA,
            "parameter_error_high": LOCATELLI_M31_CENTRAL + LOCATELLI_M31_STAT_SIGMA,
            "north_primary": LOCATELLI_M31_NORTH,
            "south_primary": LOCATELLI_M31_SOUTH,
            "conversion_check_low": LOCATELLI_DIRECT_EM_LOW,
            "conversion_check_high": LOCATELLI_DIRECT_EM_HIGH,
            "scope": "reference Combined beta=0.5 geometry evaluated at each of the 14 M31 primary fields",
            "caveat": "The primary interval uses an FWHM~80 eV response-matched O VIII emissivity closure. The pale marginal-error envelope neglects unpublished covariance; the direct-EM APEC branch before the 5 percent closure correction is retained as a conversion check.",
            "provenance": "Table 1 reference geometry projected field by field; Locatelli O VIII-band emissivity compared after matching the published 80 eV smoothing, then converted with target APEC and HI4PI absorption",
        }
    )

    zhang_luminosity_density = 1.7252986024966184e36
    zhang_primary = zhang_luminosity_density * LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY
    rows.append(
        {
            "prior_id": "zhang2024_m31_mass_nominal_20kpc",
            "component": "M31",
            "label": "Zhang et al. 2024 M31-mass stack",
            "reference": "Zhang et al. (2024)",
            "doi": "10.1051/0004-6361/202449412",
            "original_quantity": "nominal beta-profile luminosity surface density at 20 kpc",
            "original_central": zhang_luminosity_density,
            "original_low": zhang_luminosity_density,
            "original_high": zhang_luminosity_density,
            "original_units": "erg s-1 kpc-2, intrinsic 0.5-2 keV",
            "conversion_to_primary": LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY,
            "primary_central": zhang_primary,
            "primary_low": zhang_primary,
            "primary_high": zhang_primary,
            "scope": "empirical M31-stellar-mass stack; characteristic 10-30 kpc value",
            "caveat": "A nominal line, not a formal local prior; beta-parameter covariance is unavailable.",
            "provenance": "published beta profile evaluated at 20 kpc and converted with the v19 APEC template",
        }
    )

    grayson = (
        (
            "eagle_agndt9",
            "EAGLE-AGNdT9",
            4.704476594777882e35,
            4.184141227254975e35,
            5.1535212647748536e35,
        ),
        (
            "eagle",
            "EAGLE",
            9.38177618116973e35,
            8.236144961922412e35,
            1.041199465882404e36,
        ),
        (
            "eagle_noagn",
            "EAGLE-NoAGN",
            1.9709847221830294e36,
            1.6858151590706473e36,
            2.1874197384055262e36,
        ),
        (
            "simba",
            "SIMBA",
            3.405962162605925e36,
            2.838273867820679e36,
            3.982108803986934e36,
        ),
        (
            "simba_noagn",
            "SIMBA-NoAGN",
            2.2218741538234737e37,
            2.1647474653904115e37,
            2.340689955138248e37,
        ),
    )
    for identifier, label, central, low, high in grayson:
        rows.append(
            {
                "prior_id": f"grayson2025_{identifier}",
                "component": "M31",
                "label": f"Grayson et al. 2025 {label}",
                "reference": "Grayson et al. (2025)",
                "doi": "10.3847/1538-4357/ae100f",
                "original_quantity": "digitized 10-30 kpc bin in the M31-mass stellar bin",
                "original_central": central,
                "original_low": low,
                "original_high": high,
                "original_units": "erg s-1 kpc-2, intrinsic 0.5-2 keV",
                "conversion_to_primary": LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY,
                "primary_central": central * LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY,
                "primary_low": low * LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY,
                "primary_high": high * LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY,
                "scope": "simulation prediction matched to the 11<log(M*/Msun)<11.25 stack",
                "caveat": "Digitized from Fig. 5; conversion applies one v19 APEC spectrum rather than each simulation's multiphase spectrum.",
                "provenance": "author arXiv source figure; calibrated log axes and exact line/fill RGB values",
            }
        )
    ledger = pd.DataFrame(rows)
    ledger["primary_central_definition"] = (
        "converted characteristic central; see scope and caveat"
    )
    ledger["primary_interval_definition"] = (
        "converted published or digitized interval; see scope and caveat"
    )
    ledger["all_field_inverse_variance"] = np.nan
    locatelli = ledger["prior_id"].eq("locatelli2024_reference_beta0p5_m31_footprint")
    ledger.loc[locatelli, "primary_central_definition"] = (
        "side-balanced mean of North/South inverse-variance field predictions"
    )
    ledger.loc[locatelli, "primary_interval_definition"] = (
        "14-field min-max footprint variation; not an uncertainty interval"
    )
    ledger.loc[locatelli, "all_field_inverse_variance"] = LOCATELLI_M31_ALL_FIELD
    return ledger


def _save_figure(fig: plt.Figure, stem: str, column: str) -> None:
    for suffix in ("png", "pdf"):
        path = Path(OUTPUT_DIR) / f"{stem}_{column}.{suffix}"
        fig.savefig(path, dpi=400 if suffix == "png" else None, bbox_inches="tight")


def positive_constraint_curve(
    total: float,
    error: float,
    floor: float = 0.01,
    samples: int = 800,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Return a positive-quadrant x+y=total curve and its 1-sigma band."""
    if total <= 2.0 * floor:
        raise ValueError("total must exceed twice the positive plotting floor")
    x = np.geomspace(floor, total - floor, samples)
    central = total - x
    lower = np.maximum(total - error - x, floor)
    upper = total + error - x
    return x, lower, central, upper


def _constraint_plane(priors: pd.DataFrame, column: str) -> None:
    width = 3.35 if column == "1col" else 7.0
    height = 4.6 if column == "1col" else 5.0
    fig, ax = plt.subplots(figsize=(width, height))
    x = np.linspace(0.0, 2.55, 1200)

    constraints = (
        (NORTH_FLUX, NORTH_ERROR, "#2166ac", "N/NW"),
        (SOUTH_FLUX, SOUTH_ERROR, "#b2182b", "S/SE"),
    )
    for total, error, color, label in constraints:
        lower = total - error - x
        upper = total + error - x
        ax.fill_between(
            x, np.maximum(lower, 0.0), np.maximum(upper, 0.0), color=color, alpha=0.18
        )
        valid = total - x >= 0.0
        ax.plot(
            x[valid],
            total - x[valid],
            color=color,
            lw=1.5,
            label=rf"{label}: ${total:.3f}\pm{error:.3f}$",
        )
    valid = ALL_FLUX - x >= 0.0
    ax.plot(
        x[valid],
        ALL_FLUX - x[valid],
        color="black",
        lw=1.1,
        ls="--",
        label=rf"All: ${ALL_FLUX:.3f}\pm{ALL_ERROR:.3f}$",
    )

    loc = priors.loc[
        priors["prior_id"].eq("locatelli2024_reference_beta0p5_m31_footprint")
    ].iloc[0]
    ax.axvspan(
        loc.parameter_error_low,
        loc.parameter_error_high,
        facecolor="#e69f00",
        alpha=0.035,
        edgecolor="none",
        zorder=0,
    )
    ax.axvspan(
        loc.primary_low,
        loc.primary_high,
        facecolor="#e69f00",
        alpha=0.11,
        hatch="///",
        edgecolor="#e69f00",
        lw=0.0,
        zorder=0,
    )
    ax.axvline(loc.primary_central, color="#e69f00", lw=1.1, ls="-.")

    ybar = 0.105
    ax.plot(
        [loc.parameter_error_low, loc.parameter_error_high],
        [ybar, ybar],
        color="#e69f00",
        lw=0.8,
        alpha=0.65,
    )
    ax.plot(
        [loc.primary_low, loc.primary_high],
        [ybar, ybar],
        color="#d17c00",
        lw=3.0,
        solid_capstyle="butt",
    )
    ax.plot(loc.primary_central, ybar, marker="|", color="#7a4700", ms=9)
    ax.plot(
        [loc.north_primary, loc.south_primary],
        [ybar, ybar],
        linestyle="none",
        marker="o",
        ms=2.6,
        color="#7a4700",
    )
    ax.plot(
        [loc.conversion_check_low, loc.conversion_check_high],
        [ybar - 0.020, ybar - 0.020],
        color="#b36b00",
        lw=1.1,
        ls=":",
    )
    ax.text(
        2.50,
        ybar + 0.040,
        "Locatelli+24 ref. "
        + r"$\beta=0.5$"
        + "\nM31 footprint (O VIII response matched)",
        fontsize=5.7 if column == "1col" else 7.2,
        color="#8a5000",
        ha="right",
        va="bottom",
    )
    ax.text(
        loc.conversion_check_low,
        ybar - 0.052,
        "direct-EM conversion check",
        fontsize=5.1 if column == "1col" else 6.4,
        color="#8a5000",
    )

    colors = {
        "grayson2025_eagle_agndt9": "#008080",
        "grayson2025_eagle": "#000080",
        "grayson2025_eagle_noagn": "#56b4e9",
        "grayson2025_simba": "#8b0000",
        "zhang2024_m31_mass_nominal_20kpc": "black",
    }
    label_text = {
        "grayson2025_eagle_agndt9": "EAGLE-AGNdT9",
        "grayson2025_eagle": "EAGLE",
        "grayson2025_eagle_noagn": "EAGLE-NoAGN",
        "grayson2025_simba": "SIMBA",
        "zhang2024_m31_mass_nominal_20kpc": "Zhang+24 stack",
    }
    for _, row in priors.loc[priors["component"].eq("M31")].iterrows():
        identifier = str(row.prior_id)
        if identifier == "grayson2025_simba_noagn":
            continue
        color = colors[identifier]
        low = float(row.primary_low)
        high = float(row.primary_high)
        central = float(row.primary_central)
        if high > low:
            ax.axhspan(low, high, color=color, alpha=0.075, zorder=0)
        linestyle = "--" if identifier.startswith("zhang") else "-"
        ax.axhline(central, color=color, lw=1.0, ls=linestyle)
        ax.text(
            2.52,
            central,
            f"{label_text[identifier]}  {central:.2f}",
            color=color,
            fontsize=5.8 if column == "1col" else 7.0,
            va="center",
            ha="right",
            bbox={"facecolor": "white", "edgecolor": "none", "alpha": 0.72, "pad": 0.5},
        )
        if 0.0 <= ALL_FLUX - central <= 2.55:
            ax.plot(
                ALL_FLUX - central,
                central,
                marker="o",
                ms=3.3,
                markerfacecolor="white",
                markeredgecolor=color,
                markeredgewidth=0.9,
                zorder=6,
            )
    noagn = priors.loc[priors["prior_id"].eq("grayson2025_simba_noagn")].iloc[0]
    ax.annotate(
        f"SIMBA-NoAGN: {noagn.primary_central:.1f} (off scale)",
        xy=(2.28, 2.095),
        xytext=(1.60, 2.095),
        arrowprops={"arrowstyle": "-|>", "color": "#fa8072", "lw": 0.9},
        color="#c55449",
        fontsize=6.0 if column == "1col" else 7.2,
        va="center",
    )

    ax.set_xlim(0.0, 2.55)
    ax.set_ylim(0.0, 2.2)
    ax.set_xlabel(r"MW CGM contribution, $x$")
    ax.set_ylabel(r"M31 inner-halo contribution, $y$")
    ax.text(
        0.02,
        0.98,
        "LG set to zero; non-negative components",
        transform=ax.transAxes,
        ha="left",
        va="top",
        fontsize=5.8 if column == "1col" else 7.2,
    )
    ax.text(
        0.02,
        0.945,
        "Absorbed 0.4--1.25 keV; units $10^{-15}$ erg cm$^{-2}$ s$^{-1}$ arcmin$^{-2}$",
        transform=ax.transAxes,
        ha="left",
        va="top",
        fontsize=4.8 if column == "1col" else 6.2,
    )
    ax.legend(
        loc="upper left",
        bbox_to_anchor=(0.02, 0.90),
        fontsize=5.3 if column == "1col" else 6.8,
        frameon=True,
    )
    ax.tick_params(direction="in", top=True, right=True)
    _save_figure(fig, "m31_cgmsum_mw_m31_constraint_plane", column)
    plt.close(fig)


def _constraint_plane_loglog(priors: pd.DataFrame, column: str) -> None:
    width = 3.35 if column == "1col" else 7.0
    height = 4.6 if column == "1col" else 5.0
    floor = 0.01
    fig, ax = plt.subplots(figsize=(width, height))

    constraints = (
        (NORTH_FLUX, NORTH_ERROR, "#2166ac", "N/NW"),
        (SOUTH_FLUX, SOUTH_ERROR, "#b2182b", "S/SE"),
    )
    for total, error, color, label in constraints:
        x, lower, central, upper = positive_constraint_curve(total, error, floor=floor)
        ax.fill_between(x, lower, upper, color=color, alpha=0.16)
        ax.plot(
            x,
            central,
            color=color,
            lw=1.5,
            label=rf"{label}: ${total:.3f}\pm{error:.3f}$",
        )
    x, _, central, _ = positive_constraint_curve(ALL_FLUX, ALL_ERROR, floor=floor)
    ax.plot(
        x,
        central,
        color="black",
        lw=1.1,
        ls="--",
        label=rf"All: ${ALL_FLUX:.3f}\pm{ALL_ERROR:.3f}$",
    )

    loc = priors.loc[
        priors["prior_id"].eq("locatelli2024_reference_beta0p5_m31_footprint")
    ].iloc[0]
    ax.axvspan(
        loc.parameter_error_low,
        loc.parameter_error_high,
        facecolor="#e69f00",
        alpha=0.035,
        edgecolor="none",
        zorder=0,
    )
    ax.axvspan(
        loc.primary_low,
        loc.primary_high,
        facecolor="#e69f00",
        alpha=0.11,
        hatch="///",
        edgecolor="#e69f00",
        lw=0.0,
        zorder=0,
    )
    ax.axvline(loc.primary_central, color="#e69f00", lw=1.1, ls="-.")
    ax.plot(
        loc.primary_central,
        floor * 1.18,
        marker="v",
        ms=4.2,
        color="#8a5000",
        clip_on=False,
        zorder=7,
    )
    ax.annotate(
        "Locatelli+24 MW\nno non-negative M31 intersection\nwith all-field total",
        xy=(loc.primary_central, floor * 1.18),
        xytext=(4, 10),
        textcoords="offset points",
        fontsize=5.0 if column == "1col" else 6.5,
        color="#8a5000",
        ha="left",
        va="bottom",
    )

    colors = {
        "grayson2025_eagle_agndt9": "#008080",
        "grayson2025_eagle": "#000080",
        "grayson2025_eagle_noagn": "#56b4e9",
        "grayson2025_simba": "#8b0000",
        "grayson2025_simba_noagn": "#cc79a7",
        "zhang2024_m31_mass_nominal_20kpc": "black",
    }
    label_text = {
        "grayson2025_eagle_agndt9": "EAGLE-AGNdT9",
        "grayson2025_eagle": "EAGLE",
        "grayson2025_eagle_noagn": "EAGLE-NoAGN",
        "grayson2025_simba": "SIMBA",
        "grayson2025_simba_noagn": "SIMBA-NoAGN",
        "zhang2024_m31_mass_nominal_20kpc": "Zhang+24 stack",
    }
    label_offsets = {
        "zhang2024_m31_mass_nominal_20kpc": -6,
        "grayson2025_eagle_noagn": 6,
    }
    for _, row in priors.loc[priors["component"].eq("M31")].iterrows():
        identifier = str(row.prior_id)
        color = colors[identifier]
        low = float(row.primary_low)
        high = float(row.primary_high)
        model_central = float(row.primary_central)
        if high > low:
            ax.axhspan(low, high, color=color, alpha=0.075, zorder=0)
        linestyle = "--" if identifier.startswith("zhang") else "-"
        ax.axhline(model_central, color=color, lw=1.0, ls=linestyle)
        ax.annotate(
            f"{label_text[identifier]}  {model_central:.2f}",
            xy=(2.42, model_central),
            xytext=(0, label_offsets.get(identifier, 0)),
            textcoords="offset points",
            color=color,
            fontsize=5.0 if column == "1col" else 6.5,
            va="center",
            ha="right",
            bbox={
                "facecolor": "white",
                "edgecolor": "none",
                "alpha": 0.72,
                "pad": 0.4,
            },
        )
        residual = ALL_FLUX - model_central
        if residual >= floor:
            ax.plot(
                residual,
                model_central,
                marker="o",
                ms=3.3,
                markerfacecolor="white",
                markeredgecolor=color,
                markeredgewidth=0.9,
                zorder=6,
            )

    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.set_xlim(floor, 2.5)
    ax.set_ylim(floor, 20.0)
    ax.set_xlabel(r"MW CGM contribution, $x$")
    ax.set_ylabel(r"M31 inner-halo contribution, $y$")
    ax.set_title(
        "Positive quadrant only; log-log diagnostic\n"
        "Absorbed 0.4--1.25 keV; units "
        "$10^{-15}$ erg cm$^{-2}$ s$^{-1}$ arcmin$^{-2}$",
        loc="left",
        pad=6,
        fontsize=5.4 if column == "1col" else 7.0,
    )
    ax.legend(
        loc="upper left",
        bbox_to_anchor=(0.02, 0.90),
        fontsize=5.3 if column == "1col" else 6.8,
        frameon=True,
    )
    ax.grid(which="major", color="0.88", lw=0.5)
    ax.tick_params(direction="in", which="both", top=True, right=True)
    _save_figure(fig, "m31_cgmsum_mw_m31_constraint_plane_loglog", column)
    plt.close(fig)


def _covariance_ellipse(
    mean: np.ndarray, covariance: np.ndarray, delta_chi2: float, **kwargs: object
) -> Ellipse:
    eigenvalues, eigenvectors = np.linalg.eigh(covariance)
    order = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[order]
    eigenvectors = eigenvectors[:, order]
    angle = np.degrees(np.arctan2(eigenvectors[1, 0], eigenvectors[0, 0]))
    scale = np.sqrt(delta_chi2)
    return Ellipse(
        xy=mean,
        width=2.0 * scale * np.sqrt(eigenvalues[0]),
        height=2.0 * scale * np.sqrt(eigenvalues[1]),
        angle=angle,
        **kwargs,
    )


def _residual_significance(
    slope: float, total: float, difference: float, covariance: np.ndarray
) -> float:
    residual = difference - slope * total
    variance = (
        covariance[1, 1] + slope**2 * covariance[0, 0] - 2.0 * slope * covariance[0, 1]
    )
    return residual / np.sqrt(variance)


def _predictive_plane(column: str) -> pd.DataFrame:
    width = 3.35 if column == "1col" else 7.0
    height = 3.35 if column == "1col" else 4.5
    fig, ax = plt.subplots(figsize=(width, height))
    total, difference, covariance = measurement_total_difference()
    slopes = model_contrast_slopes()
    m31_slope = float(slopes["Zhang 2024 symmetric M31"])
    mw_slope = float(slopes["Locatelli 2024 MW reference"])
    mw_low, mw_high = slopes["Locatelli 2024 MW sensitivity family"]
    mw_low = float(mw_low)
    mw_high = float(mw_high)

    xx = np.linspace(0.0, 1.65, 500)
    ax.fill_between(
        xx,
        m31_slope * xx,
        mw_high * xx,
        color="#999999",
        alpha=0.12,
        label="Non-negative MW/M31 mixtures",
    )
    ax.plot(
        xx,
        m31_slope * xx,
        color="#009e73",
        lw=1.5,
        label=rf"Zhang+24 symmetric M31 ($D/T={m31_slope:+.3f}$)",
    )
    ax.fill_between(xx, mw_low * xx, mw_high * xx, color="#e69f00", alpha=0.18)
    ax.plot(
        xx,
        mw_slope * xx,
        color="#d17c00",
        lw=1.5,
        label=rf"Locatelli+24 MW ($D/T={mw_slope:+.3f}$)",
    )
    ax.axhline(0.0, color="0.35", lw=1.0, ls="--", label="Uniform component")

    mean = np.array([total, difference])
    ax.add_patch(
        _covariance_ellipse(
            mean, covariance, 5.991, facecolor="#2166ac", edgecolor="none", alpha=0.10
        )
    )
    ax.add_patch(
        _covariance_ellipse(
            mean,
            covariance,
            2.279,
            facecolor="#2166ac",
            edgecolor="#2166ac",
            lw=1.0,
            alpha=0.22,
        )
    )
    ax.errorbar(
        total,
        difference,
        fmt="o",
        color="#08306b",
        ms=4.2,
        zorder=5,
        label="Observed North-South estimator",
    )

    model_rows = []
    for name, slope, color in (
        ("Uniform", 0.0, "0.25"),
        ("Zhang 2024 symmetric M31", m31_slope, "#009e73"),
        ("Locatelli 2024 MW reference", mw_slope, "#d17c00"),
    ):
        predicted = slope * total
        significance = _residual_significance(slope, total, difference, covariance)
        ax.plot(total, predicted, marker="s", ms=4.0, color=color, zorder=4)
        model_rows.append(
            {
                "model": name,
                "contrast_slope_D_over_T": slope,
                "observed_total": total,
                "predicted_difference_at_observed_total": predicted,
                "observed_difference": difference,
                "residual_sigma_with_TD_covariance": significance,
            }
        )

    ax.text(
        0.03,
        0.97,
        "Absorbed 0.4--1.25 keV; units $10^{-15}$ erg cm$^{-2}$ s$^{-1}$ arcmin$^{-2}$\nObserved: "
        + rf"$T={total:.3f}$, $D={difference:.3f}$; "
        + rf"M31-only offset={model_rows[1]['residual_sigma_with_TD_covariance']:.2f}$\sigma$",
        transform=ax.transAxes,
        ha="left",
        va="top",
        fontsize=5.5 if column == "1col" else 7.2,
    )
    ax.text(
        0.98,
        0.96,
        "68/95% covariance",
        transform=ax.transAxes,
        ha="right",
        va="top",
        color="#2166ac",
        fontsize=5.1 if column == "1col" else 6.8,
    )
    ax.set_xlim(0.0, 1.65)
    ax.set_ylim(-0.48, 0.70)
    ax.set_xlabel(r"Side-symmetrized total, $T=(F_N+F_S)/2$")
    ax.set_ylabel(r"North--South difference, $D=F_N-F_S$")
    handles, _ = ax.get_legend_handles_labels()
    labels = [
        "Non-negative mixtures",
        rf"Symmetric M31 ($D/T={m31_slope:+.3f}$)",
        rf"MW model ($D/T={mw_slope:+.3f}$)",
        "Uniform component",
    ]
    ax.legend(
        handles[:4],
        labels,
        loc="lower right",
        fontsize=5.2 if column == "1col" else 6.8,
        frameon=True,
    )
    ax.tick_params(direction="in", top=True, right=True)
    _save_figure(fig, "m31_cgmsum_total_difference_predictive_plane", column)
    plt.close(fig)
    return pd.DataFrame(model_rows)


def main() -> None:
    plt.style.use(["science", "no-latex"])
    priors = build_prior_ledger()
    priors.to_csv(PRIOR_LEDGER_CSV, index=False)
    for column in ("1col", "2col"):
        _constraint_plane(priors, column)
        _constraint_plane_loglog(priors, column)
    predictions = _predictive_plane("1col")
    _predictive_plane("2col")
    predictions.to_csv(PREDICTIONS_CSV, index=False)
    print(f"Wrote {PRIOR_LEDGER_CSV}")
    print(f"Wrote {PREDICTIONS_CSV}")
    for path in expected_output_paths():
        print(path)


if __name__ == "__main__":
    main()