EOS Package

def _compute_paleos_dtdp(
    pressure, temperature, mat_PALEOS, solidus_func, liquidus_func, interpolation_functions
):
    """Compute dT/dP from PALEOS nabla_ad with phase-aware weighting.

    Uses solid nabla_ad below solidus, liquid nabla_ad above liquidus,
    and melt-fraction-weighted average in the mushy zone.

    Parameters
    ----------
    pressure : float
        Pressure in Pa.
    temperature : float
        Temperature in K.
    mat_PALEOS : dict
        PALEOS material properties dict.
    solidus_func : callable or None
        Solidus melting curve interpolation function.
    liquidus_func : callable or None
        Liquidus melting curve interpolation function.
    interpolation_functions : dict
        Shared interpolation cache.

    Returns
    -------
    float or None
        dT/dP in K/Pa, or None if lookup fails.
    """
    if pressure <= 0 or temperature <= 0:
        return None

    # Determine phase
    T_sol = solidus_func(pressure) if solidus_func is not None else np.nan
    T_liq = liquidus_func(pressure) if liquidus_func is not None else np.nan

    if np.isnan(T_sol) or np.isnan(T_liq) or T_liq <= T_sol:
        # Outside melting curve range or degenerate: use solid table
        nabla = _get_paleos_nabla_ad(
            pressure, temperature, mat_PALEOS, 'solid_mantle', interpolation_functions
        )
    elif temperature <= T_sol:
        # Solid phase
        nabla = _get_paleos_nabla_ad(
            pressure, temperature, mat_PALEOS, 'solid_mantle', interpolation_functions
        )
    elif temperature >= T_liq:
        # Liquid phase
        nabla = _get_paleos_nabla_ad(
            pressure, temperature, mat_PALEOS, 'melted_mantle', interpolation_functions
        )
    else:
        # Mixed phase: melt-fraction-weighted nabla_ad
        phi = (temperature - T_sol) / (T_liq - T_sol)
        nabla_solid = _get_paleos_nabla_ad(
            pressure, temperature, mat_PALEOS, 'solid_mantle', interpolation_functions
        )
        nabla_liquid = _get_paleos_nabla_ad(
            pressure, temperature, mat_PALEOS, 'melted_mantle', interpolation_functions
        )
        if (
            nabla_solid is None and nabla_liquid is None
        ):  # pragma: no cover - both phases undefined; defensive
            return None
        if (
            nabla_solid is None
        ):  # pragma: no cover - solid-side phase undefined at table edge; defensive
            nabla = nabla_liquid
        elif (
            nabla_liquid is None
        ):  # pragma: no cover - liquid-side phase undefined at table edge; defensive
            nabla = nabla_solid
        else:
            nabla = (1.0 - phi) * nabla_solid + phi * nabla_liquid

    if nabla is None or nabla <= 0:  # pragma: no cover - non-positive nabla_ad; defensive
        return None

    # Convert: dT/dP = nabla_ad * T / P
    return nabla * temperature / pressure

def _ensure_unified_cache(eos_file, interpolation_functions):
    """Ensure a unified PALEOS table is loaded into the interpolation cache.

    Tries loading from a binary pickle cache first (fast), then falls back
    to the text table (slow). Saves a pickle cache after the first text load.

    Parameters
    ----------
    eos_file : str
        Path to the unified PALEOS table file.
    interpolation_functions : dict
        Shared interpolation cache.

    Returns
    -------
    dict
        The cache entry for this file.
    """
    if eos_file not in interpolation_functions:
        import pickle

        cache_path = eos_file.replace('.dat', '.pkl')
        try:
            with open(cache_path, 'rb') as f:
                interpolation_functions[eos_file] = pickle.load(f)
            logger.debug('Loaded PALEOS cache from %s', cache_path)
        except (FileNotFoundError, EOFError, pickle.UnpicklingError):
            logger.info('Loading PALEOS table from text: %s', eos_file)
            interpolation_functions[eos_file] = load_paleos_unified_table(eos_file)
            # Save binary cache for next time
            try:
                with open(cache_path, 'wb') as f:
                    pickle.dump(interpolation_functions[eos_file], f, protocol=4)
                logger.info('Saved PALEOS binary cache: %s', cache_path)
            except OSError:
                logger.debug('Could not save cache to %s', cache_path)

    return interpolation_functions[eos_file]

def _fast_bilinear(log_p, log_t, grid, cached):
    """O(1) bilinear interpolation on a log-uniform grid (scalar).

    Optimized for the hot path: called ~2M times per structure solve.
    Uses direct floor indexing (no rounding, no max/min builtins).

    Parameters
    ----------
    log_p : float
        log10 of pressure, already clamped to grid bounds.
    log_t : float
        log10 of temperature, already clamped to valid range.
    grid : numpy.ndarray
        2D array of values, shape (n_p, n_t).
    cached : dict
        Cache entry with grid metadata.

    Returns
    -------
    float
        Interpolated value. NaN if any corner is NaN.
    """
    n_p_m2 = cached['n_p'] - 2
    n_t_m2 = cached['n_t'] - 2

    if n_p_m2 < 0 or n_t_m2 < 0:
        return grid[0, 0]

    # O(1) lower-bound index via floor division
    fp = (log_p - cached['logp_min']) / cached['dlog_p']
    ft = (log_t - cached['logt_min']) / cached['dlog_t']

    ip = int(fp)
    it = int(ft)

    # Clamp to valid range without max/min builtins
    if ip < 0:
        ip = 0
    elif ip > n_p_m2:
        ip = n_p_m2
    if it < 0:
        it = 0
    elif it > n_t_m2:
        it = n_t_m2

    # Fractional position within cell
    dp = fp - ip
    dt = ft - it
    if dp < 0.0:
        dp = 0.0
    elif dp > 1.0:
        dp = 1.0
    if dt < 0.0:
        dt = 0.0
    elif dt > 1.0:
        dt = 1.0

    # Bilinear blend
    omdp = 1.0 - dp
    omdt = 1.0 - dt
    return (
        grid[ip, it] * omdp * omdt
        + grid[ip, it + 1] * omdp * dt
        + grid[ip + 1, it] * dp * omdt
        + grid[ip + 1, it + 1] * dp * dt
    )

