aragog.eos

def dTdPs(self, P: npt.NDArray | float, S: npt.NDArray | float) -> npt.NDArray:
    """Adiabatic temperature gradient dT/dP|_S (P, S) [K/Pa].

    This is the SPIDER convention: dT/dP along the adiabat, NOT
    the dimensionless nabla_ad = d ln T / d ln P.
    """
    return self._lookup_phase_weighted('dTdPs', P, S)

def density(self, P: npt.NDArray | float, S: npt.NDArray | float) -> npt.NDArray:
    """Density rho(P, S) [kg/m^3].

    In the mushy zone (0 < phi < 1): harmonic mean of end-member
    densities evaluated at phase boundary entropies (SPIDER
    eos_composite.c:236-237).

    Outside the mushy zone (phi=0 or phi=1): single-phase table
    evaluated at the actual S (clamped to table range), matching
    SPIDER's combine_matprop(smth=0, mixed, single) path.
    """
    P = np.asarray(P, dtype=float)
    S = np.asarray(S, dtype=float)
    phi = self.melt_fraction(P, S)
    mushy = (phi > 0) & (phi < 1)

    solid_table = self._tables['density_solid']
    melt_table = self._tables['density_melt']

    # Mushy zone: evaluate at phase boundaries (Lever Rule)
    rho_sol_boundary = self._lookup_at_phase_boundary('density', P, 'solid')
    rho_liq_boundary = self._lookup_at_phase_boundary('density', P, 'melt')
    inv_rho_mushy = phi / np.maximum(rho_liq_boundary, 1.0) + (1.0 - phi) / np.maximum(
        rho_sol_boundary, 1.0
    )
    rho_mushy = 1.0 / np.maximum(inv_rho_mushy, 1e-30)

    # Single-phase: evaluate at actual S (clamped to table range)
    S_solid_clamped = np.clip(S, solid_table['S'][0], solid_table['S'][-1])
    S_melt_clamped = np.clip(S, melt_table['S'][0], melt_table['S'][-1])
    P_solid_clamped = np.clip(P, solid_table['P'][0], solid_table['P'][-1])
    P_melt_clamped = np.clip(P, melt_table['P'][0], melt_table['P'][-1])
    pts_solid = np.column_stack([P_solid_clamped.ravel(), S_solid_clamped.ravel()])
    pts_melt = np.column_stack([P_melt_clamped.ravel(), S_melt_clamped.ravel()])
    rho_solid_single = solid_table['interp'](pts_solid).reshape(P.shape)
    rho_melt_single = melt_table['interp'](pts_melt).reshape(P.shape)
    # 0.5 here is a binary table-selector for the non-mushy fallback
    # (phi is essentially 0 or 1 outside the mushy band); it is NOT
    # the rheological critical melt fraction. The RCMF is
    # ``EntropyPhaseEvaluator._phi_rheo`` and drives the viscosity
    # tanh blend separately.
    rho_single = np.where(phi >= 0.5, rho_melt_single, rho_solid_single)

    return np.where(mushy, rho_mushy, rho_single)

def heat_capacity(self, P: npt.NDArray | float, S: npt.NDArray | float) -> npt.NDArray:
    """Specific heat capacity Cp(P, S) [J/kg/K].

    Linear blend of solid and melt P-S tables by melt fraction
    (the ``cp_blend = 'linear'`` convention). Use
    ``heat_capacity_latent_blend`` for SPIDER-parity mushy-zone
    Cp that includes the latent-heat contribution.
    """
    return self._lookup_phase_weighted('heat_capacity', P, S)

