Zalmoxis.melting curves

def get_melting_curve_function(curve_id):
    """Return a callable melting curve f(P [Pa]) -> T [K] for the given identifier.

    Parameters
    ----------
    curve_id : str
        Melting curve identifier. See module docstring for valid values.

    Returns
    -------
    callable
        Function accepting pressure in Pa (scalar or array) and returning
        temperature in K.

    Raises
    ------
    ValueError
        If ``curve_id`` is not recognized.
    """
    if curve_id == 'Monteux16-solidus':
        return monteux16_solidus

    elif curve_id == 'Monteux16-liquidus-A-chondritic':

        def _liq(P):
            return monteux16_liquidus(P, model='A-chondritic')

        return _liq

    elif curve_id == 'Monteux16-liquidus-F-peridotitic':

        def _liq(P):
            return monteux16_liquidus(P, model='F-peridotitic')

        return _liq

    elif curve_id == 'Stixrude14-solidus':
        return stixrude14_solidus

    elif curve_id == 'Stixrude14-liquidus':
        return stixrude14_liquidus

    elif curve_id == 'PALEOS-liquidus':
        return paleos_liquidus

    elif curve_id == 'Monteux600-solidus-tabulated':
        return _load_tabulated_curve(
            os.path.join(ZALMOXIS_ROOT, 'data', 'melting_curves_Monteux-600', 'solidus.dat')
        )

    elif curve_id == 'Monteux600-liquidus-tabulated':
        return _load_tabulated_curve(
            os.path.join(ZALMOXIS_ROOT, 'data', 'melting_curves_Monteux-600', 'liquidus.dat')
        )

    else:
        all_valid = sorted(VALID_SOLIDUS | VALID_LIQUIDUS)
        raise ValueError(f"Unknown melting curve '{curve_id}'. Valid values: {all_valid}")

def get_solidus_liquidus_functions(
    solidus_id='Stixrude14-solidus',
    liquidus_id='Stixrude14-liquidus',
):
    """Load solidus and liquidus functions by config identifier.

    Parameters
    ----------
    solidus_id : str
        Solidus curve identifier.
    liquidus_id : str
        Liquidus curve identifier.

    Returns
    -------
    tuple of callable
        ``(solidus_func, liquidus_func)``
    """
    return get_melting_curve_function(solidus_id), get_melting_curve_function(liquidus_id)

def monteux16_liquidus(P, model='A-chondritic'):
    """Monteux+2016 liquidus temperature.

    Piecewise: Eq. 11 (P <= 20 GPa) and Eq. 13 (P > 20 GPa).

    Parameters
    ----------
    P : float or array-like
        Pressure in Pa (must be >= 0).
    model : str
        ``'A-chondritic'`` or ``'F-peridotitic'``.

    Returns
    -------
    float or ndarray
        Liquidus temperature in K.
    """
    P_arr = np.atleast_1d(np.asarray(P, dtype=float))
    T = np.empty_like(P_arr)
    lo = P_arr <= _P_TRANSITION
    hi = ~lo
    T[lo] = _liquidus_low(P_arr[lo])

    if model == 'A-chondritic':
        T[hi] = _liquidus_high_A(P_arr[hi])
    elif model == 'F-peridotitic':
        T[hi] = _liquidus_high_F(P_arr[hi])
    else:
        raise ValueError(f"Unknown model '{model}'. Use 'A-chondritic' or 'F-peridotitic'.")

    return float(T[0]) if np.ndim(P) == 0 else T

def monteux16_solidus(P):
    """Monteux+2016 solidus temperature.

    Piecewise: Eq. 10 (P <= 20 GPa) and Eq. 12 (P > 20 GPa).
    Both composition models share the same solidus.

    Parameters
    ----------
    P : float or array-like
        Pressure in Pa (must be >= 0).

    Returns
    -------
    float or ndarray
        Solidus temperature in K.
    """
    P_arr = np.atleast_1d(np.asarray(P, dtype=float))
    T = np.empty_like(P_arr)
    lo = P_arr <= _P_TRANSITION
    hi = ~lo
    T[lo] = _solidus_low(P_arr[lo])
    T[hi] = _solidus_high(P_arr[hi])
    return float(T[0]) if np.ndim(P) == 0 else T

def paleos_liquidus(P):
    """PALEOS MgSiO3 liquidus from Belonoshko+2005 / Fei+2021.

    Piecewise Simon-Glatzel fit:

    - P < 2.55 GPa: T = 1831 * (1 + P/4.6)^0.33  (Belonoshko+2005)
    - P >= 2.55 GPa: T = 6000 * (P/140)^0.26      (Fei+2021)

    This is the melting curve used internally by the PALEOS unified EOS
    tables. Using it for the mushy zone calculation ensures consistency
    between the table's phase boundaries and the derived solidus.

    Parameters
    ----------
    P : float or array-like
        Pressure in Pa (must be >= 0).

    Returns
    -------
    float or ndarray
        Liquidus temperature in K.
    """
    P_arr = np.atleast_1d(np.asarray(P, dtype=float))
    P_GPa = P_arr * 1e-9
    T = np.where(
        P_GPa < _PALEOS_P0_GPA,
        1831.0 * (1.0 + P_GPa / 4.6) ** 0.33,
        6000.0 * (P_GPa / 140.0) ** 0.26,
    )
    # Guard P=0: avoid 0^0.26 = NaN
    T = np.where(P_arr > 0, T, 0.0)
    return float(T[0]) if np.ndim(P) == 0 else T

def stixrude14_liquidus(P):
    """Silicate (MgSiO3) liquidus from Stixrude (2014) Eq. 1.9.

    Simon-like power law fitted to Lindemann melting theory:
    T_rock = 5400 K * (P / 140 GPa)^0.480

    Defined for all P > 0. At P = 0, returns 0 K (extrapolation).

    Parameters
    ----------
    P : float or array-like
        Pressure in Pa (must be >= 0).

    Returns
    -------
    float or ndarray
        Liquidus temperature in K.
    """
    P_arr = np.atleast_1d(np.asarray(P, dtype=float))
    # Guard P=0: the power law gives 0 K at P=0
    T = np.where(
        P_arr > 0,
        _STIX14_T_REF * (P_arr / _STIX14_P_REF) ** _STIX14_EXPONENT,
        0.0,
    )
    return float(T[0]) if np.ndim(P) == 0 else T

def stixrude14_solidus(P):
    """Silicate solidus from Stixrude (2014) Eqs. 1.9 + 1.10.

    The solidus is the pure-substance liquidus (Eq. 1.9) depressed by
    the cryoscopic equation (Eq. 1.10) with x0 = 0.79 (mole fraction
    of pure MgSiO3 in an Earth-like mantle composition):

    T_solidus = T_liquidus * (1 - ln(x0))^{-1}

    With x0 = 0.79 this gives ~4370 K at 140 GPa, approximately
    consistent with the experimentally measured solidus of an Earth-like
    mantle composition (~4100 K at 140 GPa, Fiquet et al. 2010). The
    ~270 K offset is within the uncertainty of the cryoscopic model.

    Parameters
    ----------
    P : float or array-like
        Pressure in Pa (must be >= 0).

    Returns
    -------
    float or ndarray
        Solidus temperature in K.
    """
    return stixrude14_liquidus(P) * _STIX14_CRYO_FACTOR