proteus.interior_struct.zalmoxis¶

def get_zalmoxis_output_filepath(outdir: str):
    """Returns the output file path for Zalmoxis data.
    Args:
        outdir (str): Output directory.
    Returns:
        str: Path to the output file.
    """
    return os.path.join(outdir, 'data', 'zalmoxis_output.dat')

def validate_zalmoxis_output_schema(
    output_path: str,
    hf_row: dict,
    rtol_radius: float = 1e-6,
    rtol_mass: float = 5e-2,
    mantle_mass_ref: float | None = None,
) -> None:
    """Verify zalmoxis_output.dat is consistent with hf_row scalars.

    The file is the contract Aragog reads inside ``solver.reset()``
    (eos_method=2). This check confirms the file's last r matches
    ``hf_row['R_int']``, and that the mantle mass matches
    ``hf_row['M_int'] - hf_row['M_core']``.
    Catches file I/O corruption, column-order mistakes, truncation,
    encoding drift, and Aragog/Zalmoxis schema desync at the file
    handover boundary.

    The mantle mass is taken from ``mantle_mass_ref`` when the caller supplies
    it: the structure's accumulator total ``mass_enclosed[-1]`` minus the exact
    core-mass target ``cmb_mass``. Re-integrating the coarse output nodes with a
    grid trapezoid instead diverges from the sub-grid-substepped RK45 integral
    across the steep interior density profile, reaching ~10% at high planet mass,
    which would false-reject a structure that actually conserves mass. With the
    reference supplied, part (b) reduces to a check that the core-mantle split
    point is well resolved: the residual it measures is the CMB-node snap
    overshoot ``mass_enclosed[cmb_index] - cmb_mass`` relative to the mantle
    mass, a fraction of one boundary shell. File-density-column corruption is
    caught downstream by Aragog's EOS-vs-mesh consistency at ``solver.reset()``;
    the radius check (part (a)) catches truncation and column swaps directly.
    When no reference is supplied the check falls back to the grid-trapezoidal
    shell-sum.

    Parameters
    ----------
    output_path : str
        Path to the just-written ``zalmoxis_output.dat``.
    hf_row : dict
        PROTEUS hf_row holding the scalar truth (R_int, M_int, M_core).
    rtol_radius : float
        Relative tolerance for the top-of-mantle vs. R_int check.
        Default 1e-6: the file's last r and ``hf_row['R_int']`` come
        from the same variable in zalmoxis_solver, so equality is
        exact modulo float-string round-trip noise. Tight tolerance
        catches truncation, last-line corruption, and column-swap
        bugs at the bit level.
    rtol_mass : float
        Relative tolerance for the integrated mantle mass vs.
        ``M_int - M_core``. Default 5e-2 (5%) reflects two stacked
        sources of legitimate mismatch:
        (a) integrator-method difference: Zalmoxis' ``mass_enclosed``
        is the ODE state from RK45 with sub-grid substepping, while
        the schema check re-integrates via a grid-trapezoidal
        shell-sum. On stiff CHILI density profiles this shows
        about 0.8 to 2.0 %.
        (b) ``blend_mesh_files`` post-write modifies the file in
        place but does not update hf_row scalars. When blending fires
        with alpha < 1 (capping large R-shifts), the file's integrated
        mass can drift up to about 5% from the unblended
        hf_row['M_int'], so this check has to tolerate it.
        5e-2 keeps a >2x margin over the worst legitimate noise
        while still catching gross corruption (column swap,
        truncation, byte-flip) at >>5 %. **The tight mass-conservation
        contract (<0.1 %) lives in the wrapper-level mass-anchor
        check on hf_row['M_int'] / hf_row['M_int_target']**, not here.

    Raises
    ------
    RuntimeError
        On any violation (file unreadable, wrong shape, radius
        mismatch, mass mismatch). Caller (the wrapper's except-block)
        catches this and routes through the fall-back path.
    """
    try:
        data = np.loadtxt(output_path)
    except Exception as exc:
        raise RuntimeError(
            'zalmoxis_output.dat schema violation: could not reload '
            f'the just-written file ({output_path}): {exc}'
        )
    if data.size == 0 or data.ndim != 2 or data.shape[1] != 5:
        raise RuntimeError(
            'zalmoxis_output.dat schema violation: unexpected shape '
            f'{data.shape if data.size else "empty"} '
            f'(expected (N, 5)) for {output_path}'
        )

    r_file = data[:, 0]
    rho_file = data[:, 2]

    # (a) top-of-mantle radius == hf_row['R_int']
    r_top = float(r_file[-1])
    r_int_hf = float(hf_row.get('R_int', 0.0))
    if r_int_hf > 0:
        r_rel = abs(r_top / r_int_hf - 1.0)
        if r_rel > rtol_radius:
            raise RuntimeError(
                'zalmoxis_output.dat schema violation: top-of-mantle '
                f'r={r_top:.6e} from file differs from '
                f'hf_row[R_int]={r_int_hf:.6e} '
                f'(rel={r_rel:.3e} > {rtol_radius:.1e})'
            )

    # (b) mantle integrated mass == hf_row['M_int'] - hf_row['M_core'].
    # Prefer the structure's own RK45 accumulator mantle mass when the
    # caller supplies it (mantle_mass_ref). That value is the sub-grid
    # substepped ODE integral, exact to the solver tolerance. Re-integrating
    # the coarse output nodes with a grid trapezoid instead diverges from it
    # across the steep core-mantle density jump (the CMB-node snap attributes
    # a whole boundary shell to one side), reaching ~10% at high planet mass,
    # which would false-reject a structure that actually conserves mass. Only
    # when no reference is available do we fall back to the grid-trapezoidal
    # shell-sum (4/3 pi (r2^3 - r1^3) * rho_avg).
    if mantle_mass_ref is not None:
        mantle_mass = float(mantle_mass_ref)
        mass_source = 'accumulator'
    else:
        shells = (
            (4.0 / 3.0)
            * np.pi
            * (r_file[1:] ** 3 - r_file[:-1] ** 3)
            * 0.5
            * (rho_file[1:] + rho_file[:-1])
        )
        mantle_mass = float(np.sum(shells))
        mass_source = 'trapezoid'
    M_int_hf = float(hf_row.get('M_int', 0.0))
    M_core_hf = float(hf_row.get('M_core', 0.0))
    expected_mantle = M_int_hf - M_core_hf
    if expected_mantle > 0:
        m_rel = abs(mantle_mass / expected_mantle - 1.0)
        log.debug(
            'zalmoxis mass-closure: mantle split m_rel=%.3e (%s, rtol_mass=%.1e)',
            m_rel,
            mass_source,
            rtol_mass,
        )
        if m_rel > rtol_mass:
            raise RuntimeError(
                'zalmoxis_output.dat schema violation: '
                f'mantle mass ({mass_source})={mantle_mass:.6e} kg '
                'differs from hf_row[M_int - M_core]='
                f'{expected_mantle:.6e} kg '
                f'(rel={m_rel:.3e} > {rtol_mass:.1e})'
            )