def _get_paleos_nabla_ad(pressure, temperature, material_dict, phase, interpolation_functions):
    """Look up nabla_ad from a PALEOS cache entry.

    Parameters
    ----------
    pressure : float
        Pressure in Pa.
    temperature : float
        Temperature in K.
    material_dict : dict
        Material properties dict (e.g. ``material_properties_iron_PALEOS_silicate_planets``).
    phase : str
        ``'solid_mantle'`` or ``'melted_mantle'``.
    interpolation_functions : dict
        Shared interpolation cache.

    Returns
    -------
    float or None
        Dimensionless adiabatic gradient nabla_ad, or None if lookup fails.
    """
    props = material_dict[phase]
    eos_file = props['eos_file']

    # Ensure the PALEOS table is loaded into the cache
    if eos_file not in interpolation_functions:
        interpolation_functions[eos_file] = load_paleos_table(eos_file)

    cached = interpolation_functions[eos_file]
    p_min, p_max = cached['p_min'], cached['p_max']

    # Clamp pressure to global bounds
    p_clamped = np.clip(pressure, p_min, p_max)
    log_p = np.log10(p_clamped)
    log_t = np.log10(max(temperature, 1.0))

    # Per-cell clamping: restrict T to the valid data range at this P
    log_t_clamped, was_clamped = _paleos_clamp_temperature(log_p, log_t, cached)
    if was_clamped and eos_file not in _paleos_clamp_warned:
        _paleos_clamp_warned.add(eos_file)
        logger.warning(
            f'PALEOS per-cell clamping active for nabla_ad in '
            f'{os.path.basename(eos_file)}: '
            f'T={temperature:.0f} K clamped to {10.0**log_t_clamped:.0f} K '
            f'at P={pressure:.2e} Pa.'
        )

    val = float(cached['nabla_ad_interp']((log_p, log_t_clamped)))

    # Nearest-neighbor fallback for ragged boundary
    if not np.isfinite(val):
        val = float(cached['nabla_ad_nn']((log_p, log_t_clamped)))

    if np.isfinite(val):
        return val
    return None

def _get_paleos_unified_nabla_ad(pressure, temperature, material_dict, interpolation_functions):
    """Look up nabla_ad from a unified PALEOS cache entry.

    Parameters
    ----------
    pressure : float
        Pressure in Pa.
    temperature : float
        Temperature in K.
    material_dict : dict
        Material properties dict with 'eos_file' key.
    interpolation_functions : dict
        Shared interpolation cache.

    Returns
    -------
    float or None
        Dimensionless adiabatic gradient, or None if lookup fails.
    """
    eos_file = material_dict['eos_file']
    cached = _ensure_unified_cache(eos_file, interpolation_functions)

    p_clamped = np.clip(pressure, cached['p_min'], cached['p_max'])
    log_p = np.log10(p_clamped)
    log_t = np.log10(max(temperature, 1.0))

    log_t_clamped, was_clamped = _paleos_clamp_temperature(log_p, log_t, cached)
    if was_clamped and eos_file not in _paleos_clamp_warned:
        _paleos_clamp_warned.add(eos_file)
        logger.warning(
            f'PALEOS unified per-cell clamping active for nabla_ad in '
            f'{os.path.basename(eos_file)}: '
            f'T={temperature:.0f} K clamped to {10.0**log_t_clamped:.0f} K '
            f'at P={pressure:.2e} Pa.'
        )

    val = _fast_bilinear(log_p, log_t_clamped, cached['nabla_ad_grid'], cached)
    if not np.isfinite(val):
        val = float(cached['nabla_ad_nn']((log_p, log_t_clamped)))

    return val if np.isfinite(val) else None

def _paleos_clamp_temperature(log_p, log_t, cached):
    """Clamp log10(T) to the per-pressure valid range of a PALEOS table.

    Uses O(1) index computation on the log-uniform pressure grid instead
    of np.interp (which does binary search on 1000+ elements). This is
    called millions of times in the scalar ODE path.

    Parameters
    ----------
    log_p : float
        log10 of pressure in Pa (already clamped to table bounds).
    log_t : float
        log10 of temperature in K.
    cached : dict
        PALEOS cache entry.

    Returns
    -------
    float
        Clamped log10(T).
    bool
        True if clamping was applied.
    """
    lt_min = cached['logt_valid_min']
    lt_max = cached['logt_valid_max']

    # O(1) path for log-uniform grids (has dlog_p metadata)
    if 'dlog_p' in cached:
        fp = (log_p - cached['logp_min']) / cached['dlog_p']
        n_p_m1 = cached['n_p'] - 2
        ip = int(fp)
        if ip < 0:
            ip = 0
        elif ip > n_p_m1:
            ip = n_p_m1
        frac = fp - ip
        if frac < 0.0:
            frac = 0.0
        elif frac > 1.0:
            frac = 1.0
        local_tmin = lt_min[ip] + frac * (lt_min[ip + 1] - lt_min[ip])
        local_tmax = lt_max[ip] + frac * (lt_max[ip + 1] - lt_max[ip])
    else:
        # Fallback for legacy/test cache dicts without grid metadata
        ulp = cached['unique_log_p']
        local_tmin = float(np.interp(log_p, ulp, lt_min))
        local_tmax = float(np.interp(log_p, ulp, lt_max))

    # Guard: NaN bounds near table edges
    if not (local_tmin == local_tmin and local_tmax == local_tmax):  # fast NaN check
        return log_t, False

    if log_t < local_tmin:
        return local_tmin, True
    elif log_t > local_tmax:
        return local_tmax, True
    return log_t, False

