Light-element core

def alloy_density(rho_fe, w, pressure, temperature):
    """Density of the light-element core at one (P, T) point.

    Parameters
    ----------
    rho_fe : float
        Density of the iron base at (P, T), kg/m^3.
    w : dict
        Light-element mass fractions from :func:`split_alloy`.
    pressure : float
        Pressure in Pa.
    temperature : float
        Temperature in K.

    Returns
    -------
    float or None
        Density in kg/m^3; None when the iron base or the Luo members have no value
        at (P, T).

    Raises
    ------
    CoreAlloyError
        A carrier that leaves no alloy mass, a composition outside the bounds of
        :func:`_linear_factor`, or a core without positive finite density.
    """
    w_c, lin = _linear_factor(w)
    if rho_fe is None or not np.isfinite(rho_fe) or rho_fe <= 0:
        return None
    rho_lin = rho_fe * (1.0 - lin)
    if w_c <= 0:
        return rho_lin
    dv_c = carrier_excess_volume(pressure, temperature)
    if dv_c is None:
        return None
    rho = 1.0 / ((1.0 - w_c) / rho_lin + w_c * dv_c)
    if not np.isfinite(rho) or rho <= 0:
        raise CoreAlloyError(f'core alloy density {rho!r} kg/m^3 for {w!r}')
    return rho

def alloy_density_batch(rho_fe, w, pressures, temperatures):
    """Vectorised :func:`alloy_density`; NaN where a point has no density."""
    rho = (alloy_density(r, w, p, t) for r, p, t in zip(rho_fe, pressures, temperatures))
    return np.array(list(rho), dtype=float)  # None becomes NaN

def carrier_excess_volume(pressure, temperature):
    """Excess specific volume of the H2O carrier in liquid iron, m^3 per kg of H2O.

    Evaluated at (P [Pa], T [K]) clamped into the Luo fit box; a clamped point is counted
    in :data:`LUO_OUT_OF_RANGE`. Returns None when a Luo member has no volume there.
    """
    p = min(max(pressure, LUO_P_RANGE[0]), LUO_P_RANGE[1])
    t = min(max(temperature, LUO_T_RANGE[0]), LUO_T_RANGE[1])
    if (p, t) != (pressure, temperature):
        if LUO_OUT_OF_RANGE['count'] == 0:
            logger.warning(
                'Core Fe-H-O carrier at P = %.3e Pa, T = %.0f K is outside the Luo et al. '
                '(2024) fit range (8000-14000 K, 50-1300 GPa); its volume is taken at the edge',
                pressure,
                temperature,
            )
        LUO_OUT_OF_RANGE['count'] += 1
    v_dry = luo_cell_volume(p, t, LUO_DRY)
    v_wet = luo_cell_volume(p, t, LUO_WET)
    if v_dry is None or v_wet is None:
        return None
    return (v_wet - v_dry) * _A3 / _M_CARRIER_CELL

def is_alloy_mixture(mixture) -> bool:
    """Return True when a layer mixture has a ``CoreAlloy:`` component."""
    return any(c.startswith(ALLOY_PREFIX) for c in mixture.components)

def luo_cell_volume(pressure, temperature, member):
    """Cell volume in A^3 at (P [Pa], T [K]) from Eq. 14, or None without a root.

    The root is bracketed in [0.15, 1.5] V0, which holds 50-1300 GPa at 8000-14000 K
    with a wide margin.
    """
    v0 = member['V0']
    lo, hi = 0.15 * v0, 1.5 * v0
    f_lo = luo_pressure(lo, temperature, member) - pressure
    f_hi = luo_pressure(hi, temperature, member) - pressure
    if not (np.isfinite(f_lo) and np.isfinite(f_hi)) or f_lo * f_hi > 0:
        return None
    return brentq(
        lambda v: luo_pressure(v, temperature, member) - pressure, lo, hi, xtol=1e-12 * v0
    )

def luo_density(pressure, temperature, member):
    """Density in kg/m^3 of a Luo et al. (2024) member at (P [Pa], T [K]), or None."""
    v = luo_cell_volume(pressure, temperature, member)
    if v is None:
        return None
    return member['cell_mass'] / (v * _A3)

def luo_pressure(volume, temperature, member):
    """Pressure of a Luo et al. (2024) liquid cell, Eq. 14.

    Parameters
    ----------
    volume : float
        Cell volume in A^3 (64 Fe, or 64 Fe + 8 H2O).
    temperature : float
        Temperature in K.
    member : dict
        :data:`LUO_DRY` or :data:`LUO_WET`.

    Returns
    -------
    float
        Pressure in Pa.
    """
    x = member['V0'] / volume
    x23 = x ** (2.0 / 3.0)
    cold = 1.5 * member['K'] * (x ** (7.0 / 3.0) - x ** (5.0 / 3.0))
    cold *= 1.0 + 0.75 * (member['Kp'] - 4.0) * (x23 - 1.0)
    thermal = (
        (temperature - LUO_T0) / 1000.0 * (member['a'] + member['b'] * x + member['c'] * x * x)
    )
    return (cold + thermal) * 1e9

def split_alloy(mixture):
    """Split a core mixture into its iron base and its light-element fractions.

    Parameters
    ----------
    mixture : LayerMixture
        Core mixture with one base component and ``CoreAlloy:<element>`` components.

    Returns
    -------
    base : str
        EOS name of the iron base.
    w : dict
        Mass fraction per element key of :data:`ALLOY_ELEMENTS` (absent keys are 0).

    Raises
    ------
    CoreAlloyError
        Unknown element, element given twice, both H forms, not exactly one base,
        a negative or non-finite fraction, or light elements that leave no iron.
    """
    base, w = [], {}
    for comp, frac in zip(mixture.components, mixture.fractions):
        if not math.isfinite(frac) or frac < 0:
            raise CoreAlloyError(f'core component {comp!r} has mass fraction {frac!r}')
        if not comp.startswith(ALLOY_PREFIX):
            base.append(comp)
            continue
        el = comp[len(ALLOY_PREFIX) :]
        if el not in ALLOY_ELEMENTS:
            raise CoreAlloyError(f'unknown core light element {el!r} in {comp!r}')
        if el in w:
            raise CoreAlloyError(f'core light element {el!r} given twice')
        w[el] = float(frac)
    if len(base) != 1:
        raise CoreAlloyError(f'a light-element core needs exactly one iron base, got {base!r}')
    if base[0] != 'PALEOS:iron':
        raise CoreAlloyError(
            f'a light-element core needs the PALEOS:iron base, got {base[0]!r}'
        )
    if 'H' in w and 'H-linear' in w:
        raise CoreAlloyError('CoreAlloy:H and CoreAlloy:H-linear exclude each other')
    if sum(w.values()) >= 1.0:
        raise CoreAlloyError(f'light-element mass fractions sum to {sum(w.values())!r} >= 1')
    w_c, _ = _linear_factor(w)
    if w_c > LUO_CARRIER_FRACTION and not _CELL_WARNED[0]:
        _CELL_WARNED[0] = True
        logger.warning(
            'Core Fe-H-O carrier mass fraction %.4f is above that of the Luo et al. (2024) '
            'wet cell (%.4f); its excess volume is used beyond the cell composition',
            w_c,
            LUO_CARRIER_FRACTION,
        )
    return base[0], w