def build_volatile_profile(hf_row: dict, mantle_eos: str):
    """Build a VolatileProfile from helpfile volatile masses.

    Computes per-phase (liquid/solid) mass fractions for dissolved volatiles
    that have Zalmoxis EOS tables. Returns None if no volatiles are dissolved
    or if the mantle liquid/solid masses are unavailable.

    Ownership of the liquid/solid split: the coupled PROTEUS path takes
    the per-phase masses from the outgassing chemistry (CALLIOPE/Aragog
    equilibrium via the helpfile), which already solved solubility.
    Zalmoxis's own ``partition_rule`` hook is a structure-side
    idealization for standalone Zalmoxis runs with no chemistry coupled;
    it is deliberately not used here.

    Parameters
    ----------
    hf_row : dict
        Current helpfile row with volatile mass keys (e.g. ``H2O_kg_liquid``,
        ``H2O_kg_solid``) and mantle mass keys: ``M_mantle`` with
        ``Phi_global`` (preferred, structural), or ``M_mantle_liquid`` and
        ``M_mantle_solid`` as a fallback.
    mantle_eos : str
        Primary mantle EOS identifier (e.g. ``'PALEOS:MgSiO3'``).

    Returns
    -------
    VolatileProfile or None
        Profile with per-phase fractions, or None if not applicable.
    """
    from zalmoxis.mixing import VolatileProfile

    # Denominators: the melt concentration must be taken against the same
    # reservoir the chemistry used and the structure carries. CALLIOPE
    # dissolves into M_mantle * Phi_global, and the structure integrates
    # the (structural) mantle mass, so split the structural M_mantle by
    # Phi_global. Aragog's M_mantle_liquid/M_mantle_solid use its PALEOS
    # EOS cell masses instead, ~7% denser than the wet structural mantle,
    # which diluted w_liquid and lost the same fraction of dissolved
    # water from the structure. Fall back to the per-phase masses when
    # the structural fields are absent (e.g. minimal test rows).
    M_mantle = float(hf_row.get('M_mantle', 0.0))
    Phi_global = float(hf_row.get('Phi_global', float('nan')))
    if M_mantle > 0 and np.isfinite(Phi_global):
        Phi_global = min(max(Phi_global, 0.0), 1.0)
        M_liq = Phi_global * M_mantle
        M_sol = (1.0 - Phi_global) * M_mantle
    else:
        M_liq = float(hf_row.get('M_mantle_liquid', 0.0))
        M_sol = float(hf_row.get('M_mantle_solid', 0.0))

    # Need mantle mass data to compute fractions
    if M_liq + M_sol <= 0:
        return None

    w_liquid = {}
    w_solid = {}
    has_nonzero = False

    for species, eos_name in _VOLATILE_EOS_MAP.items():
        kg_liq = float(hf_row.get(f'{species}_kg_liquid', 0.0))
        kg_sol = float(hf_row.get(f'{species}_kg_solid', 0.0))

        # Only include species with meaningful dissolved mass
        if kg_liq + kg_sol <= 0:
            continue

        # Mass fraction in liquid phase
        w_l = kg_liq / M_liq if M_liq > 0 else 0.0
        # Mass fraction in solid phase
        w_s = kg_sol / M_sol if M_sol > 0 else 0.0

        w_liquid[eos_name] = w_l
        w_solid[eos_name] = w_s
        has_nonzero = True

    if not has_nonzero:
        return None

    # Normalize: total volatile fraction in each phase must not exceed 0.95
    # (at least 5% silicate). Clamp proportionally if sum exceeds the limit.
    for w_dict in (w_liquid, w_solid):
        total = sum(w_dict.values())
        max_volatile_frac = 0.95
        if total > max_volatile_frac:
            scale = max_volatile_frac / total
            for k in w_dict:
                w_dict[k] *= scale

    log.info(
        'Built VolatileProfile: liquid=%s, solid=%s',
        {k: f'{v:.4f}' for k, v in w_liquid.items()},
        {k: f'{v:.4f}' for k, v in w_solid.items()},
    )

    return VolatileProfile(
        w_liquid=w_liquid,
        w_solid=w_solid,
        primary_component=mantle_eos,
    )

def extend_mantle_eos_with_volatiles(mantle_eos: str, volatile_profile) -> str:
    """Extend a single-component mantle EOS string with volatile components.

    If a VolatileProfile is provided and the mantle EOS is a single component,
    this adds the volatile EOS components with small placeholder fractions.
    The actual fractions are overridden at each ODE step by the VolatileProfile.

    Parameters
    ----------
    mantle_eos : str
        Base mantle EOS string (e.g. ``'PALEOS:MgSiO3'``).
    volatile_profile : VolatileProfile or None
        Profile containing volatile EOS component names.

    Returns
    -------
    str
        Extended EOS string (e.g.
        ``'PALEOS:MgSiO3:0.98+PALEOS:H2O:0.01+Chabrier:H:0.01'``),
        or the original string if no extension needed.
    """
    if volatile_profile is None:
        return mantle_eos

    # Don't modify if already multi-component
    if '+' in mantle_eos:
        return mantle_eos

    # Collect all volatile EOS components from the profile
    all_vol_components = set()
    for d in (volatile_profile.w_liquid, volatile_profile.w_solid):
        all_vol_components.update(d.keys())

    if not all_vol_components:
        return mantle_eos

    # Build extended string with small placeholder fractions
    # (actual fractions set by VolatileProfile at each radius)
    n_vol = len(all_vol_components)
    placeholder = 0.01  # 1% each
    primary_frac = max(0.5, 1.0 - n_vol * placeholder)
    parts = [f'{mantle_eos}:{primary_frac:.4f}']
    for comp in sorted(all_vol_components):
        parts.append(f'{comp}:{placeholder:.4f}')

    extended = '+'.join(parts)
    log.info('Extended mantle EOS: %s -> %s', mantle_eos, extended)
    return extended