def heat_capacity_latent_blend(
    self, P: npt.NDArray | float, S: npt.NDArray | float, width: float = 0.01
) -> npt.NDArray:
    """Latent-heat-augmented Cp(P, S), matching SPIDER convention.

    Mirrors SPIDER's ``eos_composite.c:226-232`` formula:

        Cp_mix = (S_liq - S_sol) / (T_liq - T_sol)
               * (T_sol + 0.5 * (T_liq - T_sol))

    which captures the latent-heat contribution to apparent heat
    capacity in the mushy zone. Outside the phase boundary the
    result reduces to the pure-phase Cp via tanh smoothing.

    The plain linear blend used by ``heat_capacity`` above misses
    the latent-heat term and underestimates Cp by up to a factor
    of 6 in the deep mushy regime relative to SPIDER's internal
    cp_s. This formula closes that gap.

    Parameters
    ----------
    P : array or float
        Pressure [Pa].
    S : array or float
        Entropy [J/kg/K].
    width : float
        Tanh smoothing width across the phase boundary, in units
        of melt fraction. SPIDER default is 0.01.

    Returns
    -------
    ndarray
        Heat capacity [J/kg/K] with the same shape as the inputs.
    """
    P = np.asarray(P, dtype=float)
    S = np.asarray(S, dtype=float)

    # Phase boundaries at this pressure
    S_sol = np.asarray(self.solidus_entropy(P)).reshape(P.shape)
    S_liq = np.asarray(self.liquidus_entropy(P)).reshape(P.shape)
    T_sol = np.asarray(self._lookup_at_phase_boundary('temperature', P, 'solid')).reshape(
        P.shape
    )
    T_liq = np.asarray(self._lookup_at_phase_boundary('temperature', P, 'melt')).reshape(
        P.shape
    )

    dS = S_liq - S_sol
    dT = T_liq - T_sol

    # Mixed-phase Cp from the latent-heat formula
    # SPIDER: Cp_mix = (S_liq - S_sol) / (T_liq - T_sol) * T_mid
    safe_dT = np.where(np.abs(dT) > 1e-10, dT, 1e-10)
    T_mid = T_sol + 0.5 * dT
    Cp_mix = (dS / safe_dT) * T_mid

    # Pure-phase Cp from the linear-blend tables
    Cp_pure = self._lookup_phase_weighted('heat_capacity', P, S)

    # Un-truncated phase fraction (SPIDER's gphi convention)
    safe_dS = np.where(np.abs(dS) > 1e-10, dS, 1e-10)
    gphi = (S - S_sol) / safe_dS

    # SPIDER-parity smoothing (util.c:245-270, get_smoothing).
    # Bell shape centred on the mushy zone: rises at gphi=0 with
    # tanh(gphi/w), falls at gphi=1 with tanh((1-gphi)/w). The
    # width parameter w controls the shoulder into the single-
    # phase regime.
    #
    # IMPORTANT: SPIDER uses tanh(x/w), not expit(x/w) which
    # equals 0.5*(1+tanh(x/(2w))). The previous expit formula
    # had an effective width of w/2, making the transition twice
    # as narrow and producing 16% less latent-heat Cp at
    # crystallization onset (gphi=0.994).
    w = max(width, 1e-6)
    smooth = np.where(
        gphi > 0.5,
        0.5 * (1.0 + np.tanh((1.0 - gphi) / w)),  # fall at liquidus
        0.5 * (1.0 + np.tanh(gphi / w)),  # rise at solidus
    )

    return smooth * Cp_mix + (1.0 - smooth) * Cp_pure

def invert_temperature(self, P: float, T_target: float) -> float:
    """Find entropy S such that T(P, S) = T_target.

    Uses Brent root-finding with pure-Python bilinear interpolation
    on the P-S temperature tables. All arithmetic inside the brentq
    callback uses only Python floats, avoiding numpy scalar
    conversions that break on numpy >= 2.4 (where ndim > 0 arrays
    cannot be implicitly converted to scalars).

    Parameters
    ----------
    P : float
        Pressure [Pa].
    T_target : float
        Target temperature [K].

    Returns
    -------
    float
        Entropy [J/kg/K] such that T(P, S) ~ T_target.

    Raises
    ------
    ValueError
        If T_target is outside the range of T(P, S) for this P.
    """
    from scipy.optimize import brentq

    # All values must be pure Python floats for numpy 2.4 safety.
    P_clamped = max(self.P_min, min(float(P), self.P_max))
    T_tgt = float(T_target)

    def residual(S_cand: float) -> float:
        # Pure-Python path: no numpy arrays created or converted.
        return self.temperature_scalar(P_clamped, S_cand) - T_tgt

    S_lo = self.S_min  # already Python float from __init__
    S_hi = self.S_max
    f_lo = residual(S_lo)
    f_hi = residual(S_hi)

    if f_lo * f_hi > 0:
        raise ValueError(
            f'Cannot invert T={T_target:.1f} K at P={P:.2e} Pa: '
            f'T(S_min={S_lo:.0f})={T_tgt + f_lo:.0f} K, '
            f'T(S_max={S_hi:.0f})={T_tgt + f_hi:.0f} K. '
            f'Target outside table range.'
        )

    S_root = brentq(residual, S_lo, S_hi, xtol=0.1, rtol=1e-10)
    # brentq returns a Python float (from C double), but be safe
    return float(S_root)