def calculate_density(
    pressure,
    material_dictionaries,
    layer_eos,
    temperature,
    solidus_func,
    liquidus_func,
    interpolation_functions=None,
    mushy_zone_factor=1.0,
):
    """Calculate density for a single layer given its EOS identifier.

    Parameters
    ----------
    pressure : float
        Pressure at which to evaluate the EOS, in Pa.
    material_dictionaries : dict
        EOS registry dict keyed by EOS identifier string
        (from ``eos_properties.EOS_REGISTRY``).
    layer_eos : str
        Per-layer EOS identifier, for example ``"Seager2007:iron"``,
        ``"WolfBower2018:MgSiO3"``, ``"PALEOS:iron"``, or ``"Analytic:iron"``.
    temperature : float
        Temperature at which to evaluate the EOS, in K.
    solidus_func : callable or None
        Interpolation function for the solidus melting curve.
    liquidus_func : callable or None
        Interpolation function for the liquidus melting curve.
    interpolation_functions : dict, optional
        Cache of interpolation functions used to avoid redundant loading.
    mushy_zone_factor : float, optional
        Cryoscopic depression factor for unified PALEOS tables. 1.0 = no
        mushy zone (sharp boundary from table). Default 1.0.

    Returns
    -------
    float or None
        Density in kg/m^3, or ``None`` on failure.

    Raises
    ------
    ValueError
        If ``layer_eos`` is not recognized.
    """
    if interpolation_functions is None:
        interpolation_functions = {}

    # Analytic EOS: no material dict needed
    if layer_eos.startswith('Analytic:'):
        material_key = layer_eos.split(':', 1)[1]
        return get_analytic_density(pressure, material_key)

    # Vinet (Rose-Vinet) EOS: no material dict needed
    if layer_eos.startswith('Vinet:'):
        material_key = layer_eos.split(':', 1)[1]
        return get_vinet_density(pressure, material_key)

    # Look up material properties from the registry
    mat = material_dictionaries.get(layer_eos)
    if mat is None:
        raise ValueError(f"Unknown layer EOS '{layer_eos}'.")

    # PALEOS-API live tabulation: rewrite the entry in place on first access
    # so the remaining dispatch goes through the normal paleos_unified / paleos
    # code paths. A sticky `_api_resolved` flag short-circuits the full
    # `_is_paleos_api(mat)` dispatch chain on every subsequent RHS call;
    # without it the dispatch check itself dominated the standalone solve
    # (~25M calls per main(), ~6 % of self time in cProfile).
    if '_api_resolved' not in mat:
        if _is_paleos_api(mat):
            from .paleos_api_cache import resolve_registry_entry

            resolve_registry_entry(mat)
        mat['_api_resolved'] = True

    # Unified PALEOS tables (single file per material, all phases included)
    if mat.get('format') == 'paleos_unified':
        return get_paleos_unified_density(
            pressure, temperature, mat, mushy_zone_factor, interpolation_functions
        )

    # T-dependent EOS with separate solid/liquid tables (WB2018, RTPress, PALEOS-2phase)
    if 'melted_mantle' in mat:
        return get_Tdep_density(
            pressure, temperature, mat, solidus_func, liquidus_func, interpolation_functions
        )

    # Seager2007 static EOS (1D P-rho tables)
    # Determine the layer key from the material dict
    for layer_key in ('core', 'mantle', 'ice_layer'):
        if layer_key in mat:
            return get_tabulated_eos(
                pressure, mat, layer_key, interpolation_functions=interpolation_functions
            )

    raise ValueError(f"Cannot determine layer key for EOS '{layer_eos}'.")

def calculate_density_batch(
    pressures,
    temperatures,
    material_dictionaries,
    layer_eos,
    solidus_func,
    liquidus_func,
    interpolation_functions,
    mushy_zone_factor=1.0,
):
    """Vectorized density lookup for a batch of (P, T) points sharing one EOS.

    For unified PALEOS tables, uses the vectorized interpolator path.
    For other EOS types, falls back to scalar calculate_density per point.

    Parameters
    ----------
    pressures : numpy.ndarray
        1D array of pressures in Pa.
    temperatures : numpy.ndarray
        1D array of temperatures in K.
    material_dictionaries : dict
        EOS registry.
    layer_eos : str
        EOS identifier string.
    solidus_func : callable or None
        Solidus melting curve function.
    liquidus_func : callable or None
        Liquidus melting curve function.
    interpolation_functions : dict
        Interpolation cache.
    mushy_zone_factor : float
        Cryoscopic depression factor.

    Returns
    -------
    numpy.ndarray
        1D array of densities in kg/m^3. NaN where lookup fails.
    """
    if interpolation_functions is None:
        interpolation_functions = {}

    mat = material_dictionaries.get(layer_eos)
    # See calculate_density for the _api_resolved sticky-flag rationale.
    if mat is not None and '_api_resolved' not in mat:
        if _is_paleos_api(mat):
            from .paleos_api_cache import resolve_registry_entry

            resolve_registry_entry(mat)
        mat['_api_resolved'] = True
    if mat is not None and mat.get('format') == 'paleos_unified':
        return get_paleos_unified_density_batch(
            pressures, temperatures, mat, mushy_zone_factor, interpolation_functions
        )

    # Fallback: scalar loop for non-unified EOS types
    n = len(pressures)
    result = np.full(n, np.nan)
    for i in range(n):
        val = calculate_density(
            pressures[i],
            material_dictionaries,
            layer_eos,
            temperatures[i],
            solidus_func,
            liquidus_func,
            interpolation_functions,
            mushy_zone_factor,
        )
        result[i] = val if val is not None else np.nan
    return result