def solve_superliquidus_adiabat(config: Config, hf_row: dict | None) -> dict:
    """Solve for the coolest fully molten adiabat with a controlled superheat.

    The ``liquidus_super`` initial condition starts the mantle on a single
    adiabat (uniform specific entropy) that lies a minimum of
    ``config.planet.delta_T_super`` Kelvin above the configured liquidus,
    evaluated at the most-constraining mantle depth. Solving for the surface
    temperature that achieves this margin fixes the whole isentropic profile.

    The superheat is checked against whatever solidus/liquidus parameterisation
    is configured, and the surface-temperature search window is anchored to the
    surface liquidus, so the initial condition adapts to the melting curve in
    use rather than relying on a fixed surface temperature or entropy value.

    When the binding (minimum-superheat) depth is shallow, as it is for the
    PALEOS MgSiO3 liquidus, the solved entropy is essentially independent of
    planet mass, which keeps a mass grid on a common initial adiabat. That is a
    property of the adiabat-vs-liquidus slope ordering, not a guarantee: for a
    steeper liquidus the binding can migrate toward the core-mantle boundary, in
    which case the solve logs a warning (the surface anchor is then weakly
    constrained and approaches a CMB-liquidus anchor in the melting curve's
    extrapolated regime).

    Parameters
    ----------
    config : Config
        PROTEUS configuration. Uses ``planet.delta_T_super`` (required superheat
        in K), ``planet.mass_tot``, ``interior_struct.core_frac`` and
        ``interior_struct.zalmoxis.mantle_eos``.
    hf_row : dict or None
        Helpfile row. ``hf_row['P_cmb']`` is used when populated; otherwise a
        Noack & Lasbleis (2020) mass-aware estimate is used on the first call.

    Returns
    -------
    dict
        ``surface_T`` [K], ``S_target`` [J/(kg K)], ``cmb_T`` [K],
        ``achieved_superheat`` [K], ``binding_P`` [Pa] and ``P_cmb`` [Pa].

    Raises
    ------
    RuntimeError
        If the requested superheat cannot be reached before the deep adiabat
        exhausts the EOS table (the mantle is too deep to be molten with that
        much superheat at this mass). The message reports the largest
        achievable superheat so the user can lower ``delta_T_super``.
    """
    try:
        from zalmoxis.eos_export import compute_entropy_adiabat
        from zalmoxis.melting_curves import paleos_liquidus
    except (ImportError, ModuleNotFoundError) as e:
        raise RuntimeError(
            'liquidus_super mode requires Zalmoxis '
            '(zalmoxis.eos_export.compute_entropy_adiabat and '
            'zalmoxis.melting_curves.paleos_liquidus); import failed: '
            f'{e}'
        )

    delta = float(config.planet.delta_T_super)

    P_cmb = hf_row.get('P_cmb') if isinstance(hf_row, dict) else None
    if not P_cmb or P_cmb <= 0:
        from proteus.utils.structure_estimate import estimate_P_cmb_NL20

        P_cmb = estimate_P_cmb_NL20(
            float(config.planet.mass_tot),
            float(config.interior_struct.core_frac),
            str(config.interior_struct.core_frac_mode),
        )
        log.warning(
            'liquidus_super: hf_row["P_cmb"] not yet populated; using '
            'Noack & Lasbleis (2020) mass-aware fallback P_cmb=%.1f GPa '
            '(mass_tot=%.2f M_Earth). The energetics initial condition is '
            're-derived against the converged Zalmoxis P_cmb on the next '
            'iteration.',
            P_cmb / 1e9,
            float(config.planet.mass_tot),
        )
    P_cmb = float(P_cmb)

    mantle_eos = config.interior_struct.zalmoxis.mantle_eos
    P_surface = 1e5  # 1 bar surface anchor for the adiabat
    global _SUPERLIQ_LAST_ANCHOR
    _cache_key = (round(P_cmb / 1e6), round(delta, 3), str(mantle_eos))
    if _cache_key in _SUPERLIQ_CACHE:
        cached = dict(_SUPERLIQ_CACHE[_cache_key])
        _SUPERLIQ_LAST_ANCHOR = float(cached['cmb_T'])
        return cached

    mat_dicts = load_zalmoxis_material_dictionaries()
    solid_eos, liquid_eos = resolve_2phase_mgsio3_paths(mantle_eos, mat_dicts)
    eos_file = mat_dicts.get(mantle_eos, {}).get('eos_file', '') or solid_eos or ''
    melt_funcs = load_zalmoxis_solidus_liquidus_functions(mantle_eos, config)
    if melt_funcs is not None:
        sol_func, liq_func = melt_funcs
    else:
        liq_func = paleos_liquidus

        def sol_func(P):
            return _SUPERLIQ_DEFAULT_MUSHY * np.asarray(paleos_liquidus(P))

    def _probe(T_surf: float) -> dict:
        """Adiabat from this surface T: minimum superheat over depth + validity.

        The adiabat is isentropic by construction, so ``valid`` requires a
        uniform ``S_profile``: a bracket failure (the deep adiabat exhausting
        the EOS table) plateaus the profile and breaks isentropy, which a
        non-uniform ``S_profile``, a NaN, or a cooling-with-depth segment all
        flag. This catches the plateau that a monotonicity-only test misses.
        """
        result = compute_entropy_adiabat(
            eos_file=eos_file,
            T_surface=float(T_surf),
            P_surface=P_surface,
            P_cmb=P_cmb,
            n_points=_SUPERLIQ_N_POINTS,
            solidus_func=sol_func,
            liquidus_func=liq_func,
            solid_eos_file=solid_eos,
            liquid_eos_file=liquid_eos,
        )
        P = np.asarray(result['P'], dtype=float)
        T = np.asarray(result['T'], dtype=float)
        S_target = result['S_target']
        S_prof = np.asarray(result['S_profile'], dtype=float)
        order = np.argsort(P)
        P, T = P[order], T[order]
        finite = bool(
            np.isfinite(T).all() and np.isfinite(S_target) and np.isfinite(S_prof).all()
        )
        s_drift = (
            (float(np.max(S_prof)) - float(np.min(S_prof))) / max(abs(float(S_target)), 1.0)
            if finite
            else np.inf
        )
        valid = bool(
            finite
            and s_drift < _SUPERLIQ_MAX_S_DRIFT
            and np.all(np.diff(T) > -1.0)  # no gross cooling-with-depth
        )
        liq = np.asarray(liq_func(P), dtype=float)
        i = int(np.argmin(T - liq))
        return {
            'superheat': float(T[i] - liq[i]),
            'binding_P': float(P[i]),
            'S_target': float(S_target),
            'cmb_T': float(T[-1]),
            'valid': valid,
        }

    # Scan surface temperature over a window anchored to the SURFACE LIQUIDUS,
    # so the search adapts to the configured melting curve instead of a fixed
    # Kelvin band. The coolest valid molten adiabat brackets the delta crossing
    # from below, so delta_T_super is honoured as a true minimum (no hot floor).
    T_liq_surf = float(np.asarray(liq_func(P_surface)).reshape(-1)[0])
    scan = []
    for T_surf in np.linspace(
        T_liq_surf, T_liq_surf + _SUPERLIQ_SCAN_SPAN_K, _SUPERLIQ_SCAN_STEPS
    ):
        d = _probe(float(T_surf))
        if d['valid'] and np.isfinite(d['superheat']):
            scan.append((float(T_surf), d))
            # Stop once the increasing branch has bracketed the delta crossing;
            # scanning the full span is only needed to locate the table ceiling
            # in the unreachable-superheat (error) case.
            if d['superheat'] >= delta:
                break
    if not scan:
        raise RuntimeError(
            'liquidus_super: no valid molten adiabat found while solving the '
            f'initial condition (P_cmb={P_cmb / 1e9:.0f} GPa, surface '
            f'liquidus={T_liq_surf:.0f} K). The EOS table may not support a '
            'molten mantle at this pressure.'
        )

    # Superheat rises with surface temperature until the deep adiabat hits the
    # EOS-table ceiling; restrict the solve to that increasing branch.
    superheats = [d['superheat'] for _, d in scan]
    branch = scan[: int(np.argmax(superheats)) + 1]
    if branch[-1][1]['superheat'] < delta:
        ceil_T, ceil_d = branch[-1]
        raise RuntimeError(
            'liquidus_super: cannot initialise a fully molten mantle with '
            f'delta_T_super={delta:.0f} K at P_cmb={P_cmb / 1e9:.0f} GPa. The '
            f'largest achievable superheat is {ceil_d["superheat"]:.0f} K (at '
            f'surface T={ceil_T:.0f} K) before the deep adiabat exhausts the '
            'EOS table. Lower delta_T_super or the planet mass, or extend the '
            'EOS table to higher temperature.'
        )

    # Bracket the delta crossing on the increasing branch, then bisection-refine.
    if branch[0][1]['superheat'] >= delta:
        # The coolest in-table adiabat already meets the margin (e.g. delta=0,
        # or the cool end of the valid band is itself table-limited): it is the
        # coolest fully molten adiabat available, so return it.
        T_solved = branch[0][0]
    else:
        T_lo, T_hi = branch[0][0], branch[-1][0]
        for k in range(1, len(branch)):
            if branch[k][1]['superheat'] >= delta:
                T_lo, T_hi = branch[k - 1][0], branch[k][0]
                break
        for _ in range(_SUPERLIQ_N_BISECT):
            mid = 0.5 * (T_lo + T_hi)
            dm = _probe(mid)
            if dm['valid'] and dm['superheat'] >= delta:
                T_hi = mid
            else:
                T_lo = mid
        T_solved = T_hi

    final = _probe(T_solved)
    if not final['valid']:
        raise RuntimeError(
            f'liquidus_super: solved surface T={T_solved:.0f} K yielded an '
            f'out-of-table adiabat (P_cmb={P_cmb / 1e9:.0f} GPa); the EOS table '
            'is exhausted at this mass.'
        )

    binding_frac = final['binding_P'] / P_cmb
    if final['binding_P'] > FEI2021_LIQUIDUS_P_CALIB_PA:
        # The tightest superheat is set where the liquidus is an extrapolation
        # (beyond its calibration). The margin there is uncertain and the
        # construction approaches a CMB-liquidus anchor, which gives a
        # cold-surface IC; warn so a melting-curve swap that moves the binding
        # into the extrapolated regime does not silently land there. (A binding
        # that is fractionally deep but below the calibration, as for a
        # low-mass planet whose whole mantle is shallow, is fine.)
        log.warning(
            'liquidus_super: the minimum-superheat depth (binding P=%.0f GPa) '
            'is beyond the liquidus calibration (~%.0f GPa), so the superheat '
            'margin there is set against an extrapolated liquidus. Verify the '
            'liquidus parameterisation is appropriate for this EOS and mass.',
            final['binding_P'] / 1e9,
            FEI2021_LIQUIDUS_P_CALIB_PA / 1e9,
        )

    log.info(
        'liquidus_super: surface T=%.0f K gives a fully molten adiabat at least '
        '%.0f K above the liquidus (achieved %.0f K at P=%.0f GPa = %.0f%% of '
        'P_cmb); T_cmb=%.0f K, S=%.1f J/(kg K), P_cmb=%.0f GPa.',
        T_solved,
        delta,
        final['superheat'],
        final['binding_P'] / 1e9,
        binding_frac * 100.0,
        final['cmb_T'],
        final['S_target'],
        P_cmb / 1e9,
    )
    out = {
        'surface_T': float(T_solved),
        'S_target': float(final['S_target']),
        'cmb_T': float(final['cmb_T']),
        'achieved_superheat': float(final['superheat']),
        'binding_P': float(final['binding_P']),
        'P_cmb': P_cmb,
    }
    _SUPERLIQ_CACHE[_cache_key] = dict(out)
    _SUPERLIQ_LAST_ANCHOR = float(out['cmb_T'])
    return out