def latent_heat(self, P: npt.NDArray | float) -> npt.NDArray:
    """Latent heat L(P) = T_fus × (S_liq - S_sol) [J/kg].

    P-dependent, following SPIDER convention. T_fus is the average
    of solidus and liquidus temperatures at the given pressure.
    """
    P = np.asarray(P, dtype=float)
    S_sol = self.solidus_entropy(P)
    S_liq = self.liquidus_entropy(P)
    T_sol = self._lookup_at_phase_boundary('temperature', P, 'solid')
    T_liq = self._lookup_at_phase_boundary('temperature', P, 'melt')
    T_fus = 0.5 * (T_sol + T_liq)
    return T_fus * np.maximum(S_liq - S_sol, 1.0)

def liquidus_entropy(self, P: npt.NDArray | float) -> npt.NDArray:
    """Liquidus entropy S_liq(P) [J/kg/K]."""
    return np.asarray(self._liquidus['interp'](P), dtype=float)

def liquidus_entropy_dP(self, P: npt.NDArray | float) -> npt.NDArray:
    """dS_liq/dP at the given pressure(s), in J/(kg·K·Pa)."""
    return np.asarray(self._liquidus['dinterp'](P), dtype=float)

def melt_fraction(self, P: npt.NDArray | float, S: npt.NDArray | float) -> npt.NDArray:
    """Melt fraction phi from position between solidus and liquidus.

    phi = 0 for S <= S_sol, phi = 1 for S >= S_liq, linear between.
    """
    P = np.asarray(P, dtype=float)
    S = np.asarray(S, dtype=float)
    S_sol = self.solidus_entropy(P)
    S_liq = self.liquidus_entropy(P)
    dS = np.maximum(S_liq - S_sol, 1e-10)
    phi = np.clip((S - S_sol) / dS, 0.0, 1.0)
    return phi

def solidus_entropy(self, P: npt.NDArray | float) -> npt.NDArray:
    """Solidus entropy S_sol(P) [J/kg/K]."""
    return np.asarray(self._solidus['interp'](P), dtype=float)

def solidus_entropy_dP(self, P: npt.NDArray | float) -> npt.NDArray:
    """dS_sol/dP at the given pressure(s), in J/(kg·K·Pa).

    Needed by the SPIDER-parity mixing-flux formula in
    ``entropy_state.update``. SPIDER computes Jmix via a bracket
    expression ``dS/dr − [φ dS_liq/dP + (1−φ) dS_sol/dP] dP/dr``
    (``energy.c::GetMixingHeatFlux`` lines 307-309). Exposing the
    phase-boundary P-derivatives here lets Aragog match SPIDER's
    formula exactly without resorting to un-truncated gphi
    arithmetic.
    """
    return np.asarray(self._solidus['dinterp'](P), dtype=float)

def specific_enthalpy(
    self,
    P: npt.NDArray | float,
    S: npt.NDArray | float,
) -> npt.NDArray:
    """EOS-consistent specific enthalpy h(P, S) [J/kg].

    Bilinear interpolation on the precomputed table built in
    ``_build_enthalpy_table``. Inputs outside the table P or S
    range are clamped to the nearest table edge so the lookup
    never returns NaN; callers querying out-of-range states
    receive the boundary value instead of an exception.

    Parameters
    ----------
    P : array or float
        Pressure [Pa].
    S : array or float
        Entropy [J/kg/K].

    Returns
    -------
    ndarray
        Specific enthalpy [J/kg], same shape as broadcast(P, S).
    """
    P = np.asarray(P, dtype=float)
    S = np.asarray(S, dtype=float)
    P_clamped = np.clip(P, self._h_grid_P[0], self._h_grid_P[-1])
    S_clamped = np.clip(S, self._h_grid_S[0], self._h_grid_S[-1])
    pts = np.column_stack([P_clamped.ravel(), S_clamped.ravel()])
    return self._h_interp(pts).reshape(P_clamped.shape)

def temperature(self, P: npt.NDArray | float, S: npt.NDArray | float) -> npt.NDArray:
    """Temperature T(P, S) [K]."""
    return self._lookup_phase_weighted('temperature', P, S)

