zephyrus.boiloff

def lambda_restricted(M_p: float, R_p: float, T_eq: float, mu_kg: float) -> float:
    """Restricted Jeans parameter Lambda = G M_p mu / (kB T_eq R_p).

    Dimensionless (Fossati et al. 2017). ``R_p`` [m] is the radius of the
    level it is evaluated at and ``mu_kg`` [kg] the mean molecular mass
    there; the dispatcher passes the launch level for both. For isothermal
    gas this equals ``2 R_B / R_p`` with the Bondi radius at ``T_eq``.
    """
    return G * M_p * mu_kg / (kb * T_eq * R_p)

def parker_mach(x: float) -> float:
    """Photospheric Mach number of the transonic isothermal Parker wind.

    Exact Lambert-W form of Owen & Wu (2016): ``Mach = sqrt(-W0(-f))`` with
    ``f = x^-4 exp(3 - 4/x)`` and ``x = R_launch / R_B <= 1``. At ``x = 1``
    the launch level sits at the sonic (Bondi) radius and the Mach number
    is 1; for small ``x`` it shuts off exponentially.

    Raises
    ------
    ValueError
        If ``x`` is outside ``(0, 1]``; the caller clamps inflated
        configurations to 1 with a flag before calling.
    """
    if not 0.0 < x <= 1.0:
        raise ValueError('parker_mach needs 0 < x <= 1')
    f = x**-4 * math.exp(3.0 - 4.0 / x)
    # f -> 1/e as x -> 1 (sonic point at the Bondi radius, Mach 1); clamp at
    # the W0 branch point against floating-point overshoot.
    if f >= math.exp(-1.0) * (1.0 - 1e-12):
        return 1.0
    w = lambertw(-f, 0)
    return math.sqrt(max(-w.real, 0.0))