def load_zalmoxis_configuration(
    config: Config,
    hf_row: dict,
    temperature_mode_override: str | None = None,
    external_temperature_source: bool = False,
):
    """Loads the model configuration for Zalmoxis and calculates the dry mass of the planet based on the total mass and the mass of volatiles.
    Args:
        config (Config): The configuration object containing the Zalmoxis parameters.
        hf_row (dict): A dictionary containing the mass of volatiles and other parameters.
        temperature_mode_override: Optional local override for
            ``config.planet.temperature_mode``. Lets callers force a different
            structure-solve mode (e.g. 'adiabatic' for SPIDER coupling with
            T-dependent mantle EOS) without mutating the shared Config object.
            When None, falls back to ``config.planet.temperature_mode``.
        external_temperature_source: True when the caller passes an external
            temperature function or array so the structure solve follows an
            evolved T(r) instead of the internal temperature-mode dispatch. In
            'liquidus_super' mode this lets ``cmb_temperature`` reuse the last
            solved super-liquidus anchor rather than re-solving it, since the
            anchor is not consumed when an external source drives the solve.
    Returns:
        dict: A dictionary containing the Zalmoxis configuration parameters.
    """

    # The Zalmoxis solver consumes core_frac as a mass fraction and does
    # not read core_frac_mode. Warn if a radius fraction was requested, so
    # the user is not surprised by the mass-fraction interpretation.
    if config.interior_struct.core_frac_mode == 'radius':
        log.warning(
            'interior_struct.core_frac_mode = "radius" has no effect with the '
            'zalmoxis module: core_frac (%.3f) is interpreted as a mass '
            'fraction. Use core_frac_mode = "mass", or switch to the dummy or '
            'spider module if a radius fraction is intended.',
            config.interior_struct.core_frac,
        )

    # Setup target planet mass (input parameter) as the total mass of the planet (dry mass + volatiles) [kg]
    total_planet_mass = config.planet.mass_tot * M_earth

    log.debug(
        'Total target planet mass (dry mass + volatiles): %s kg '
        'with EOS: core=%s, mantle=%s, ice=%s',
        total_planet_mass,
        config.interior_struct.zalmoxis.core_eos,
        config.interior_struct.zalmoxis.mantle_eos,
        config.interior_struct.zalmoxis.ice_layer_eos or 'none',
    )

    # Calculate the volatile mass excluded from the structure target.
    # Whole-planet oxygen accounting (issue #677): atmospheric+dissolved O
    # is summed alongside H/C/N/S so the dry-mass target passed to
    # Zalmoxis correctly reserves space for the O that CALLIOPE places in
    # atmospheric H2O, CO2, SO2, etc. Mantle FeO-bound O remains in M_int
    # implicitly via the PALEOS density tables; we do not double-count it.
    # With dry_mantle the full inventory is excluded (the mantle EOS
    # represents bare silicate). When the mantle EOS carries dissolved
    # volatiles (dry_mantle = false), only the atmospheric inventory may
    # be excluded: the dissolved mass is already part of the wet-mantle
    # EOS, and subtracting it again would remove it twice. Note this
    # keeps ALL dissolved species inside the target, while the EOS blend
    # only represents those in _VOLATILE_EOS_MAP (H2O in practice); the
    # dissolved mass of unmapped species (CO2, CH4, N2, S2) lands in the
    # silicate remainder fraction and is packed at silicate density.
    # Mass stays exact, compressibility is misrepresented; the radius
    # error scales as (dissolved mass fraction) x (density contrast) and
    # is negligible at typical inventories. Escaped mass is already
    # debited from the *_kg_* inventories and needs no separate term.
    # Defensive .get(): some pre-IC paths invoke Zalmoxis before
    # calc_target_elemental_inventories has populated all element columns.
    dry_mantle = config.interior_struct.zalmoxis.dry_mantle
    M_volatiles = 0.0
    for e in element_list:
        if dry_mantle:
            M_volatiles += float(hf_row.get(e + '_kg_total', 0.0))
        else:
            M_volatiles += float(hf_row.get(e + '_kg_atm', 0.0))

    log.debug(f'Volatile mass: {M_volatiles} kg')
    log.debug(
        'Mass budget: total=%.6e kg (%.4f M_earth), volatiles=%.6e kg (%.2f%%)',
        total_planet_mass,
        config.planet.mass_tot,
        M_volatiles,
        100.0 * M_volatiles / total_planet_mass if total_planet_mass > 0 else 0,
    )

    # Calculate the target planet mass (dry mass) by subtracting the mass of volatiles from the total planet mass
    planet_mass = total_planet_mass - M_volatiles

    log.debug(f'Target planet mass (dry mass): {planet_mass} kg ')

    # Build per-layer EOS config dict from PROTEUS config fields
    layer_eos_config = {
        'core': config.interior_struct.zalmoxis.core_eos,
        'mantle': config.interior_struct.zalmoxis.mantle_eos,
    }
    if config.interior_struct.zalmoxis.ice_layer_eos is not None:
        layer_eos_config['ice_layer'] = config.interior_struct.zalmoxis.ice_layer_eos

    # Mushy zone factor: controls width of partially molten region in PALEOS
    # unified EOS. Applied as T_solidus = T_liquidus * mushy_zone_factor.
    mzf = config.interior_struct.zalmoxis.mushy_zone_factor
    mushy_zone_factors = {
        'PALEOS:iron': mzf,
        'PALEOS:MgSiO3': mzf,
        'PALEOS:H2O': mzf,
    }

    zc = config.interior_struct.zalmoxis
    log.debug(
        'Zalmoxis config: num_layers=%d, outer_solver=%s, tol_outer=%.1e, '
        'tol_inner=%.1e, use_jax=%s, use_anderson=%s',
        zc.num_levels,
        zc.outer_solver,
        zc.solver_tol_outer,
        zc.solver_tol_inner,
        zc.use_jax,
        zc.use_anderson,
    )

    # Core fraction. The Zalmoxis solver consumes `core_mass_fraction`
    # literally as a mass fraction and does not read `core_frac_mode`, so
    # `core_frac` is always interpreted as a mass fraction here regardless
    # of the mode. `core_frac_mode = "radius"` is only honoured by the
    # dummy and spider structure modules; a warning is emitted above when
    # it is combined with zalmoxis.
    return {
        'planet_mass': planet_mass,
        'core_mass_fraction': config.interior_struct.core_frac,
        'core_frac_mode': config.interior_struct.core_frac_mode,
        'mantle_mass_fraction': config.interior_struct.zalmoxis.mantle_mass_fraction,
        # For the structure solve, 'accretion' and 'isentropic' both reduce
        # to 'adiabatic' inside Zalmoxis. Zalmoxis only solves the structure
        # (M-R via hydrostatic + EOS); the entropy IC for Aragog/SPIDER is
        # set independently from config.planet.ini_entropy. 'accretion'
        # delays White+Li T-profile until after structure converges;
        # 'isentropic' (CHILI protocol) means the energetics solver
        # consumes ini_entropy, not the Zalmoxis T-profile, so the choice
        # of structure-solve T is decoupled from the IC. 'adiabatic_from_cmb'
        # is passed through to Zalmoxis with the CMB-anchor temperature so
        # the structure-side T(r) integrates upward from T_cmb at R_cmb,
        # matching the entropy that the energetics solver receives via
        # compute_initial_entropy. 'liquidus_super' maps to
        # 'adiabatic_from_cmb' here, with cmb_temperature derived from the
        # Fei+2021 liquidus at the converged P_cmb (or a Noack & Lasbleis
        # (2020) mass-aware P_cmb estimate on the very first call before
        # Zalmoxis has populated P_cmb) plus delta_T_super. The energetics
        # IC step recomputes this
        # exact same anchor against the converged P_cmb, so the structure
        # solve and the entropy IC stay in agreement after the first
        # round-trip. temperature_mode_override lets SPIDER coupling force
        # adiabatic without mutating the shared Config object (see proteus
        # rules §"Config mutability").
        'temperature_mode': _resolve_zalmoxis_temperature_mode(
            temperature_mode_override or config.planet.temperature_mode
        ),
        'surface_temperature': config.planet.tsurf_init,
        'cmb_temperature': _resolve_zalmoxis_cmb_temperature(
            config,
            hf_row,
            temperature_mode_override or config.planet.temperature_mode,
            external_temperature_source=external_temperature_source,
        ),
        'center_temperature': config.planet.tcenter_init,
        'temp_profile_file': None,
        'layer_eos_config': layer_eos_config,
        'mushy_zone_factor': mzf,
        'mushy_zone_factors': mushy_zone_factors,
        'num_layers': config.interior_struct.zalmoxis.num_levels,
        'target_surface_pressure': _get_target_surface_pressure(config, hf_row),
        # Solver tolerances and iteration limits
        'tolerance_outer': config.interior_struct.zalmoxis.solver_tol_outer,
        'tolerance_inner': config.interior_struct.zalmoxis.solver_tol_inner,
        'max_iterations_outer': config.interior_struct.zalmoxis.solver_max_iter_outer,
        'max_iterations_inner': config.interior_struct.zalmoxis.solver_max_iter_inner,
        # JAX+diffrax structure path and Anderson Picard acceleration,
        # both opt-in and defaulting off. See `Zalmoxis.use_jax` /
        # `Zalmoxis.use_anderson` in proteus.config._struct.
        'use_jax': config.interior_struct.zalmoxis.use_jax,
        'use_anderson': config.interior_struct.zalmoxis.use_anderson,
        # outer mass-radius solver dispatch ('newton' default |
        # 'picard'). When 'newton', Zalmoxis uses Newton + brentq
        # bracketing on f(R) = M(R) - M_target instead of the
        # damped-Picard fixed-point loop. Newton requires tight
        # integrator tolerances; we auto-apply newton_relative_tolerance
        # / newton_absolute_tolerance when the Newton path is selected.
        'outer_solver': config.interior_struct.zalmoxis.outer_solver,
        # Newton-specific knobs (newton_max_iter, newton_tol) AND
        # tightened integrator tolerances are passed ONLY when the
        # Newton path is selected. Picard runs see the dict without
        # any Newton keys and without tightened tolerances, so a
        # future Zalmoxis guard against unknown keys would not break
        # Picard callers.
        **(
            {
                'newton_max_iter': (config.interior_struct.zalmoxis.newton_max_iter),
                'newton_tol': config.interior_struct.zalmoxis.newton_tol,
                'relative_tolerance': (
                    config.interior_struct.zalmoxis.newton_relative_tolerance
                ),
                'absolute_tolerance': (
                    config.interior_struct.zalmoxis.newton_absolute_tolerance
                ),
            }
            if config.interior_struct.zalmoxis.outer_solver == 'newton'
            else {}
        ),
    }

