Zalmoxis.binodal

def mass_to_mole_fraction(w_1, w_2, mu_1, mu_2):
    """Convert mass fractions to mole fraction of component 1.

    Parameters
    ----------
    w_1 : float
        Mass fraction of component 1.
    w_2 : float
        Mass fraction of component 2.
    mu_1 : float
        Molar mass of component 1 (kg/mol).
    mu_2 : float
        Molar mass of component 2 (kg/mol).

    Returns
    -------
    float
        Mole fraction of component 1.
    """
    n1 = w_1 / mu_1
    n2 = w_2 / mu_2
    total = n1 + n2
    if total <= 0:
        return 0.0
    return n1 / total

def rogers2025_binodal_temperature(x_H2, P_GPa):
    """Analytic binodal temperature for the H2-MgSiO3 system.

    Piecewise generalized logistic fit from Rogers+2025 Eqs. A1-A11.

    Parameters
    ----------
    x_H2 : float
        H2 mole fraction in the H2-MgSiO3 binary system.
    P_GPa : float
        Pressure in GPa.

    Returns
    -------
    float
        Binodal temperature in K. Returns 0 for x_H2 <= 0 or x_H2 >= 1.
    """
    if x_H2 <= 0.0 or x_H2 >= 1.0:
        return 0.0

    T_c = _rogers2025_T_critical(P_GPa)
    if T_c <= 0:
        # At P > 35 GPa, T_c < 0: always miscible
        return 0.0

    # Evaluate both branches and take the minimum (lower envelope).
    # The ascending branch rises from x=0 toward T_c; the descending
    # branch falls from T_c toward x=1. Their natural crossing defines
    # the binodal peak. Using min() avoids the artificial kink that a
    # hard if/else cutoff at x_c would produce.
    arg_asc = _R25_ALPHA_3 * (x_H2 - _R25_ALPHA_4)
    arg_asc = max(min(arg_asc, 500.0), -500.0)
    denom_asc = (1.0 + _R25_ALPHA_2 * math.exp(-arg_asc)) ** (1.0 / _R25_ALPHA_5)
    T_asc = T_c / denom_asc

    arg_desc = _R25_BETA_3 * (x_H2 - _R25_BETA_4)
    arg_desc = max(min(arg_desc, 500.0), -500.0)
    denom_desc = (1.0 + _R25_BETA_2 * math.exp(-arg_desc)) ** (1.0 / _R25_BETA_5)
    T_desc = T_c / denom_desc

    return min(T_asc, T_desc)

def rogers2025_suppression_weight(P_Pa, T_K, w_H2, w_sil, T_scale=50.0):
    """Sigmoid suppression weight for H2-MgSiO3 miscibility.

    Returns ~1 when the system is above the binodal (H2 miscible with
    silicate) and ~0 when below (H2 immiscible, should be suppressed).

    Parameters
    ----------
    P_Pa : float
        Pressure in Pa.
    T_K : float
        Temperature in K.
    w_H2 : float
        Mass fraction of H2 in the binary system.
    w_sil : float
        Mass fraction of silicate (MgSiO3) in the binary system.
    T_scale : float
        Sigmoid width in K. Default 50 K gives a ~200 K transition zone.

    Returns
    -------
    float
        Weight in [0, 1].
    """
    if w_H2 <= 0 or w_sil <= 0:
        return 1.0

    x_H2 = mass_to_mole_fraction(w_H2, w_sil, MU_H2, MU_MGSIO3)
    P_GPa = P_Pa * 1e-9

    T_binodal = rogers2025_binodal_temperature(x_H2, P_GPa)

    if T_binodal <= 0:
        # No binodal (pure endmember or P > 35 GPa): always miscible
        return 1.0

    if T_scale <= 0:
        # Hard cutoff
        return 1.0 if T_K >= T_binodal else 0.0

    arg = (T_K - T_binodal) / T_scale
    arg = max(min(arg, 500.0), -500.0)
    return 1.0 / (1.0 + math.exp(-arg))