def temperature_scalar(self, P: float, S: float) -> float:
    """Temperature T(P, S) for a single point, pure Python.

    Uses the same Lever Rule blending as ``temperature()`` but
    with pure-Python arithmetic throughout. Safe for use inside
    scipy.optimize callbacks on numpy >= 2.4.

    Parameters
    ----------
    P : float
        Pressure [Pa].
    S : float
        Entropy [J/kg/K].

    Returns
    -------
    float
        Temperature [K].
    """
    phi = self._melt_fraction_scalar(P, S)

    solid_table = self._tables['temperature_solid']
    melt_table = self._tables['temperature_melt']

    # Lever Rule: in the mushy zone, evaluate at phase boundary
    # entropies; outside, evaluate at actual S.
    if 0.0 < phi < 1.0:
        S_for_solid = self._solidus_entropy_scalar(P)
        S_for_melt = self._liquidus_entropy_scalar(P)
    else:
        S_for_solid = S
        S_for_melt = S

    # Clamp to each table's own range
    S_sol_clamped = max(
        solid_table['S_list'][0], min(S_for_solid, solid_table['S_list'][-1])
    )
    S_melt_clamped = max(melt_table['S_list'][0], min(S_for_melt, melt_table['S_list'][-1]))
    P_sol_clamped = max(solid_table['P_list'][0], min(P, solid_table['P_list'][-1]))
    P_melt_clamped = max(melt_table['P_list'][0], min(P, melt_table['P_list'][-1]))

    val_solid = _interp_bilinear_scalar(
        solid_table['P_list'],
        solid_table['S_list'],
        solid_table['values'],
        P_sol_clamped,
        S_sol_clamped,
    )
    val_melt = _interp_bilinear_scalar(
        melt_table['P_list'],
        melt_table['S_list'],
        melt_table['values'],
        P_melt_clamped,
        S_melt_clamped,
    )

    # NaN-safe phase-weighted blend
    result = 0.0
    if phi > 0.0:
        result += phi * val_melt
    if phi < 1.0:
        result += (1.0 - phi) * val_solid
    return result

def thermal_expansivity(
    self, P: npt.NDArray | float, S: npt.NDArray | float
) -> npt.NDArray:
    """Thermal expansivity alpha(P, S) [1/K].

    When P-S thermal_exp tables are available (generated alongside
    SPIDER tables), uses them directly for exact parity with SPIDER.
    Otherwise, derives alpha from the thermodynamic identity:
    alpha = rho * Cp * |dT/dP|_S| / T.
    """
    if self._has_alpha_tables:
        return self._lookup_phase_weighted('thermal_exp', P, S)

    T = self.temperature(P, S)
    rho = self.density(P, S)
    Cp = self.heat_capacity(P, S)
    dTdPs_val = self.dTdPs(P, S)
    alpha = rho * Cp * np.abs(dTdPs_val) / np.maximum(T, 1.0)
    return alpha

def thermal_expansivity_composite_blend(
    self,
    P: npt.NDArray | float,
    S: npt.NDArray | float,
    width: float = 0.01,
) -> npt.NDArray:
    """Composite two-phase thermal expansivity, matching SPIDER.

    Mirrors SPIDER's ``eos_composite.c:246`` formula:

        alpha_mix = (rho_solid - rho_melt) / (T_liq - T_sol) / rho

    where ``rho`` is the harmonic-mean composite density. This
    effective alpha captures the density change from the phase
    transition (Clausius-Clapeyron contribution), which dominates
    over the single-phase alpha by a factor of ~5 in the deep
    mushy zone. Outside the mushy band the result reduces to the
    pure-phase alpha from the P-S tables via tanh smoothing
    (SPIDER ``eos_composite.c:279``).

    Parameters
    ----------
    P : array or float
        Pressure [Pa].
    S : array or float
        Entropy [J/kg/K].
    width : float
        Tanh smoothing width across the phase boundary.

    Returns
    -------
    ndarray
        Thermal expansivity [1/K].
    """
    P = np.asarray(P, dtype=float)
    S = np.asarray(S, dtype=float)

    # Phase-boundary quantities at this pressure
    S_sol = np.asarray(self.solidus_entropy(P)).reshape(P.shape)
    S_liq = np.asarray(self.liquidus_entropy(P)).reshape(P.shape)
    T_sol = np.asarray(self._lookup_at_phase_boundary('temperature', P, 'solid')).reshape(
        P.shape
    )
    T_liq = np.asarray(self._lookup_at_phase_boundary('temperature', P, 'melt')).reshape(
        P.shape
    )
    rho_sol = np.asarray(self._lookup_at_phase_boundary('density', P, 'solid')).reshape(
        P.shape
    )
    rho_liq = np.asarray(self._lookup_at_phase_boundary('density', P, 'melt')).reshape(
        P.shape
    )

    # Composite density (harmonic mean, SPIDER eos_composite.c:236)
    phi = self.melt_fraction(P, S)
    inv_rho = phi / np.maximum(rho_liq, 1.0) + (1.0 - phi) / np.maximum(rho_sol, 1.0)
    rho_comp = 1.0 / np.maximum(inv_rho, 1e-30)

    # Mixed-phase alpha (SPIDER eos_composite.c:246)
    dT = T_liq - T_sol
    safe_dT = np.where(np.abs(dT) > 1e-10, dT, 1e-10)
    alpha_mix = (rho_sol - rho_liq) / safe_dT / rho_comp

    # Pure-phase alpha from the linear-blend tables
    alpha_pure = self._lookup_phase_weighted('thermal_exp', P, S)

    # SPIDER-parity tanh smoothing (same as heat_capacity_latent_blend)
    dS = S_liq - S_sol
    safe_dS = np.where(np.abs(dS) > 1e-10, dS, 1e-10)
    gphi = (S - S_sol) / safe_dS

    w = max(width, 1e-6)
    smooth = np.where(
        gphi > 0.5,
        0.5 * (1.0 + np.tanh((1.0 - gphi) / w)),
        0.5 * (1.0 + np.tanh(gphi / w)),
    )

    return smooth * alpha_mix + (1.0 - smooth) * alpha_pure