def calculate_temperature_profile(
    radii,
    temperature_mode,
    surface_temperature,
    center_temperature,
    input_dir,
    temp_profile_file,
    cmb_temperature=None,
):
    """Return a callable temperature profile for a planetary interior model.

    Parameters
    ----------
    radii : array_like
        Radial grid of the planet, in m.
    temperature_mode : {"isothermal", "linear", "prescribed", "adiabatic", "adiabatic_from_cmb"}
        Temperature profile mode.

        - ``"isothermal"``: constant temperature equal to
          ``surface_temperature``.
        - ``"linear"``: linear profile from ``center_temperature`` at
          ``r = 0`` to ``surface_temperature`` at the surface.
        - ``"prescribed"``: read the temperature profile from a text file.
        - ``"adiabatic"``: return a linear profile as an initial guess; the
          actual adiabat is computed elsewhere in the main iteration loop
          and is anchored at ``surface_temperature``.
        - ``"adiabatic_from_cmb"``: same as ``"adiabatic"`` but the actual
          adiabat (computed elsewhere) is anchored at ``cmb_temperature``
          at the core-mantle boundary and integrated outward to the
          surface. The initial guess returned here is still linear, with
          ``cmb_temperature`` used as the deep anchor in place of
          ``center_temperature`` so the first density iteration sees a
          reasonable mantle T(r).
    surface_temperature : float
        Temperature at the surface, in K.
    center_temperature : float
        Temperature at the center, in K. Used for ``"linear"`` and
        ``"adiabatic"`` modes.
    cmb_temperature : float, optional
        Temperature at the core-mantle boundary, in K. Required for
        ``"adiabatic_from_cmb"`` mode; ignored otherwise.
    input_dir : str or path-like
        Directory containing the prescribed temperature profile file.
    temp_profile_file : str
        Name of the file containing the prescribed temperature profile from
        center to surface. Required when ``temperature_mode="prescribed"``.

    Returns
    -------
    callable
        Function of radius returning temperature in K. The callable accepts a
        scalar or array-like radius and returns a float or NumPy array.

    Raises
    ------
    ValueError
        If ``temperature_mode="prescribed"`` and the file does not exist.
    ValueError
        If the prescribed temperature profile length does not match
        ``radii``.
    ValueError
        If ``temperature_mode`` is not recognized.
    """
    radii = np.array(radii)

    if temperature_mode == 'isothermal':
        return lambda r: np.full_like(r, surface_temperature, dtype=float)

    elif temperature_mode == 'linear':
        return lambda r: (
            surface_temperature
            + (center_temperature - surface_temperature) * (1 - np.array(r) / radii[-1])
        )

    elif temperature_mode == 'prescribed':
        temp_profile_path = os.path.join(input_dir, temp_profile_file)
        if not os.path.exists(temp_profile_path):
            raise ValueError(
                "Temperature profile file must be provided and exist for 'prescribed' temperature mode."
            )
        temp_profile = np.loadtxt(temp_profile_path)
        if len(temp_profile) != len(radii):
            raise ValueError('Temperature profile length does not match radii length.')
        # Vectorized interpolation for arbitrary radius points
        return lambda r: np.interp(np.array(r), radii, temp_profile)

    elif temperature_mode == 'adiabatic':
        # Return linear profile as initial guess for the first outer iteration.
        # The actual adiabat is computed in main() using P(r), g(r) from the solver.
        return lambda r: (
            surface_temperature
            + (center_temperature - surface_temperature) * (1 - np.array(r) / radii[-1])
        )

    elif temperature_mode == 'adiabatic_from_cmb':
        # Initial guess: linear from surface to a deep anchor that uses
        # cmb_temperature when provided (otherwise falls back to
        # center_temperature). The CMB index is unknown at this point in
        # the solve, so anchor at r=0 with the CMB value as a safe
        # over-estimate; the actual upward adiabat from CMB is computed in
        # main() once the converged structure exposes the CMB index.
        deep_anchor = (
            cmb_temperature
            if cmb_temperature is not None and cmb_temperature > 0
            else center_temperature
        )
        return lambda r: (
            surface_temperature
            + (deep_anchor - surface_temperature) * (1 - np.array(r) / radii[-1])
        )

    else:
        raise ValueError(
            f"Unknown temperature mode '{temperature_mode}'. "
            f"Valid options: 'isothermal', 'linear', 'prescribed', 'adiabatic', "
            f"'adiabatic_from_cmb'."
        )

def compute_adiabatic_temperature(
    radii,
    pressure,
    mass_enclosed,
    surface_temperature,
    cmb_mass,
    core_mantle_mass,
    layer_mixtures,
    material_dictionaries,
    interpolation_functions=None,
    solidus_func=None,
    liquidus_func=None,
    mushy_zone_factors=None,
    condensed_rho_min=CONDENSED_RHO_MIN_DEFAULT,
    condensed_rho_scale=CONDENSED_RHO_SCALE_DEFAULT,
    binodal_T_scale=50.0,
    anchor='surface',
    cmb_temperature=None,
):
    """Compute an adiabatic temperature profile using native EOS gradient tables.

    Supports single-material and multi-material (volume-additive) layers.
    For multi-material layers, nabla_ad is mass-fraction-weighted across
    components.

    Parameters
    ----------
    radii : numpy.ndarray
        Radial grid in ascending order from center to surface, in m.
    pressure : numpy.ndarray
        Pressure at each radius, in Pa.
    mass_enclosed : numpy.ndarray
        Enclosed mass at each radius, in kg.
    surface_temperature : float
        Temperature at the surface, in K. Used as the anchor when
        ``anchor='surface'``; carried through as-is for the core when
        ``anchor='cmb'`` (no integration is done below the CMB).
    cmb_mass : float
        Core-mantle boundary mass, in kg.
    core_mantle_mass : float
        Total core-plus-mantle mass, in kg.
    layer_mixtures : dict
        Per-layer LayerMixture objects.
    material_dictionaries : dict
        EOS registry dict keyed by EOS identifier string.
    interpolation_functions : dict, optional
        Cache of interpolation functions used to avoid redundant loading.
    solidus_func : callable, optional
        Interpolation function for the solidus melting curve.
    liquidus_func : callable, optional
        Interpolation function for the liquidus melting curve.
    mushy_zone_factors : dict or float or None, optional
        Per-EOS mushy zone factors. Dict keyed by EOS name, a single
        float (applied to all), or None (default 1.0 for all).
    condensed_rho_min : float, optional
        Sigmoid center for phase-aware suppression (kg/m^3). Default 322.
    condensed_rho_scale : float, optional
        Sigmoid width for phase-aware suppression (kg/m^3). Default 50.
    binodal_T_scale : float, optional
        Binodal sigmoid width in K for H2 miscibility suppression.
        Default 50.
    anchor : {'surface', 'cmb'}, optional
        Where the adiabat is anchored. ``'surface'`` (default) anchors
        ``T[n-1] = surface_temperature`` and integrates inward to the
        center, the original behaviour. ``'cmb'`` anchors
        ``T[i_cmb] = cmb_temperature`` at the first mantle shell and
        integrates outward to the surface (mantle only); shells below
        the CMB carry the surface anchor through, since the core
        adiabat is decoupled and the energetics solver computes T_core
        independently.
    cmb_temperature : float, optional
        Temperature at the core-mantle boundary, in K. Required when
        ``anchor='cmb'``.

    Returns
    -------
    numpy.ndarray
        Temperature at each radial point, in K.
    """
    from ..mixing import get_mixed_nabla_ad
    from ..structure_model import get_layer_mixture

    if interpolation_functions is None:
        interpolation_functions = {}

    from ..mixing import any_component_is_tdep

    if not any_component_is_tdep(layer_mixtures):
        raise ValueError(
            'Adiabatic temperature mode requires at least one T-dependent EOS '
            'layer, but none found. '
            "Use 'linear' or 'isothermal' instead."
        )

    n = len(radii)
    T = np.zeros(n)

    if anchor == 'cmb':
        if cmb_temperature is None or cmb_temperature <= 0:
            raise ValueError(
                f"anchor='cmb' requires a positive cmb_temperature, got {cmb_temperature}."
            )

        # Find the first mantle shell (index where mass_enclosed >= cmb_mass).
        cmb_index = int(np.searchsorted(mass_enclosed, cmb_mass))
        cmb_index = max(1, min(cmb_index, n - 1))

        # Carry the surface anchor through the core (no integration);
        # the energetics solver handles T_core via core_heatcap and the
        # Bower+2018 adiabatic ratio. Anchor the mantle at CMB.
        T[:cmb_index] = surface_temperature
        T[cmb_index] = cmb_temperature

        # Integrate upward from CMB to surface (i = cmb_index .. n-1).
        for i in range(cmb_index + 1, n):
            mixture = get_layer_mixture(
                mass_enclosed[i],
                cmb_mass,
                core_mantle_mass,
                layer_mixtures,
            )

            if not mixture.has_tdep():
                T[i] = T[i - 1]
                continue

            P_eval = pressure[i - 1]
            T_eval = T[i - 1]
            dP = pressure[i] - pressure[i - 1]  # negative going outward

            if P_eval < 1e5 or pressure[i] < 1e5:
                T[i] = T[i - 1]
                continue

            nabla = get_mixed_nabla_ad(
                P_eval,
                T_eval,
                mixture,
                material_dictionaries,
                interpolation_functions,
                solidus_func,
                liquidus_func,
                mushy_zone_factors,
                condensed_rho_min,
                condensed_rho_scale,
                binodal_T_scale,
            )

            if nabla is not None and nabla > 0 and P_eval > 0 and T_eval > 0 and dP < 0:
                dtdp = nabla * T_eval / P_eval
                T_new = T_eval + dtdp * dP  # dP<0 => T cools outward
                T[i] = max(min(T_new, 100000.0), 100.0)
            else:
                T[i] = T_eval

        return T

    # anchor == 'surface' (original behaviour): integrate inward from
    # T[n-1] = surface_temperature to the center.
    T[n - 1] = surface_temperature

    for i in range(n - 2, -1, -1):
        mixture = get_layer_mixture(
            mass_enclosed[i],
            cmb_mass,
            core_mantle_mass,
            layer_mixtures,
        )

        if not mixture.has_tdep():
            T[i] = T[i + 1]
            continue

        P_eval = pressure[i + 1]
        T_eval = T[i + 1]
        dP = pressure[i] - pressure[i + 1]

        # Skip at very low P to prevent T/P divergence
        if P_eval < 1e5 or pressure[i] < 1e5:
            T[i] = T[i + 1]
            continue

        # Get nabla_ad (handles single and multi-component mixtures)
        nabla = get_mixed_nabla_ad(
            P_eval,
            T_eval,
            mixture,
            material_dictionaries,
            interpolation_functions,
            solidus_func,
            liquidus_func,
            mushy_zone_factors,
            condensed_rho_min,
            condensed_rho_scale,
            binodal_T_scale,
        )

        if nabla is not None and nabla > 0 and P_eval > 0 and T_eval > 0 and dP > 0:
            dtdp = nabla * T_eval / P_eval
            T_new = T_eval + dtdp * dP
            # Cap temperature to the PALEOS table maximum (100,000 K).
            # Without this cap, numerical noise at low P or in mixed
            # compositions can cause runaway T during the adiabat
            # integration, which then pollutes the T(P) interpolator
            # and prevents the Brent solver from bracketing.
            T[i] = min(T_new, 100000.0)
        else:
            T[i] = T_eval

    return T