def gupta2025_gibbs_mixing(x_H2, T, P_GPa):
    """Gibbs free energy of H2-H2O mixing per mole (Eq. A2).

    G_mix = RT [y ln(y) + (1-y) ln(1-y)] + W * y * (1-y)

    where y = x / (x + lambda * (1-x)) is the transformed composition
    (Eq. A3) accounting for asymmetric mixing.

    Parameters
    ----------
    x_H2 : float
        H2 mole fraction.
    T : float
        Temperature in K.
    P_GPa : float
        Pressure in GPa.

    Returns
    -------
    float
        Gibbs free energy of mixing in J/mol.
    """
    if x_H2 <= 0 or x_H2 >= 1:
        return 0.0

    lam = _gupta2025_lambda(T)
    # Transformed composition (Eq. A3)
    denom = x_H2 + lam * (1.0 - x_H2)
    if denom <= 0:
        return 0.0
    y = x_H2 / denom

    if y <= 0 or y >= 1:
        return 0.0

    W = _gupta2025_W(T, P_GPa)

    # Ideal mixing + asymmetric Margules
    G = R_GAS * T * (y * np.log(y) + (1.0 - y) * np.log(1.0 - y)) + W * y * (1.0 - y)
    return G

def gupta2025_critical_pressure(T):
    """Critical pressure at temperature T (Eq. A8).

    P_c = [2RT - (W_H - T * W_S)] / W_V(T)

    Above this pressure, the system is in a single phase for all compositions.

    Parameters
    ----------
    T : float
        Temperature in K.

    Returns
    -------
    float
        Critical pressure in GPa. May be negative (system is always
        single-phase at this T for all P > 0).
    """
    W_V = _gupta2025_W_V(T)
    if abs(W_V) < 1e-30:
        return np.inf
    numerator = 2.0 * R_GAS * T - (_G25_W_H - T * _G25_W_S)
    return numerator / W_V

def gupta2025_critical_temperature(P_GPa, T_bounds=(300.0, 6000.0)):
    """Critical temperature at given pressure (inverse of P_c(T)).

    Uses a precomputed interpolation table for speed. Falls back to
    brentq for pressures outside the table range.

    Parameters
    ----------
    P_GPa : float
        Pressure in GPa.
    T_bounds : tuple of float
        Search interval for brentq fallback. Default (300, 6000).

    Returns
    -------
    float or None
        Critical temperature in K, or None if no root exists.
    """
    if P_GPa <= 0:
        return None
    log_p = math.log10(P_GPa)
    if _G25_TCRIT_LOG_P[0] <= log_p <= _G25_TCRIT_LOG_P[-1]:
        T_c = float(np.interp(log_p, _G25_TCRIT_LOG_P, _G25_TCRIT_VALS))
        if np.isfinite(T_c):
            return T_c
    # Fallback for out-of-range pressures
    return _gupta2025_critical_temperature_brentq(P_GPa, T_bounds)