def load_zalmoxis_material_dictionaries():
    """Build an EOS registry dict with file paths pointing to FWL_DATA.

    Returns the same dict format as Zalmoxis ``EOS_REGISTRY``, but with
    every ``eos_file`` path resolved under ``FWL_DATA/zalmoxis_eos/``
    instead of ``ZALMOXIS_ROOT/data/``.  This ensures that Zalmoxis,
    when called from PROTEUS, reads EOS data from the central FWL_DATA
    location managed by ``download_zalmoxis_eos()``.

    Returns
    -------
    dict
        Flat dict keyed by EOS identifier string (e.g.
        ``"Seager2007:iron"``, ``"PALEOS:MgSiO3"``, ``"Chabrier:H"``).
    """
    eos_base = get_zalmoxis_eos_dir()

    # Seager2007 paths (also in the EOS_material_properties location)
    seager_dir = eos_base / 'EOS_Seager2007'
    if not seager_dir.exists():
        seager_dir = FWL_DATA_DIR / 'EOS_material_properties' / 'EOS_Seager2007'

    _seager_iron = {'eos_file': str(seager_dir / 'eos_seager07_iron.txt')}
    _seager_silicate = {'eos_file': str(seager_dir / 'eos_seager07_silicate.txt')}
    _seager_water = {'eos_file': str(seager_dir / 'eos_seager07_water.txt')}

    # Wolf & Bower 2018
    wb_dir = eos_base / 'EOS_WolfBower2018_1TPa'
    _wb_melted = {
        'eos_file': str(wb_dir / 'density_melt.dat'),
        'adiabat_grad_file': str(wb_dir / 'adiabat_temp_grad_melt.dat'),
    }
    _wb_solid = {'eos_file': str(wb_dir / 'density_solid.dat')}

    # RTPress 100 TPa
    rt_dir = eos_base / 'EOS_RTPress_melt_100TPa'
    _rt_melted = {
        'eos_file': str(rt_dir / 'density_melt.dat'),
        'adiabat_grad_file': str(rt_dir / 'adiabat_temp_grad_melt.dat'),
    }

    # PALEOS 2-phase MgSiO3 (separate solid/liquid, Zenodo 19680050).
    # 150 pts/decade (default) and 600 pts/decade (highres) variants.
    paleos2ph_dir = eos_base / 'EOS_PALEOS_MgSiO3'
    _paleos2ph_melted = {
        'eos_file': str(paleos2ph_dir / 'paleos_mgsio3_tables_pt_proteus_liquid.dat'),
        'format': 'paleos',
    }
    _paleos2ph_solid = {
        'eos_file': str(paleos2ph_dir / 'paleos_mgsio3_tables_pt_proteus_solid.dat'),
        'format': 'paleos',
    }
    _paleos2ph_melted_highres = {
        'eos_file': str(paleos2ph_dir / 'paleos_mgsio3_tables_pt_proteus_liquid_highres.dat'),
        'format': 'paleos',
    }
    _paleos2ph_solid_highres = {
        'eos_file': str(paleos2ph_dir / 'paleos_mgsio3_tables_pt_proteus_solid_highres.dat'),
        'format': 'paleos',
    }

    # PALEOS unified tables
    _paleos_iron = {
        'eos_file': str(eos_base / 'EOS_PALEOS_iron' / 'paleos_iron_eos_table_pt.dat'),
        'format': 'paleos_unified',
    }
    _paleos_mgsio3 = {
        'eos_file': str(
            eos_base / 'EOS_PALEOS_MgSiO3_unified' / 'paleos_mgsio3_eos_table_pt.dat'
        ),
        'format': 'paleos_unified',
    }
    _paleos_h2o = {
        'eos_file': str(eos_base / 'EOS_PALEOS_H2O' / 'paleos_water_eos_table_pt.dat'),
        'format': 'paleos_unified',
    }

    # Chabrier H/He
    _chabrier_h = {
        'eos_file': str(eos_base / 'EOS_Chabrier2021_HHe' / 'chabrier2021_H.dat'),
        'format': 'paleos_unified',
    }

    # PALEOS-API live tabulation entries. These carry only a GridSpec at
    # build time; the dispatch layer (``zalmoxis.eos.paleos_api_cache``)
    # mutates them in place on first density query to populate ``eos_file``
    # and rewrite ``format`` to the downstream value (``paleos_unified`` or
    # ``paleos``). Cache keys are SHA+grid-hash under
    # ``$ZALMOXIS_ROOT/data/EOS_PALEOS_API/``; cold-cache cost is a one-time
    # generator run (see ``zalmoxis.eos.paleos_api``).
    from zalmoxis.eos.paleos_api import (
        make_default_grid_h2o,
        make_default_grid_iron,
        make_default_grid_mgsio3,
    )

    _paleos_api_iron = {
        'format': 'paleos_api',
        'material': 'iron',
        'grid_spec': make_default_grid_iron(),
    }
    _paleos_api_mgsio3 = {
        'format': 'paleos_api',
        'material': 'mgsio3',
        'grid_spec': make_default_grid_mgsio3(),
    }
    _paleos_api_h2o = {
        'format': 'paleos_api',
        'material': 'h2o',
        'grid_spec': make_default_grid_h2o(),
        'h2o_table_path': None,
    }
    _paleos_api_2ph_mgsio3_melted = {
        'format': 'paleos_api_2phase',
        'material': 'mgsio3',
        'side': 'liquid',
        'grid_spec': make_default_grid_mgsio3(),
    }
    _paleos_api_2ph_mgsio3_solid = {
        'format': 'paleos_api_2phase',
        'material': 'mgsio3',
        'side': 'solid',
        'grid_spec': make_default_grid_mgsio3(),
    }

    return {
        # Seager2007 static
        'Seager2007:iron': {'core': _seager_iron},
        'Seager2007:MgSiO3': {'mantle': _seager_silicate},
        'Seager2007:H2O': {'ice_layer': _seager_water},
        # Wolf & Bower 2018 T-dependent
        'WolfBower2018:MgSiO3': {
            'core': _seager_iron,
            'melted_mantle': _wb_melted,
            'solid_mantle': _wb_solid,
        },
        # RTPress 100 TPa extended melt + WB2018 solid
        'RTPress100TPa:MgSiO3': {
            'core': _seager_iron,
            'melted_mantle': _rt_melted,
            'solid_mantle': _wb_solid,
        },
        # PALEOS 2-phase MgSiO3 (Zenodo 19680050; 150 pts/decade default,
        # 600 pts/decade as -highres variant for sensitivity tests).
        'PALEOS-2phase:MgSiO3': {
            'core': _seager_iron,
            'melted_mantle': _paleos2ph_melted,
            'solid_mantle': _paleos2ph_solid,
        },
        'PALEOS-2phase:MgSiO3-highres': {
            'core': _seager_iron,
            'melted_mantle': _paleos2ph_melted_highres,
            'solid_mantle': _paleos2ph_solid_highres,
        },
        # PALEOS unified
        'PALEOS:iron': _paleos_iron,
        'PALEOS:MgSiO3': _paleos_mgsio3,
        'PALEOS:H2O': _paleos_h2o,
        # PALEOS-API live-tabulated (dispatch populates eos_file on demand)
        'PALEOS-API:iron': _paleos_api_iron,
        'PALEOS-API:MgSiO3': _paleos_api_mgsio3,
        'PALEOS-API:H2O': _paleos_api_h2o,
        'PALEOS-API-2phase:MgSiO3': {
            'core': _seager_iron,
            'melted_mantle': _paleos_api_2ph_mgsio3_melted,
            'solid_mantle': _paleos_api_2ph_mgsio3_solid,
        },
        # Chabrier H/He
        'Chabrier:H': _chabrier_h,
    }