def create_pressure_density_files(
    outer_iter, inner_iter, pressure_iter, radii, pressure, density
):
    """Append pressure and density profiles to output files for a given iteration.

    Parameters
    ----------
    outer_iter : int
        Current outer iteration index.
    inner_iter : int
        Current inner iteration index.
    pressure_iter : int
        Current pressure iteration index.
    radii : numpy.ndarray
        Radial positions, in m.
    pressure : numpy.ndarray
        Pressure values corresponding to ``radii``, in Pa.
    density : numpy.ndarray
        Density values corresponding to ``radii``, in kg/m^3.

    Returns
    -------
    None
    """

    output_dir = os.path.join(get_zalmoxis_root(), 'output')
    os.makedirs(output_dir, exist_ok=True)
    pressure_file = os.path.join(output_dir, 'pressure_profiles.txt')
    density_file = os.path.join(output_dir, 'density_profiles.txt')

    # Only delete the files once at the beginning of the run
    if outer_iter == 0 and inner_iter == 0 and pressure_iter == 0:
        for file_path in [pressure_file, density_file]:
            if os.path.exists(file_path):
                os.remove(file_path)

    # Append current iteration's pressure profile to file
    with open(pressure_file, 'a') as f:
        f.write(f'# Pressure iteration {pressure_iter}\n')
        np.savetxt(f, np.column_stack((radii, pressure)), header='radius pressure', comments='')
        f.write('\n')

    # Append current iteration's density profile to file
    with open(density_file, 'a') as f:
        f.write(f'# Pressure iteration {pressure_iter}\n')
        np.savetxt(f, np.column_stack((radii, density)), header='radius density', comments='')
        f.write('\n')