def gupta2025_coexistence_compositions(T, P_GPa):
    """Coexisting H2 mole fractions from common-tangent construction (Eq. A7).

    At temperatures below the critical curve, the system separates into
    two phases: an H2-rich phase and an H2O-rich phase. This function
    finds the compositions of both phases.

    Parameters
    ----------
    T : float
        Temperature in K.
    P_GPa : float
        Pressure in GPa.

    Returns
    -------
    tuple of (float, float) or None
        (x_H2O_rich, x_H2_rich) mole fractions of H2 in each phase,
        or None if the system is single-phase (above critical curve).
    """
    P_c = gupta2025_critical_pressure(T)
    if P_GPa >= P_c:
        # Above critical pressure: single phase
        return None

    # Compute Gibbs energy on a fine grid to locate the two minima
    # of the common tangent
    n_pts = 500
    x_arr = np.linspace(0.001, 0.999, n_pts)
    G_arr = np.array([gupta2025_gibbs_mixing(x, T, P_GPa) for x in x_arr])

    # Find local minima of G
    minima = []
    for i in range(1, n_pts - 1):
        if G_arr[i] < G_arr[i - 1] and G_arr[i] < G_arr[i + 1]:
            minima.append(i)

    if len(minima) < 2:
        # No two-phase region
        return None

    # Refine: find the common tangent between first and last minimum
    i1, i2 = minima[0], minima[-1]
    x1, x2 = x_arr[i1], x_arr[i2]

    # Common tangent condition: the line from (x1, G1) to (x2, G2)
    # must be tangent to G(x) at both endpoints. Use Brent on:
    # f(x1, x2) = dG/dx(x1) - (G(x2)-G(x1))/(x2-x1) = 0
    # and similarly for x2.

    def _dG_dx(x):
        """Numerical derivative of G_mix at x."""
        dx = 1e-6
        return (
            gupta2025_gibbs_mixing(x + dx, T, P_GPa) - gupta2025_gibbs_mixing(x - dx, T, P_GPa)
        ) / (2 * dx)

    def _tangent_residual(x_pair):
        """Residual for common tangent: slope at x1 = slope at x2 = chord slope."""
        xa, xb = x_pair
        if xa >= xb or xa <= 0 or xb >= 1:
            return [1e10, 1e10]
        Ga = gupta2025_gibbs_mixing(xa, T, P_GPa)
        Gb = gupta2025_gibbs_mixing(xb, T, P_GPa)
        chord = (Gb - Ga) / (xb - xa)
        return [_dG_dx(xa) - chord, _dG_dx(xb) - chord]

    # Simple iterative refinement from the grid minima
    from scipy.optimize import fsolve

    try:
        sol = fsolve(_tangent_residual, [x1, x2], full_output=True)
        x_sol = sol[0]
        if x_sol[0] < x_sol[1] and 0 < x_sol[0] < 1 and 0 < x_sol[1] < 1:
            return (float(x_sol[0]), float(x_sol[1]))
    except Exception:
        pass

    # Fallback: return grid minima
    if 0 < x1 < x2 < 1:
        return (float(x1), float(x2))

    return None

def gupta2025_suppression_weight(P_Pa, T_K, w_H2, w_H2O, T_scale=50.0):
    """Sigmoid suppression weight for H2-H2O miscibility.

    Returns ~1 when the system is above the critical curve (H2 miscible
    with H2O) and ~0 when below (immiscible, H2 should be suppressed).

    Parameters
    ----------
    P_Pa : float
        Pressure in Pa.
    T_K : float
        Temperature in K.
    w_H2 : float
        Mass fraction of H2.
    w_H2O : float
        Mass fraction of H2O.
    T_scale : float
        Sigmoid width in K. Default 50 K.

    Returns
    -------
    float
        Weight in [0, 1].
    """
    if w_H2 <= 0 or w_H2O <= 0:
        return 1.0

    # The Gupta+2025 model is parameterized for T >= 300 K.
    # At lower T, the W_V and lambda terms diverge (1/T^2 dependence),
    # producing unphysical critical temperatures. Return 0 (immiscible)
    # as a safe default: at T < 300 K, H2 and H2O are certainly phase-
    # separated.
    if T_K < 300.0:
        return 0.0

    P_GPa = P_Pa * 1e-9

    # The critical temperature at this pressure (fast interpolated lookup)
    T_crit = gupta2025_critical_temperature(P_GPa)

    if T_crit is None:
        # Could not find a critical temperature: assume miscible
        return 1.0

    # Floor at the H2O critical temperature (647 K). Below this, H2O
    # condenses to liquid/ice and the Margules Gibbs model (fitted to
    # supercritical DFT-MD data) is not valid. H2 and condensed H2O
    # are always immiscible, so the floor ensures correct suppression
    # for cold sub-Neptune models with condensed water layers.
    T_crit = max(T_crit, 647.0)

    if T_scale <= 0:
        return 1.0 if T_K >= T_crit else 0.0

    arg = (T_K - T_crit) / T_scale
    arg = max(min(arg, 500.0), -500.0)
    return 1.0 / (1.0 + math.exp(-arg))