def capacitance(self) -> npt.NDArray:
    """Capacitance for the entropy equation: rho * T [kg K / m^3].

    The entropy equation is rho * T * dS/dt = div(F).
    """
    return self._density * self._temperature

def dTdPs(self) -> npt.NDArray:
    """Adiabatic temperature gradient dT/dP|_S [K/Pa]."""
    return self._dTdPs_val

def dTdrs(self) -> npt.NDArray:
    """Adiabatic temperature gradient dT/dr|_S [K/m].

    dT/dr|_S = -g * alpha * T / Cp
    """
    return -self._g * self._thermal_expansivity * self._temperature / self._heat_capacity

def delta_specific_volume(self) -> FloatOrArray:
    """Specific volume difference between solid and liquid [m^3/kg]."""
    if self._const_properties:
        return np.zeros_like(self._density)
    rho_s = self._eos._lookup_at_phase_boundary('density', self.pressure, 'solid')
    rho_l = self._eos._lookup_at_phase_boundary('density', self.pressure, 'melt')
    return 1.0 / np.maximum(rho_l, 1.0) - 1.0 / np.maximum(rho_s, 1.0)

def latent_heat(self) -> FloatOrArray:
    """Latent heat L(P) = T_fus × (S_liq - S_sol) [J/kg].

    P-dependent from EOS tables, matching SPIDER convention.
    """
    return self._latent_heat_val