def get_Tdep_density(
    pressure,
    temperature,
    material_properties_iron_Tdep_silicate_planets,
    solidus_func,
    liquidus_func,
    interpolation_functions=None,
):
    """Compute mantle density for a temperature-dependent EOS with phase changes.

    Parameters
    ----------
    pressure : float
        Pressure at which to evaluate the EOS, in Pa.
    temperature : float
        Temperature at which to evaluate the EOS, in K.
    material_properties_iron_Tdep_silicate_planets : dict
        Dictionary containing temperature-dependent material properties for
        the MgSiO3 EOS.
    solidus_func : callable
        Interpolation function for the solidus melting curve.
    liquidus_func : callable
        Interpolation function for the liquidus melting curve.
    interpolation_functions : dict, optional
        Cache of interpolation functions used to avoid redundant loading of
        EOS tables.

    Returns
    -------
    float or None
        Density in kg/m^3. Returns ``None`` if the density cannot be evaluated.

    Raises
    ------
    ValueError
        If ``solidus_func`` or ``liquidus_func`` is not provided.

    Notes
    -----
    The mantle phase is determined by comparing the input temperature to the
    solidus and liquidus temperatures at the given pressure.

    In the mixed-phase region, the density is computed using linear melt
    fraction and volume additivity.
    """

    if interpolation_functions is None:
        interpolation_functions = {}

    if solidus_func is None or liquidus_func is None:
        raise ValueError(
            'solidus_func and liquidus_func must be provided for WolfBower2018:MgSiO3 EOS.'
        )

    T_sol = solidus_func(pressure)
    T_liq = liquidus_func(pressure)

    # Pressure outside melting curve range -- default to solid phase
    if np.isnan(T_sol) or np.isnan(T_liq):
        logger.debug(
            f'Melting curve undefined at P={pressure:.2e} Pa. Defaulting to solid phase.'
        )
        return get_tabulated_eos(
            pressure,
            material_properties_iron_Tdep_silicate_planets,
            'solid_mantle',
            temperature,
            interpolation_functions,
        )

    if temperature <= T_sol:
        # Solid phase
        rho = get_tabulated_eos(
            pressure,
            material_properties_iron_Tdep_silicate_planets,
            'solid_mantle',
            temperature,
            interpolation_functions,
        )
        return rho

    elif temperature >= T_liq:
        # Liquid phase
        rho = get_tabulated_eos(
            pressure,
            material_properties_iron_Tdep_silicate_planets,
            'melted_mantle',
            temperature,
            interpolation_functions,
        )
        return rho

    else:
        # Mixed phase: linear melt fraction between solidus and liquidus.
        # Guard against degenerate melting curves where T_liq == T_sol.
        if T_liq <= T_sol:
            return get_tabulated_eos(
                pressure,
                material_properties_iron_Tdep_silicate_planets,
                'melted_mantle',
                temperature,
                interpolation_functions,
            )
        # Smoothstep ramp s = x*x*(3 - 2*x) replaces the linear melt
        # fraction so d(1/rho)/dT vanishes at T=T_sol and T=T_liq,
        # eliminating the lever-rule kink that drives inner Picard
        # plateau on hot mushy profiles. Midpoint s(0.5)=0.5 = linear,
        # so the in-mushy density profile is preserved away from the
        # boundaries. Mirror change in jax_eos/tdep.py for parity.
        frac_melt_raw = (temperature - T_sol) / (T_liq - T_sol)
        x = max(0.0, min(1.0, frac_melt_raw))
        frac_melt = x * x * (3.0 - 2.0 * x)
        rho_solid = get_tabulated_eos(
            pressure,
            material_properties_iron_Tdep_silicate_planets,
            'solid_mantle',
            temperature,
            interpolation_functions,
        )
        rho_liquid = get_tabulated_eos(
            pressure,
            material_properties_iron_Tdep_silicate_planets,
            'melted_mantle',
            temperature,
            interpolation_functions,
        )
        # Guard against out-of-bounds pressure returning None
        if rho_solid is None or rho_liquid is None:
            return None
        # Calculate mixed density by volume additivity
        specific_volume_mixed = frac_melt * (1 / rho_liquid) + (1 - frac_melt) * (1 / rho_solid)
        rho_mixed = 1 / specific_volume_mixed
        return rho_mixed

def get_Tdep_material(pressure, temperature, solidus_func, liquidus_func):
    """Determine the mantle phase for a temperature-dependent EOS.

    Parameters
    ----------
    pressure : float or array_like
        Pressure in Pa. May be a scalar or an array.
    temperature : float or array_like
        Temperature in K. May be a scalar or an array.
    solidus_func : callable
        Interpolation function for the solidus melting curve.
    liquidus_func : callable
        Interpolation function for the liquidus melting curve.

    Returns
    -------
    str or numpy.ndarray
        Material phase label(s). Possible values are ``"solid_mantle"``,
        ``"mixed_mantle"``, and ``"melted_mantle"``. Returns a string for
        scalar inputs and a NumPy array of strings for array inputs.
    """

    # Define per-point evaluation
    def evaluate_phase(P, T):
        T_sol = solidus_func(P)
        T_liq = liquidus_func(P)
        # Guard against degenerate melting curves where T_liq == T_sol
        if T_liq <= T_sol:
            return 'melted_mantle' if T >= T_sol else 'solid_mantle'
        frac_melt = (T - T_sol) / (T_liq - T_sol)
        if frac_melt < 0:
            return 'solid_mantle'
        elif frac_melt <= 1.0:
            return 'mixed_mantle'
        else:
            return 'melted_mantle'

    # Vectorize function for array support
    vectorized_eval = np.vectorize(evaluate_phase, otypes=[str])

    # Apply depending on input type
    if np.isscalar(pressure) and np.isscalar(temperature):
        return evaluate_phase(pressure, temperature)
    else:
        return vectorized_eval(pressure, temperature)