def check_zalmoxis_eos_files(layer_eos_config: dict, mat_dicts: dict) -> None:
    """Fail fast when a selected EOS table file is missing on disk.

    Walks the registry entries selected by ``layer_eos_config`` and
    collects every referenced table path that does not exist, then
    raises one actionable error. Without this check the solver emits
    one read error per shell and ends in a non-convergence failure that
    hides the real cause. Registry entries without an ``eos_file``
    (PALEOS-API live tabulation) generate their tables on demand and
    are skipped.

    Parameters
    ----------
    layer_eos_config : dict
        Per-layer EOS identifier strings (``'core'``, ``'mantle'``, ...).
    mat_dicts : dict
        EOS registry from :func:`load_zalmoxis_material_dictionaries`.

    Raises
    ------
    RuntimeError
        If any selected EOS table file is missing, naming every missing
        path and the command that downloads them.
    """
    missing: set[str] = set()
    for identifier in layer_eos_config.values():
        for component in str(identifier).split('+'):
            entry = mat_dicts.get(_strip_fraction_tokens(component))
            if entry is None:
                # Unknown identifiers fail later with a registry error.
                continue
            # Flat entries carry 'eos_file' directly; nested entries map
            # layer roles (core / melted_mantle / ...) to flat entries.
            subentries = [entry] if 'eos_file' in entry else list(entry.values())
            for sub in subentries:
                if not isinstance(sub, dict):
                    continue
                for field in ('eos_file', 'adiabat_grad_file'):
                    path = sub.get(field)
                    if path and not os.path.isfile(path):
                        missing.add(path)
    if missing:
        listing = '\n  '.join(sorted(missing))
        raise RuntimeError(
            f'Interior EOS table file(s) not found:\n  {listing}\n'
            'Download them with '
            '`proteus get interiordata --config-path <config.toml>`, '
            'or run `proteus start` once without --offline.'
        )

def resolve_2phase_mgsio3_paths(mantle_eos: str, mat_dicts: dict):
    """Return (solid_eos_path, liquid_eos_path) for the 2-phase MgSiO3 tables.

    Selects the API key (``PALEOS-API-2phase:MgSiO3``) when ``mantle_eos``
    is from the PALEOS-API family, otherwise the shipped key
    (``PALEOS-2phase:MgSiO3``). For API entries, calls
    :func:`zalmoxis.eos.paleos_api_cache.resolve_registry_entry` to
    materialise cached ``.dat`` paths in place; this is required because
    :func:`load_zalmoxis_material_dictionaries` rebuilds a fresh registry
    each call, so any earlier in-place mutation is lost.

    Returns
    -------
    tuple[str | None, str | None]
        Absolute filesystem paths if both tables exist, ``(None, None)``
        otherwise. Caller is responsible for treating ``None`` as "no
        2-phase tables available, fall back".
    """
    use_api = mantle_eos.startswith(('PALEOS-API:', 'PALEOS-API-2phase:'))
    if use_api:
        twophase_key = 'PALEOS-API-2phase:MgSiO3'
    elif mantle_eos == 'PALEOS-2phase:MgSiO3-highres':
        twophase_key = 'PALEOS-2phase:MgSiO3-highres'
    else:
        twophase_key = 'PALEOS-2phase:MgSiO3'
    twophase = mat_dicts.get(twophase_key, {})
    if not twophase:
        log.warning(
            'resolve_2phase_mgsio3_paths: registry has no entry for %s '
            '(mantle_eos=%s); 2-phase fallback disabled. This is a '
            'silent-wrong landmine if the caller continues with empty paths.',
            twophase_key,
            mantle_eos,
        )
        return None, None
    if use_api:
        try:
            from zalmoxis.eos.paleos_api_cache import resolve_registry_entry

            resolve_registry_entry(twophase)
        except (ImportError, ModuleNotFoundError) as e:
            log.warning(
                'resolve_2phase_mgsio3_paths: PALEOS-API resolver unavailable '
                '(%s); cannot materialise %s tables.',
                e,
                twophase_key,
            )
            return None, None
    solid_eos = twophase.get('solid_mantle', {}).get('eos_file', '')
    liquid_eos = twophase.get('melted_mantle', {}).get('eos_file', '')
    solid_eos = solid_eos if solid_eos and os.path.isfile(solid_eos) else None
    liquid_eos = liquid_eos if liquid_eos and os.path.isfile(liquid_eos) else None
    return solid_eos, liquid_eos

def load_zalmoxis_solidus_liquidus_functions(mantle_eos: str, config: Config):
    """Loads the solidus and liquidus functions for Zalmoxis based on the mantle EOS.

    Melting curves are needed for two purposes:
    1. Temperature-dependent density in the mushy zone (WolfBower2018, RTPress).
    2. phi(r) blending in VolatileProfile (any EOS with dissolved volatiles).

    For WolfBower2018/RTPress100TPa, loads SPIDER-format P-T files from FWL_DATA.
    For PALEOS unified and PALEOS-2phase, the liquidus comes from the analytic
    Belonoshko+2005 / Fei+2021 curve (Zalmoxis ``'PALEOS-liquidus'``) which is
    the basis Zalmoxis uses for MgSiO3 phase separation, and the solidus is
    derived as T_sol = T_liq * mushy_zone_factor. This keeps the curves used
    for phi-blending and 2-phase nabla_ad consistent with the unified PALEOS
    density-interpolation phase boundaries.

    Parameters
    ----------
    mantle_eos : str
        Mantle EOS string (e.g. ``"WolfBower2018:MgSiO3"``, ``"PALEOS:MgSiO3"``,
        ``"PALEOS-2phase:MgSiO3"``).
    config : Config
        PROTEUS configuration object.

    Returns
    -------
    tuple or None
        (solidus_func, liquidus_func) callable P [Pa] -> T [K], or None.
    """
    _TDEP_PREFIXES = ('WolfBower2018', 'RTPress100TPa')
    if mantle_eos.startswith(_TDEP_PREFIXES):
        return get_zalmoxis_melting_curves(config)

    # PALEOS unified and PALEOS-2phase: both use the same analytic Belonoshko+2005 /
    # Fei+2021 melting curve (`PALEOS-liquidus`) as the basis for MgSiO3 phase
    # separation. The unified path uses it for in-table density interpolation;
    # the 2-phase path uses it inside `_compute_paleos_dtdp` to weight nabla_ad
    # across the solid/liquid blend (mixing.py:_compute_paleos_dtdp). The solidus
    # is derived as `liquidus * mushy_zone_factor` so the mushy band lines up with
    # the unified PALEOS density interpolation. Without these curves, the
    # 2-phase nabla_ad call fails and Zalmoxis structure solve diverges; the
    # unified path falls back to phi=0.5 everywhere in VolatileProfile.
    if mantle_eos.startswith(
        ('PALEOS:', 'PALEOS-2phase:', 'PALEOS-API:', 'PALEOS-API-2phase:')
    ):
        try:
            from zalmoxis.melting_curves import get_solidus_liquidus_functions

            _, liquidus_func = get_solidus_liquidus_functions(
                solidus_id='Stixrude14-solidus',  # required by API but unused; solidus is built below as mushy_zone_factor * liquidus (default 0.8 = the Stixrude 2014 solidus/liquidus ratio)
                liquidus_id='PALEOS-liquidus',
            )
            mzf = config.interior_struct.zalmoxis.mushy_zone_factor
            solidus_func = _make_derived_solidus(liquidus_func, mzf)
            log.info(
                'PALEOS melting curves (%s): liquidus from PALEOS, '
                'solidus = liquidus * %.2f (mushy_zone_factor)',
                mantle_eos,
                mzf,
            )
            return solidus_func, liquidus_func
        except Exception as e:
            log.warning('Could not load PALEOS melting curves: %s', e)
            return None

    return None

def scale_temperature_profile_for_aragog(
    config: Config, mantle_radii: np.ndarray, mantle_temperature_profile: np.ndarray
):
    """Scales the temperature profile obtained from Zalmoxis to match the number of levels required by Aragog.
    Args:
        config (Config): The configuration object containing the configuration parameters.
        mantle_radii (np.ndarray): The radial positions of the mantle layers from Zalmoxis.
        mantle_temperature_profile (np.ndarray): The temperature profile of the mantle layers from Zalmoxis.
    Returns:
        np.ndarray: The scaled temperature profile matching the number of levels in Aragog.
    """

    # Number of levels in Aragog mesh
    mesh_grid_size = config.interior_energetics.num_levels - 1

    # Create new evenly spaced radial positions for Aragog
    radii_to_interpolate = np.linspace(mantle_radii[0], mantle_radii[-1], mesh_grid_size)

    # Cubic interpolation onto the Aragog radial mesh
    cubic_interp = interp1d(mantle_radii, mantle_temperature_profile, kind='cubic')
    return cubic_interp(radii_to_interpolate)