def bolometric_candidate(
    M_p: float,
    R_p: float,
    T_eq: float,
    kappa_photo: float,
    launch: dict,
    F_int: float,
    lambda_gate: float,
    lambda_crit: float,
    k_tide: float = 1.0,
    residual_mode: str = 'off',
) -> tuple[float, dict]:
    """The bolometrically driven candidate mass-loss rate, in kg/s.

    Computed at every call: while ``lambda_gate < lambda_crit`` the
    atmosphere is inflated enough to boil off and the candidate is
    ``min(Parker rate, Bondi cap)``; past the gate the same machinery stays
    alive as a residual, additionally capped by the interior luminosity
    (see the module provenance notes).

    Parameters
    ----------
    M_p, R_p : float
        Planet mass [kg] and radius [m].
    T_eq : float
        Equilibrium temperature [K]; the wind runs at ``T_eq / 2^(1/4)``.
    kappa_photo : float
        Photospheric opacity [m^2 kg^-1]; the Parker rate scales as its
        inverse.
    launch : dict
        The photospheric working level (from
        :func:`zephyrus.profiles.photospheric_level`): uses ``p``, ``r``,
        ``mmw`` (molecular, not atomized), and ``rho``.
    F_int : float
        Interior heat flux [W m^-2], for the luminosity cap.
    lambda_gate : float
        The restricted Jeans parameter of the configuration.
    lambda_crit : float
        The activation threshold (20 by default upstream; band 15 to 35).
    k_tide : float
        Erkaev tidal reduction factor of the escape barrier, evaluated at
        ``xi = R_Hill / R_p``. Divides the luminosity cap, which is the
        only term that measures a barrier. The default of 1 is the untidal
        form; a non-positive value means the barrier has vanished and the
        cap is dropped.
    residual_mode : str
        ``'off'`` or ``'luminosity_capped'``: whether the candidate competes
        for the dispatched rate past the activation gate. It changes only
        the reported ``competes``, never the rate.

    Returns
    -------
    (rate, detail)
        The candidate rate [kg/s], the minimum over the caps in force, and
        a detail dict carrying the wind state, each cap, ``binding_cap``
        (which of them set the rate), the activation state ``active``,
        whether the candidate ``competes``, and flags (``bondi_inflated``
        when the launch level sits above the Bondi radius).

    Raises
    ------
    ValueError
        If ``residual_mode`` is not a supported value.
    """
    if residual_mode not in ('off', 'luminosity_capped'):
        raise ValueError("residual_mode must be 'off' or 'luminosity_capped'")
    T_w = T_eq / 2.0**0.25
    mu = launch['mmw']
    c_s = math.sqrt(kb * T_w / mu)
    R_B = G * M_p / (2.0 * c_s**2)
    R_launch = launch['r']
    x = R_launch / R_B
    flags = {}
    if x > 1.0:
        flags['bondi_inflated'] = True  # photosphere above the sonic radius
        x = 1.0
    mach = parker_mach(x)
    mdot_parker = 4.0 * math.pi * G * M_p * mach / (kappa_photo * c_s)

    rho_launch = launch['rho']
    mdot_bondi = (
        4.0 * math.pi * R_B**2 * c_s * rho_launch * math.exp(2.0 - 2.0 * R_B / R_launch)
    )

    # Plane-parallel optical depth of the launch level, reported only: tau
    # far from 1 means the prescribed level and the opacity do not describe
    # the photosphere the Parker rate assumes.
    g_launch = G * M_p / R_launch**2
    tau_launch = kappa_photo * launch['p'] / g_launch

    active = lambda_gate < lambda_crit
    caps = {'parker': mdot_parker, 'bondi': mdot_bondi}
    mdot_lum = None
    if not active:
        L = 4.0 * math.pi * R_p**2 * F_int
        g = G * M_p / R_p**2
        barrier = g * R_p * k_tide  # J/kg to lift gas out, tides included
        mdot_lum = L / barrier if barrier > 0.0 else math.inf
        caps['luminosity'] = mdot_lum
    binding_cap = min(caps, key=caps.get)
    rate = caps[binding_cap]
    if binding_cap == 'luminosity':
        # Flagged because the cap switches on at the gate, so a state crossing
        # the threshold drops discontinuously while keeping the same label.
        flags['luminosity_capped'] = True
    return rate, dict(
        T_wind=T_w,
        c_s=c_s,
        R_B=R_B,
        x=min(R_launch / R_B, 1.0),
        mach=mach,
        mdot_parker=mdot_parker,
        mdot_bondi=mdot_bondi,
        mdot_luminosity=mdot_lum,
        k_tide=k_tide,
        active=active,
        competes=active or residual_mode == 'luminosity_capped',
        residual_mode=residual_mode,
        binding_cap=binding_cap,
        rate_kg_s=rate,
        R_sonic=R_B,
        tau_launch=tau_launch,
        p_launch=float(launch['p']),
        flags=flags,
    )

def tang_timescale_check(
    M_p: float, R_p: float, F_int: float, mdot: float, reservoirs: dict | None
) -> dict:
    """Boil-off termination diagnostic after Tang et al. (2024, Eq. 8).

    Boil-off has ended once the mass-loss timescale ``t_Mdot = M_env/Mdot``
    reaches the cooling timescale ``t_cool = G M_p M_env / (R_p L)`` with
    ``L = 4 pi R_p^2 F_int`` (their Eq. 9 with the envelope concentration
    factor at 1 and the radiative-convective boundary at ``R_p``, both
    stated approximations). Returns ``{'evaluated': False}`` when no
    reservoir masses are supplied; otherwise the two timescales and the
    verdict. Diagnostic only: it never gates the rate.
    """
    if not reservoirs:
        return {'evaluated': False}
    M_env = sum(reservoirs.values())
    if M_env <= 0.0 or mdot <= 0.0:
        return {'evaluated': False}
    L = 4.0 * math.pi * R_p**2 * F_int
    t_mdot = M_env / mdot
    t_cool = G * M_p * M_env / (R_p * L)
    return {
        'evaluated': True,
        't_mdot_s': t_mdot,
        't_cool_s': t_cool,
        'terminated': t_mdot >= t_cool,
    }