def get_paleos_unified_density(
    pressure, temperature, material_dict, mushy_zone_factor, interpolation_functions
):
    """Look up density from a unified PALEOS table.

    When ``mushy_zone_factor == 1.0`` (no mushy zone), the density is read
    directly from the table (the stable phase at each (P, T) is already
    encoded). When ``mushy_zone_factor < 1.0``, a synthetic solidus is
    derived as ``T_sol = T_liq * mushy_zone_factor`` and the density in the
    mushy zone is volume-averaged between the solid-side and liquid-side
    table values.

    Parameters
    ----------
    pressure : float
        Pressure in Pa.
    temperature : float
        Temperature in K.
    material_dict : dict
        Material properties dict with 'eos_file' and 'format' keys.
    mushy_zone_factor : float
        Cryoscopic depression factor. 1.0 = no mushy zone (sharp boundary).
        < 1.0 = solidus at this fraction of the extracted liquidus.
    interpolation_functions : dict
        Shared interpolation cache.

    Returns
    -------
    float or None
        Density in kg/m^3, or None on failure.
    """
    eos_file = material_dict['eos_file']
    try:
        cached = _ensure_unified_cache(eos_file, interpolation_functions)

        p_min, p_max = cached['p_min'], cached['p_max']
        if pressure < p_min:
            pressure = p_min
        elif pressure > p_max:
            pressure = p_max

        # W5: math.log10 for scalar input avoids ~1 us/call of numpy-dispatch
        # overhead vs np.log10. Same IEEE-754 result for a Python float input
        # (both ultimately call the same libc log10).
        log_p = math.log10(pressure)
        log_t = math.log10(temperature if temperature > 1.0 else 1.0)

        # Per-cell clamping
        log_t_clamped, was_clamped = _paleos_clamp_temperature(log_p, log_t, cached)
        if was_clamped and eos_file not in _paleos_clamp_warned:
            _paleos_clamp_warned.add(eos_file)
            logger.warning(
                f'PALEOS unified per-cell clamping active for '
                f'{os.path.basename(eos_file)}: '
                f'T={temperature:.0f} K clamped to {10.0**log_t_clamped:.0f} K '
                f'at P={pressure:.2e} Pa.'
            )

        if mushy_zone_factor >= 1.0 or len(cached['liquidus_log_p']) == 0:
            # Direct lookup: no mushy zone (fast bilinear path)
            density = _fast_bilinear(log_p, log_t_clamped, cached['density_grid'], cached)
            if density != density:  # fast NaN check (NaN != NaN)
                density = float(cached['density_nn']((log_p, log_t_clamped)))
            return density if density == density else None

        # Mushy zone: interpolate liquidus T at this P.
        # If query pressure is outside the liquidus coverage (e.g. at
        # pressures where no liquid phase exists), fall back to direct lookup.
        liq_lp = cached['liquidus_log_p']
        if log_p < liq_lp[0] or log_p > liq_lp[-1]:
            density = _fast_bilinear(log_p, log_t_clamped, cached['density_grid'], cached)
            if not np.isfinite(density):
                density = float(cached['density_nn']((log_p, log_t_clamped)))
            return density if np.isfinite(density) else None

        # PALEOS's own melting curve at this pressure
        log_t_melt = float(np.interp(log_p, liq_lp, cached['liquidus_log_t']))
        T_melt = 10.0**log_t_melt

        # Derive mushy zone boundaries (currently from PALEOS liquidus,
        # but may come from external melting curves in future).
        T_liq = T_melt
        T_sol = T_liq * mushy_zone_factor

        # Clamp endpoints against PALEOS's internal phase boundary so that
        # solid-side queries never land on the liquid side and vice versa.
        # The guard offset (_DT_PHASE_GUARD) always shifts T_liq up and may
        # shift T_sol down; only warn when the clamp corrects a genuine
        # cross-boundary incursion (T_sol above T_melt or T_liq below it).
        sol_crossed = T_sol > T_melt
        liq_crossed = T_liq < T_melt
        T_sol = min(T_sol, T_melt - _DT_PHASE_GUARD)
        T_liq = max(T_liq, T_melt + _DT_PHASE_GUARD)

        if sol_crossed or liq_crossed:
            if eos_file not in _paleos_phase_guard_warned:
                _paleos_phase_guard_warned.add(eos_file)
                logger.warning(
                    'Mushy zone endpoints crossed PALEOS melting curve '
                    'for %s at P=%.2e Pa (T_melt=%.1f K): '
                    'T_sol=%.1f K %s, T_liq=%.1f K %s. '
                    'Clamped to safe side.',
                    os.path.basename(eos_file),
                    pressure,
                    T_melt,
                    T_sol,
                    '(was above T_melt)' if sol_crossed else '(ok)',
                    T_liq,
                    '(was below T_melt)' if liq_crossed else '(ok)',
                )
        log_t_sol = math.log10(max(T_sol, 1.0))
        log_t_liq = math.log10(T_liq)

        if temperature >= T_liq:
            # Above liquidus: direct lookup
            density = _fast_bilinear(log_p, log_t_clamped, cached['density_grid'], cached)
            if not np.isfinite(density):
                density = float(cached['density_nn']((log_p, log_t_clamped)))
            return density if np.isfinite(density) else None

        if temperature <= T_sol:
            # Below solidus: direct lookup
            density = _fast_bilinear(log_p, log_t_clamped, cached['density_grid'], cached)
            if not np.isfinite(density):
                density = float(cached['density_nn']((log_p, log_t_clamped)))
            return density if np.isfinite(density) else None

        # In mushy zone: volume-average between solid-side and liquid-side
        phi = (temperature - T_sol) / (T_liq - T_sol)

        # Compute per-P clamp bounds ONCE (inline) rather than calling
        # _paleos_clamp_temperature twice (once for T_sol, once for T_liq) —
        # both clamps share the same log_p, so ip/frac/local_tmin/local_tmax
        # are identical. This saves 2 function calls + dict-lookups per
        # mushy-zone RHS call. Math is identical to _paleos_clamp_temperature;
        # output bit-matches.
        _lt_min = cached['logt_valid_min']
        _lt_max = cached['logt_valid_max']
        if 'dlog_p' in cached:
            _fp = (log_p - cached['logp_min']) / cached['dlog_p']
            _n_p_m1 = cached['n_p'] - 2
            _ip = int(_fp)
            if _ip < 0:  # pragma: no cover - log_p clamped earlier; defensive
                _ip = 0
            elif _ip > _n_p_m1:  # pragma: no cover - log_p clamped earlier; defensive
                _ip = _n_p_m1
            _frac = _fp - _ip
            if (
                _frac < 0.0
            ):  # pragma: no cover - frac is non-negative by construction; defensive
                _frac = 0.0
            elif _frac > 1.0:  # pragma: no cover - frac is bounded by construction; defensive
                _frac = 1.0
            _local_tmin = _lt_min[_ip] + _frac * (_lt_min[_ip + 1] - _lt_min[_ip])
            _local_tmax = _lt_max[_ip] + _frac * (_lt_max[_ip + 1] - _lt_max[_ip])
        else:
            _ulp = cached['unique_log_p']
            _local_tmin = float(np.interp(log_p, _ulp, _lt_min))
            _local_tmax = float(np.interp(log_p, _ulp, _lt_max))

        # NaN-bound fallback matches _paleos_clamp_temperature's contract:
        # if bounds are NaN near table edges, return the input log_t unclamped.
        if _local_tmin == _local_tmin and _local_tmax == _local_tmax:
            if log_t_sol < _local_tmin:
                log_t_sol_c = _local_tmin
            elif log_t_sol > _local_tmax:
                log_t_sol_c = _local_tmax
            else:
                log_t_sol_c = log_t_sol
            if log_t_liq < _local_tmin:
                log_t_liq_c = _local_tmin
            elif log_t_liq > _local_tmax:
                log_t_liq_c = _local_tmax
            else:
                log_t_liq_c = log_t_liq
        else:
            log_t_sol_c = log_t_sol
            log_t_liq_c = log_t_liq

        # Solid-side: density at T_sol
        rho_sol = _fast_bilinear(log_p, log_t_sol_c, cached['density_grid'], cached)
        if not np.isfinite(rho_sol):
            rho_sol = float(cached['density_nn']((log_p, log_t_sol_c)))

        # Liquid-side: density at T_liq
        rho_liq = _fast_bilinear(log_p, log_t_liq_c, cached['density_grid'], cached)
        if not np.isfinite(rho_liq):
            rho_liq = float(cached['density_nn']((log_p, log_t_liq_c)))

        if not (np.isfinite(rho_sol) and np.isfinite(rho_liq)):
            return None

        # Volume additivity
        specific_volume = phi * (1.0 / rho_liq) + (1.0 - phi) * (1.0 / rho_sol)
        return 1.0 / specific_volume

    except Exception as e:
        logger.error(
            f'Error in PALEOS unified density at P={pressure:.2e} Pa, '
            f'T={temperature:.1f} K: {e}'
        )
        return None