def write_spider_mesh_file(
    outdir: str,
    mantle_radii: np.ndarray,
    mantle_pressure: np.ndarray,
    mantle_density: np.ndarray,
    mantle_gravity: np.ndarray,
    num_basic: int,
) -> str:
    """Write an external mesh file for SPIDER from Zalmoxis mantle profiles.

    Interpolates the Zalmoxis mantle arrays onto uniformly-spaced SPIDER
    basic and staggered nodes, then writes the mesh file in the format
    expected by SPIDER's ``SetMeshFromExternalFile()``.

    Parameters
    ----------
    outdir : str
        PROTEUS output directory (file is written to ``outdir/data/``).
    mantle_radii : np.ndarray
        Radial positions from CMB to surface, ascending [m].
    mantle_pressure : np.ndarray
        Pressure at each radius [Pa].
    mantle_density : np.ndarray
        Density at each radius [kg/m^3].
    mantle_gravity : np.ndarray
        Gravity magnitude at each radius [m/s^2] (positive).
    num_basic : int
        Number of SPIDER basic nodes (shell boundaries).

    Returns
    -------
    str
        Path to the written mesh file.
    """
    num_staggered = num_basic - 1
    R_surf = float(mantle_radii[-1])
    R_cmb = float(mantle_radii[0])

    # Basic nodes: uniform spacing from surface to CMB (descending r)
    r_b = np.linspace(R_surf, R_cmb, num_basic)
    # Staggered nodes: midpoints between consecutive basic nodes
    r_s = 0.5 * (r_b[:-1] + r_b[1:])

    # Interpolate Zalmoxis profiles onto node positions
    # mantle_radii is ascending, np.interp requires ascending xp
    P_b = np.interp(r_b, mantle_radii, mantle_pressure)
    rho_b = np.interp(r_b, mantle_radii, mantle_density)
    g_b = np.interp(r_b, mantle_radii, mantle_gravity)

    P_s = np.interp(r_s, mantle_radii, mantle_pressure)
    rho_s = np.interp(r_s, mantle_radii, mantle_density)
    g_s = np.interp(r_s, mantle_radii, mantle_gravity)

    # Negate gravity for SPIDER convention (inward-pointing, negative)
    g_b = -np.abs(g_b)
    g_s = -np.abs(g_s)

    # Write mesh file
    mesh_path = os.path.join(outdir, 'data', 'spider_mesh.dat')
    with open(mesh_path, 'w') as f:
        f.write(f'# {num_basic} {num_staggered}\n')
        for i in range(num_basic):
            f.write(f'{r_b[i]:.15e} {P_b[i]:.15e} {rho_b[i]:.15e} {g_b[i]:.15e}\n')
        for i in range(num_staggered):
            f.write(f'{r_s[i]:.15e} {P_s[i]:.15e} {rho_s[i]:.15e} {g_s[i]:.15e}\n')

    log.info(
        'Wrote SPIDER mesh file: %s (%d basic + %d staggered nodes)',
        mesh_path,
        num_basic,
        num_staggered,
    )

    return mesh_path

def read_ps_cache_pointer(outdir: str) -> str | None:
    """Return the shared PS-cache table directory recorded for `outdir`.

    Reads the pointer written by :func:`_write_ps_cache_pointer`. Used by
    the resume path to locate PROTEUS_PS_CACHE_DIR tables that live outside
    the run directory.

    Parameters
    ----------
    outdir : str
        The run output directory.

    Returns
    -------
    str or None
        The recorded table directory, or None when the pointer is absent,
        unreadable, or empty.
    """

    pointer = os.path.join(outdir, 'data', PS_CACHE_POINTER_NAME)
    if not os.path.isfile(pointer):
        return None
    try:
        with open(pointer) as f:
            cache_dir = f.read().strip()
    except OSError:
        return None
    return cache_dir or None

