Zalmoxis.eos functions

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:
            return None
        # Fall back to whichever is available if one is None
        if nabla_solid is None:
            nabla = nabla_liquid
        elif nabla_liquid is None:
            nabla = nabla_solid
        else:
            nabla = (1.0 - phi) * nabla_solid + phi * nabla_liquid

    if nabla is None or nabla <= 0:
        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.

    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:
        interpolation_functions[eos_file] = load_paleos_unified_table(eos_file)
    return interpolation_functions[eos_file]

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

    Replaces scipy's RegularGridInterpolator for scalar lookups. The
    index computation uses O(1) arithmetic with a rounding correction
    to handle floating-point spacing jitter in the grid axes, then
    computes the fractional position from the actual stored grid
    coordinates for exact interpolation.

    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 (density or nabla_ad), shape (n_p, n_t).
    cached : dict
        Cache entry with grid metadata (logp_min, dlog_p, unique_log_p,
        unique_log_t, etc.).

    Returns
    -------
    float
        Interpolated value. NaN if any corner is NaN.
    """
    # Guard: grids with fewer than 2 points cannot be interpolated
    if cached['n_p'] < 2 or cached['n_t'] < 2:
        return grid[0, 0]

    ulp = cached['unique_log_p']
    ult = cached['unique_log_t']
    n_p = cached['n_p']
    n_t = cached['n_t']

    # O(1) index estimation for nominally log-uniform grid.
    # Use round() to get the nearest node, then determine the lower
    # bounding index. This handles floating-point spacing jitter that
    # causes int() to land on the wrong cell.
    fp = (log_p - cached['logp_min']) / cached['dlog_p']
    ft = (log_t - cached['logt_min']) / cached['dlog_t']

    # Nearest node via rounding, then find lower bounding index
    ip_near = min(max(int(fp + 0.5), 0), n_p - 1)
    it_near = min(max(int(ft + 0.5), 0), n_t - 1)

    # Determine lower bounding index: if the query is below the
    # nearest node, step back by one
    if ip_near > 0 and log_p < ulp[ip_near]:
        ip = ip_near - 1
    else:
        ip = ip_near
    ip = max(0, min(ip, n_p - 2))

    if it_near > 0 and log_t < ult[it_near]:
        it = it_near - 1
    else:
        it = it_near
    it = max(0, min(it, n_t - 2))

    # Fractional parts computed from actual grid coordinates
    span_p = ulp[ip + 1] - ulp[ip]
    span_t = ult[it + 1] - ult[it]

    dp = (log_p - ulp[ip]) / span_p if span_p > 0 else 0.0
    dt = (log_t - ult[it]) / span_t if span_t > 0 else 0.0

    # Clamp to [0, 1] for boundary safety
    dp = max(0.0, min(dp, 1.0))
    dt = max(0.0, min(dt, 1.0))

    # Four corner values
    v00 = grid[ip, it]
    v01 = grid[ip, it + 1]
    v10 = grid[ip + 1, it]
    v11 = grid[ip + 1, it + 1]

    # Bilinear interpolation
    return v00 * (1 - dp) * (1 - dt) + v01 * (1 - dp) * dt + v10 * dp * (1 - dt) + v11 * 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.

    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 from ``load_paleos_table()``.

    Returns
    -------
    float
        Clamped log10(T), guaranteed to fall within the valid data region
        at the given pressure.
    bool
        True if clamping was applied, False if the original value was valid.
    """
    ulp = cached['unique_log_p']
    lt_min = cached['logt_valid_min']
    lt_max = cached['logt_valid_max']

    # Interpolate per-pressure T bounds at the query pressure
    local_tmin = float(np.interp(log_p, ulp, lt_min))
    local_tmax = float(np.interp(log_p, ulp, lt_max))

    # Guard: if the pressure is near an all-NaN row, the interpolated
    # bounds are NaN. Return unclamped and let the NN fallback handle it.
    if not (np.isfinite(local_tmin) and np.isfinite(local_tmax)):
        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)

    # 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}'.")

    # 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)
    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,
):
    """
    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"}
        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.
    surface_temperature : float
        Temperature at the surface, in K.
    center_temperature : float
        Temperature at the center, in K. Used for ``"linear"`` and
        ``"adiabatic"`` modes.
    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])
        )

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

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,
):
    """
    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.
    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.

    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)
    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
    """

    pressure_file = os.path.join(ZALMOXIS_ROOT, 'output_files', 'pressure_profiles.txt')
    density_file = os.path.join(ZALMOXIS_ROOT, 'output_files', '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,
            )
        frac_melt = (temperature - T_sol) / (T_liq - T_sol)
        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_analytic_density(pressure: float, material_key: str) -> float | None:
    """
    Compute density from the Seager et al. (2007) analytic modified polytrope.

    Parameters
    ----------
    pressure : float
        Pressure in Pa.
    material_key : str
        One of the keys in SEAGER2007_MATERIALS
        (e.g., "iron", "MgSiO3", "H2O", "graphite", "SiC", "MgFeSiO3").

    Returns
    -------
    float
        Density in kg/m^3. Returns rho_0 for non-positive pressure (allows
        the ODE solver to remain stable during pressure adjustment iterations).
        Returns None only for NaN pressure.

    Raises
    ------
    ValueError
        If material_key is not recognized.
    """
    if material_key not in SEAGER2007_MATERIALS:
        raise ValueError(
            f"Unknown material key '{material_key}'. Valid keys: {sorted(VALID_MATERIAL_KEYS)}"
        )

    rho_0, c, n = SEAGER2007_MATERIALS[material_key]

    # Guard against nonphysical pressure
    if np.isnan(pressure):
        return None
    if pressure <= 0:
        return float(rho_0)

    if pressure > P_MAX:
        logger.warning(
            f'Pressure {pressure:.2e} Pa exceeds validity limit of {P_MAX:.0e} Pa '
            f'for Seager+2007 analytic EOS. Results may be inaccurate.'
        )

    density = rho_0 + c * pressure**n

    return float(density)

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 or pressure > p_max:
            logger.debug(
                f'PALEOS unified: P={pressure:.2e} Pa out of bounds '
                f'[{p_min:.2e}, {p_max:.2e}]. Clamping.'
            )
            pressure = np.clip(pressure, p_min, p_max)

        log_p = np.log10(pressure)
        log_t = np.log10(max(temperature, 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 not np.isfinite(density):
                density = float(cached['density_nn']((log_p, log_t_clamped)))
            return density if np.isfinite(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 = np.log10(max(T_sol, 1.0))
        log_t_liq = np.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)

        # Solid-side: density at T_sol
        log_t_sol_c, _ = _paleos_clamp_temperature(log_p, log_t_sol, cached)
        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
        log_t_liq_c, _ = _paleos_clamp_temperature(log_p, log_t_liq, cached)
        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: batch interpolation
        pts = np.column_stack([log_p, log_t_clamped])
        result = cached['density_interp'](pts)
        # NN fallback for NaN entries
        nan_mask = ~np.isfinite(result)
        if np.any(nan_mask):
            result[nan_mask] = cached['density_nn'](pts[nan_mask])
        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
    pts = np.column_stack([log_p, log_t_clamped])
    result = cached['density_interp'](pts)
    nan_mask = ~np.isfinite(result)
    if np.any(nan_mask):
        result[nan_mask] = cached['density_nn'](pts[nan_mask])

    # 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)
        pts_sol = np.column_stack([log_p[m_idx], log_t_sol_c])
        rho_sol = cached['density_interp'](pts_sol)
        nn_sol = ~np.isfinite(rho_sol)
        if np.any(nn_sol):
            rho_sol[nn_sol] = cached['density_nn'](pts_sol[nn_sol])

        # 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)
        pts_liq = np.column_stack([log_p[m_idx], log_t_liq_c])
        rho_liq = cached['density_interp'](pts_liq)
        nn_liq = ~np.isfinite(rho_liq)
        if np.any(nn_liq):
            rho_liq[nn_liq] = cached['density_nn'](pts_liq[nn_liq])

        # 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