def get_paleos_unified_density_batch(
    pressures, temperatures, material_dict, mushy_zone_factor, interpolation_functions
):
    """Vectorized density lookup from a unified PALEOS table.

    Parameters
    ----------
    pressures : numpy.ndarray
        1D array of pressures in Pa.
    temperatures : numpy.ndarray
        1D array of temperatures in K.
    material_dict : dict
        Material properties dict with 'eos_file' and 'format' keys.
    mushy_zone_factor : float
        Cryoscopic depression factor. 1.0 = no mushy zone.
    interpolation_functions : dict
        Shared interpolation cache.

    Returns
    -------
    numpy.ndarray
        1D array of densities in kg/m^3. NaN where lookup fails.
    """
    eos_file = material_dict['eos_file']
    cached = _ensure_unified_cache(eos_file, interpolation_functions)

    p_clamped = np.clip(pressures, cached['p_min'], cached['p_max'])
    log_p = np.log10(p_clamped)
    log_t = np.log10(np.maximum(temperatures, 1.0))

    # Per-cell clamping (vectorized)
    ulp = cached['unique_log_p']
    lt_min = cached['logt_valid_min']
    lt_max = cached['logt_valid_max']
    local_tmin = np.interp(log_p, ulp, lt_min)
    local_tmax = np.interp(log_p, ulp, lt_max)

    valid_bounds = np.isfinite(local_tmin) & np.isfinite(local_tmax)
    log_t_clamped = log_t.copy()
    log_t_clamped = np.where(valid_bounds & (log_t < local_tmin), local_tmin, log_t_clamped)
    log_t_clamped = np.where(valid_bounds & (log_t > local_tmax), local_tmax, log_t_clamped)

    if mushy_zone_factor >= 1.0 or len(cached['liquidus_log_p']) == 0:
        # Direct lookup: fast vectorized bilinear interpolation
        result = fast_bilinear_batch(log_p, log_t_clamped, cached['density_grid'], cached)
        # NN fallback for NaN entries
        nan_mask = ~np.isfinite(result)
        if np.any(nan_mask):
            pts_nn = np.column_stack([log_p[nan_mask], log_t_clamped[nan_mask]])
            result[nan_mask] = cached['density_nn'](pts_nn)
        return result

    # Mushy zone path (vectorized).
    # Compute T_melt from the PALEOS analytic melting curve rather than
    # interpolating the extracted liquidus grid. This is faster and avoids
    # branching per-element for the "outside liquidus coverage" case.
    from ..melting_curves import paleos_liquidus

    T_melt = paleos_liquidus(pressures)

    # Derive mushy zone boundaries (currently from PALEOS liquidus,
    # but may come from external melting curves in future).
    T_liq = T_melt.copy()
    T_sol = T_liq * mushy_zone_factor

    # Clamp endpoints against PALEOS's own melting curve
    T_sol = np.minimum(T_sol, T_melt - _DT_PHASE_GUARD)
    T_liq = np.maximum(T_liq, T_melt + _DT_PHASE_GUARD)

    log_t_sol = np.log10(np.maximum(T_sol, 1.0))
    log_t_liq = np.log10(np.maximum(T_liq, 1.0))

    # Classify shells: above liquidus, below solidus, or in mushy zone
    above = temperatures >= T_liq
    below = temperatures <= T_sol
    mushy = ~above & ~below

    # Direct lookup for above-liquidus and below-solidus shells
    result = fast_bilinear_batch(log_p, log_t_clamped, cached['density_grid'], cached)
    nan_mask = ~np.isfinite(result)
    if np.any(nan_mask):
        pts_nn = np.column_stack([log_p[nan_mask], log_t_clamped[nan_mask]])
        result[nan_mask] = cached['density_nn'](pts_nn)

    # Mushy zone shells: volume-average between solid-side and liquid-side
    if np.any(mushy):
        m_idx = np.where(mushy)[0]
        phi = (temperatures[m_idx] - T_sol[m_idx]) / (T_liq[m_idx] - T_sol[m_idx])

        # Solid-side density at T_sol
        log_t_sol_c = log_t_sol[m_idx].copy()
        sol_tmin = np.interp(log_p[m_idx], ulp, lt_min)
        sol_tmax = np.interp(log_p[m_idx], ulp, lt_max)
        sol_valid = np.isfinite(sol_tmin) & np.isfinite(sol_tmax)
        log_t_sol_c = np.where(sol_valid & (log_t_sol_c < sol_tmin), sol_tmin, log_t_sol_c)
        log_t_sol_c = np.where(sol_valid & (log_t_sol_c > sol_tmax), sol_tmax, log_t_sol_c)
        rho_sol = fast_bilinear_batch(log_p[m_idx], log_t_sol_c, cached['density_grid'], cached)
        nn_sol = ~np.isfinite(rho_sol)
        if np.any(nn_sol):
            pts_sol_nn = np.column_stack([log_p[m_idx][nn_sol], log_t_sol_c[nn_sol]])
            rho_sol[nn_sol] = cached['density_nn'](pts_sol_nn)

        # Liquid-side density at T_liq
        log_t_liq_c = log_t_liq[m_idx].copy()
        liq_tmin = np.interp(log_p[m_idx], ulp, lt_min)
        liq_tmax = np.interp(log_p[m_idx], ulp, lt_max)
        liq_valid = np.isfinite(liq_tmin) & np.isfinite(liq_tmax)
        log_t_liq_c = np.where(liq_valid & (log_t_liq_c < liq_tmin), liq_tmin, log_t_liq_c)
        log_t_liq_c = np.where(liq_valid & (log_t_liq_c > liq_tmax), liq_tmax, log_t_liq_c)
        rho_liq = fast_bilinear_batch(log_p[m_idx], log_t_liq_c, cached['density_grid'], cached)
        nn_liq = ~np.isfinite(rho_liq)
        if np.any(nn_liq):
            pts_liq_nn = np.column_stack([log_p[m_idx][nn_liq], log_t_liq_c[nn_liq]])
            rho_liq[nn_liq] = cached['density_nn'](pts_liq_nn)

        # Volume additivity
        both_ok = np.isfinite(rho_sol) & np.isfinite(rho_liq) & (rho_sol > 0) & (rho_liq > 0)
        spec_vol = phi * (1.0 / np.where(both_ok, rho_liq, 1.0)) + (1.0 - phi) * (
            1.0 / np.where(both_ok, rho_sol, 1.0)
        )
        result[m_idx] = np.where(both_ok, 1.0 / spec_vol, np.nan)

    return result