def generate_spider_tables(config: Config, outdir: str):
    """Generate P-S EOS tables and phase boundaries from PALEOS data.

    Produces P-S lookup tables for density, temperature, heat capacity,
    thermal expansion, and adiabatic gradient, plus solidus/liquidus phase
    boundaries in S(P) format. These are consumed by the entropy-IC verify
    in Aragog (and by SPIDER if the structure module is SPIDER).

    Supports two PALEOS layouts:

    1. ``paleos_unified`` (e.g. ``PALEOS:MgSiO3``): the structural backbone is
       the single unified P-T table covering both phases plus mushy zone, while
       the per-phase property surfaces are built from the sibling two-phase
       solid + liquid tables when those are present (see the unified branch
       below), so the densities stay resolved across the melting-curve
       discontinuity. The solidus is derived from ``mushy_zone_factor *
       liquidus`` (default 0.8, the Stixrude 2014 solidus/liquidus ratio); the
       liquidus is the analytic PALEOS Belonoshko+2005 / Fei+2021 curve.
    2. ``PALEOS-2phase:<solid>`` (e.g. ``PALEOS-2phase:MgSiO3``): separate
       solid + liquid PALEOS tables. Phase boundaries are sampled at the
       PALEOS-liquidus temperature from each phase table directly. The
       ``mushy_zone_factor`` config value is ignored (no analytic mushy
       zone exists for 2-phase; the gap between solid-table-top and
       liquid-table-bottom defines the latent heat).

    For non-PALEOS EOS types (WolfBower2018, RTPress100TPa), returns None
    and the caller is expected to fall back on pre-existing SPIDER tables.

    Parameters
    ----------
    config : Config
        Configuration object with struct.zalmoxis settings.
    outdir : str
        Output directory. Tables are written to ``outdir/data/spider_eos/``.

    Returns
    -------
    dict or None
        Keys ``'eos_dir'``, ``'solidus_path'``, ``'liquidus_path'`` with
        absolute paths. Returns None if the mantle EOS is not PALEOS.
    """
    from zalmoxis.eos_export import generate_spider_eos_tables, generate_spider_phase_boundaries
    from zalmoxis.melting_curves import get_solidus_liquidus_functions

    mantle_eos = config.interior_struct.zalmoxis.mantle_eos

    # Use FWL_DATA paths (not ZALMOXIS_ROOT) for EOS file lookup
    mat_dicts = load_zalmoxis_material_dictionaries()
    eos_entry = mat_dicts.get(mantle_eos)

    if eos_entry is None:
        log.info(
            'Mantle EOS %s not found in material dictionary; using pre-existing SPIDER tables.',
            mantle_eos,
        )
        return None

    # PALEOS-API live tabulation: materialise cached .dat paths in place so the
    # downstream format / eos_file lookups see concrete paths. No-op for
    # non-PALEOS-API entries. First call on a cold cache triggers generation.
    from zalmoxis.eos.dispatch import _is_paleos_api

    if _is_paleos_api(eos_entry):
        from zalmoxis.eos.paleos_api_cache import resolve_registry_entry

        log.info(
            'PALEOS-API live tabulation: resolving cached tables for %s '
            '(cold-cache build may take up to ~1 h at 600 pts/decade)',
            mantle_eos,
        )
        resolve_registry_entry(eos_entry)

    # Detect format: paleos_unified vs PALEOS-2phase (nested dict).
    is_unified = eos_entry.get('format') == 'paleos_unified'
    is_twophase = (
        'melted_mantle' in eos_entry
        and 'solid_mantle' in eos_entry
        and isinstance(eos_entry.get('melted_mantle'), dict)
        and isinstance(eos_entry.get('solid_mantle'), dict)
    )

    if not (is_unified or is_twophase):
        log.info(
            'Mantle EOS %s is neither PALEOS unified nor PALEOS-2phase; '
            'using pre-existing SPIDER tables.',
            mantle_eos,
        )
        return None

    # Resolve unified file (if present) and 2-phase files (if present).
    eos_file = eos_entry.get('eos_file', '')
    eos_file = eos_file if eos_file and os.path.isfile(eos_file) else None

    if is_twophase:
        solid_eos = eos_entry['solid_mantle'].get('eos_file', '')
        liquid_eos = eos_entry['melted_mantle'].get('eos_file', '')
    else:
        # Unified mantle: also look for sibling 2-phase tables to harden
        # the property surfaces (avoids interpolation across the melting
        # curve discontinuity in the unified table). Use the API-aware
        # helper so PALEOS-API unified runs pull API 2-phase tables
        # rather than silently pulling shipped Zenodo ones.
        # Net effect for PALEOS:MgSiO3: the structure solve uses the unified
        # table, but the per-phase property/density surfaces are taken from
        # these two-phase tables when present. If they are absent the code
        # below falls back to the unified table alone (entropy near the
        # melting curve is then less reliable). The solidus stays synthetic
        # (mushy_zone_factor * liquidus) in both cases.
        solid_eos, liquid_eos = resolve_2phase_mgsio3_paths(mantle_eos, mat_dicts)

    solid_eos = solid_eos if solid_eos and os.path.isfile(solid_eos) else None
    liquid_eos = liquid_eos if liquid_eos and os.path.isfile(liquid_eos) else None

    # For the 2-phase path, the downstream Zalmoxis functions still require
    # an `eos_file` positional (used only to seed default property
    # interpolators that are immediately overridden by the 2-phase ones).
    # Pass the solid table as a sentinel: any valid PALEOS table works.
    if eos_file is None:
        if solid_eos is not None:
            eos_file = solid_eos
        else:
            log.warning(
                'No PALEOS EOS file available for %s '
                '(unified missing and 2-phase incomplete); skipping table gen.',
                mantle_eos,
            )
            return None

    if is_twophase and not (solid_eos and liquid_eos):
        log.warning(
            'PALEOS-2phase entry %s missing solid or liquid file; skipping.',
            mantle_eos,
        )
        return None

    # Phase boundaries: PALEOS-liquidus is the analytic Belonoshko+2005 /
    # Fei+2021 Simon-Glatzel curve. For PALEOS-2phase, mushy_zone_factor=1.0
    # collapses solidus = liquidus and the latent-heat gap is supplied by
    # the entropy difference between solid_table.s(P, T_liq) and
    # liquid_table.s(P, T_liq).
    _, liquidus_func = get_solidus_liquidus_functions(
        solidus_id='Stixrude14-solidus',  # unused, but API requires it
        liquidus_id='PALEOS-liquidus',
    )
    mzf = config.interior_struct.zalmoxis.mushy_zone_factor
    solidus_func = _make_derived_solidus(liquidus_func, mzf)
    if is_twophase:
        log.info(
            'PALEOS-2phase phase boundaries: solidus T = liquidus T '
            '(mushy_zone_factor=%.2f, latent heat from 2-phase tables)',
            mzf,
        )
    else:
        log.info(
            'PALEOS unified phase boundaries: solidus = liquidus * %.2f (mushy_zone_factor)',
            mzf,
        )

    # Determine pressure range from planet mass (higher mass needs wider range)
    mass_tot = config.planet.mass_tot or 1.0
    # P_max for the SPIDER P-S lookup grid. Must cover the actual P_cmb
    # of the planet; the 10 TPa cap covers very massive rocky planets
    # (mass_tot well above 2) without hitting the table edge. See
    # interior_energetics/aragog.py for the matching cap and the
    # comment on EOS / melting-curve calibration ranges.
    P_max = min(1.0e13, 150e9 * mass_tot + 200e9)

    if solid_eos and liquid_eos:
        log.info('Using PALEOS-2phase tables for entropy-IC table generation')

    # Table resolution from config
    nP = config.interior_struct.zalmoxis.lookup_nP
    nS = config.interior_struct.zalmoxis.lookup_nS

    layout = '2phase' if is_twophase else 'unified'
    cache_key = _ps_cache_key(
        P_max=P_max,
        nP=nP,
        nS=nS,
        mzf=mzf,
        layout=layout,
        mantle_eos=mantle_eos,
        eos_file=eos_file,
        solid_eos=solid_eos,
        liquid_eos=liquid_eos,
    )

    # Table location. Default: per-run output/<run>/data/spider_eos. When
    # PROTEUS_PS_CACHE_DIR is set, the directory is keyed by cache_key so that
    # independent runs with the same planet mass, table resolution, and mantle
    # EOS reuse one generated table instead of each rebuilding the slow
    # full-resolution PALEOS P-S table. The cache_key encodes everything that
    # changes the table (P_max, nP, nS, mushy_zone_factor, layout, and the
    # resolved EOS identity), so reuse is exact.
    _ps_cache_root = os.environ.get('PROTEUS_PS_CACHE_DIR')
    if _ps_cache_root:
        _safe_key = cache_key.replace('.', 'p').replace('=', '-').replace('+', '')
        spider_eos_dir = os.path.join(_ps_cache_root, _safe_key)
        os.makedirs(spider_eos_dir, exist_ok=True)
        # The shared-cache tables live outside the run directory, so a resumed
        # run cannot find them via the per-run output/<run>/data/spider_eos
        # path. Leave a pointer so resume can follow the tables to the cache.
        _write_ps_cache_pointer(outdir, spider_eos_dir)
    else:
        spider_eos_dir = os.path.join(outdir, 'data', 'spider_eos')

    # Cache check: skip regeneration if tables exist and pressure range unchanged.
    # The pressure range depends on planet mass, which doesn't change during evolution.
    cache_marker = os.path.join(spider_eos_dir, '.cache_info.txt')
    if os.path.isfile(cache_marker):
        with open(cache_marker) as f:
            existing_key = f.read().strip()
        if existing_key == cache_key:
            # Tables are up to date. File names must match the writer in
            # zalmoxis.eos_export.generate_spider_phase_boundaries, which
            # emits solidus_P-S.dat / liquidus_P-S.dat.
            solidus_path = os.path.join(spider_eos_dir, 'solidus_P-S.dat')
            liquidus_path = os.path.join(spider_eos_dir, 'liquidus_P-S.dat')
            if os.path.isfile(solidus_path) and os.path.isfile(liquidus_path):
                log.info(
                    'Reusing cached PALEOS-derived P-S entropy tables (P_max=%.2e, %dx%d)',
                    P_max,
                    nP,
                    nS,
                )
                return {
                    'eos_dir': spider_eos_dir,
                    'solidus_path': solidus_path,
                    'liquidus_path': liquidus_path,
                }

    # Choose where to generate. For a shared PROTEUS_PS_CACHE_DIR the tables are
    # written into a private staging directory on the same filesystem and then
    # published into spider_eos_dir with per-file atomic renames, so two cluster
    # runs that miss the cache marker at the same moment cannot interleave
    # partial writes into the shared table directory. For a per-run directory
    # there is no sharing, so tables are generated in place as before.
    if _ps_cache_root:
        gen_dir = tempfile.mkdtemp(prefix='.gen-', dir=os.path.dirname(spider_eos_dir))
    else:
        gen_dir = spider_eos_dir

    try:
        # Generate phase boundaries
        log.info(
            'Generating PALEOS-derived P-S phase boundaries (%d P points)...',
            nP,
        )
        generate_spider_phase_boundaries(
            solidus_func=solidus_func,
            liquidus_func=liquidus_func,
            eos_file=eos_file,
            P_range=(1e5, P_max),
            n_P=nP,
            output_dir=gen_dir,
            solid_eos_file=solid_eos,
            liquid_eos_file=liquid_eos,
        )

        # Generate full EOS tables
        log.info(
            'Generating PALEOS-derived P-S EOS tables (%d x %d)...',
            nP,
            nS,
        )
        generate_spider_eos_tables(
            eos_file=eos_file,
            solidus_func=solidus_func,
            liquidus_func=liquidus_func,
            P_range=(1e5, P_max),
            n_P=nP,
            n_S=nS,
            output_dir=gen_dir,
            solid_eos_file=solid_eos,
            liquid_eos_file=liquid_eos,
        )

        # Publish staged tables into the shared cache with atomic renames.
        if gen_dir != spider_eos_dir:
            _publish_ps_tables(gen_dir, spider_eos_dir)
    finally:
        if gen_dir != spider_eos_dir:
            shutil.rmtree(gen_dir, ignore_errors=True)

    # Write the cache marker last, atomically, so a concurrent reader only
    # trusts the directory once every table file is already in place.
    try:
        _atomic_write_text(cache_marker, cache_key)
    except OSError:
        pass

    # File names match the writer in zalmoxis.eos_export, which emits
    # solidus_P-S.dat / liquidus_P-S.dat; use the published locations.
    return {
        'eos_dir': spider_eos_dir,
        'solidus_path': os.path.join(spider_eos_dir, 'solidus_P-S.dat'),
        'liquidus_path': os.path.join(spider_eos_dir, 'liquidus_P-S.dat'),
    }