def relative_velocity(self) -> FloatOrArray:
    """Melt-solid relative velocity for gravitational separation [m/s].

    Returns zero in const_properties mode (no phase contrast).

    Uses Abe (1993) three-regime permeability model based on porosity
    (volume fraction of melt, not mass fraction), matching SPIDER's
    GetGravitationalHeatFlux in energy.c.

    Regimes (F = K/porosity, the quantity multiplying delta_rho*g/eta):
    1. Blake-Kozeny-Carman (low porosity): F = d^2 por^2 / ((1-por)^2 * 1000)
    2. Rumpf-Gupte (intermediate): F = d^2 por^4.5 * (5/7)
    3. Stokes settling (high porosity): F = d^2 * 2/9
    """
    if self._const_properties:
        return np.zeros_like(self._density)
    rho_s = self._eos._lookup_at_phase_boundary('density', self.pressure, 'solid')
    rho_l = self._eos._lookup_at_phase_boundary('density', self.pressure, 'melt')
    delta_rho = rho_l - rho_s  # typically negative (melt lighter)
    g = self._g
    d = self._grain_size
    eta_l = self._visc_liquid

    # Porosity (volume fraction of melt) from densities. Smoothed
    # with sqrt-based soft clip + soft max so the CVODE BDF
    # predictor sees a C^infty RHS (previous np.clip + np.maximum
    # had two derivative jumps that locked the solver at order 1).
    rho = self._density
    drho = rho_s - rho_l
    eps = 1.0e-3  # kg/m^3; far below any physical density contrast
    drho_smoothmax = 0.5 * (drho + 1.0 + np.sqrt((drho - 1.0) ** 2 + eps * eps))
    porosity_raw = (rho_s - rho) / drho_smoothmax
    # smooth_clip(porosity_raw, 0, 1) via two soft-max operations
    eps_p = 1.0e-3  # dimensionless; invisible except near [0,1] edges
    p_lo = 0.5 * (porosity_raw + np.sqrt(porosity_raw * porosity_raw + eps_p * eps_p))
    hi_u = 1.0 - p_lo
    porosity = 1.0 - 0.5 * (hi_u + np.sqrt(hi_u * hi_u + eps_p * eps_p))

    # Three-regime permeability / porosity (Abe 1993/1995, SPIDER convention).
    # F = permeability(porosity) / porosity. The relative velocity is
    # v = |delta_rho| * g * F / eta_liquid.
    por = np.maximum(porosity, 1e-20)
    one_m_por = np.maximum(1.0 - porosity, 1e-20)

    # Blake-Kozeny-Carman (low porosity): K/por = d^2 por^2 / ((1-por)^2 * 1000)
    F_bkc = d**2 * por**2 / (one_m_por**2 * 1000.0)
    # Rumpf-Gupte (intermediate): K/por = d^2 por^4.5 * 5/7
    F_rg = d**2 * por**4.5 * (5.0 / 7.0)
    # Stokes settling (high porosity): K/por = d^2 * 2/9
    F_stokes = d**2 * 2.0 / 9.0

    # Regime switching at critical porosities. These two values
    # (zeta_1 = 0.0769452, zeta_2 = 0.771462) are NOT material
    # parameters: they are the analytical equality points of the
    # three permeability laws (BKC = RG at zeta_1, RG = Stokes at
    # zeta_2) under the equal-density-ratio limit, derived in
    # Bower et al. 2018 Eqs. 13a-c. They depend only on the
    # functional forms of the three regime laws, not on the
    # underlying material; varying composition changes individual
    # F_bkc / F_rg / F_stokes prefactors but the equality points
    # remain fixed by construction. The tanh blend widths (0.02
    # and 0.05) are numerical-smoothing tunables, not physics.
    w_rg = tanh_weight(porosity, 0.0769452, 0.02)
    w_stokes = tanh_weight(porosity, 0.771462, 0.05)
    F = (1.0 - w_rg) * F_bkc + (w_rg - w_stokes) * F_rg + w_stokes * F_stokes
    F = np.maximum(F, 0.0)

    # Relative velocity: v = |delta_rho| * g * F / eta_liquid
    # Sign convention: positive = outward (melt rising, solid sinking)
    # Smoothed |delta_rho| via sqrt(x^2 + eps^2); eps tiny compared
    # to physical delta_rho ~ -500 kg/m^3, so bulk result unchanged.
    abs_drho = np.sqrt(delta_rho * delta_rho + 1.0e-12)
    v_rel = abs_drho * g * F / np.maximum(eta_l, 1e-10)

    return v_rel

def set_entropy(self, entropy: npt.NDArray) -> None:
    """Set the entropy profile."""
    self.entropy = np.asarray(entropy, dtype=float)

def set_pressure(self, pressure: npt.NDArray) -> None:
    """Set the pressure profile."""
    self.pressure = np.asarray(pressure, dtype=float)

def set_temperature(self, temperature: npt.NDArray) -> None:
    """Not used in entropy mode. Use set_entropy() instead."""
    raise NotImplementedError(
        'EntropyPhaseEvaluator uses set_entropy(), not set_temperature(). '
        'This evaluator cannot be used with the T-based solver.'
    )

def temperature(self) -> FloatOrArray:
    """Temperature from EOS lookup (not a state variable)."""
    return self._temperature

def update(self) -> None:
    """Recompute all cached properties from current (P, S).

    Two modes:
    - const_properties=True: analytical T(S) with constant rho, Cp,
      alpha, k, visc. Matches SPIDER's -use_const_properties. No EOS
      table lookups. phi=1 always (no phase transitions).
    - const_properties=False (default): single-pass EOS evaluation
      mirroring SPIDER's EOSEval_Composite_TwoPhase.
    """
    if self._const_properties:
        return self._update_const()
    return self._update_eos()