zephyrus.knudsen

def lar_omega_star(pair: tuple, T: float, which: str = 'omega11') -> float:
    """Reduced collision integral Omega^(l,s)* (Laricchiuta et al. 2009, Eqs. 2-3)."""
    beta, eps0_mev, _re = LARICCHIUTA_PAIRS[pair]
    a1, a2, a3, a4, a5, a6, a7 = _lar_ai(beta, which)
    x = math.log(kb * T / (eps0_mev * 1e-3 * ev2joule))
    s1 = 1.0 / (1.0 + math.exp(-2.0 * (x - a3) / a4))
    s2 = 1.0 / (1.0 + math.exp(-2.0 * (x - a6) / a7))
    return math.exp((a1 + a2 * x) * s1 + a5 * s2)

def lar_sigma2_omega(pair: tuple, T: float, which: str = 'omega11') -> float:
    """Dimensional sigma^2 Omega^(l,s)* in Angstrom^2 (Laricchiuta et al. 2009, Eq. 4)."""
    beta, _eps0, re_a = LARICCHIUTA_PAIRS[pair]
    x0 = _XI1_M6 * beta**_XI2_M6
    return (x0 * re_a) ** 2 * lar_omega_star(pair, T, which)

def lar_sigma_diff(pair: tuple, T: float) -> float:
    """Momentum-transfer cross section pi sigma^2 Omega^(1,1)*, in m^2."""
    return math.pi * lar_sigma2_omega(pair, T, 'omega11') * 1e-20

def viscosity_pure(pair: tuple, mass_gmol: float, T: float) -> float:
    """First-approximation Chapman-Enskog viscosity of a pure gas, in Pa s.

    ``eta = 2.6693e-5 sqrt(M T) / (sigma^2 Omega^(2,2)*)`` poise, with the
    molar mass in g/mol and the collision integral in Angstrom^2 (the
    standard first Chapman-Enskog approximation; see e.g., Hirschfelder,
    Curtiss & Bird 1954). Exposed as the validation route: it anchors the
    Laricchiuta transcription on measured viscosities.
    """
    s2o22 = lar_sigma2_omega(pair, T, 'omega22')
    eta_poise = 2.6693e-5 * math.sqrt(mass_gmol * T) / s2o22
    return eta_poise * 0.1

def sigma_zk90_hydrogen(species: str, T: float) -> float:
    """H-H or H2-H2 momentum-transfer cross section in m^2 (Zahnle et al. 1990, Eq. 30).

    ``sigma_c = (3 sqrt(pi) / (16 b11)) sqrt(2 k T / mu11)`` in cgs, with
    ``mu11 = m / 2`` and ``b11 = 4.96e17 T^0.75 cm^-1 s^-1`` for H2-H2,
    scaled by 1.91 for H-H. Temperature dependence T^-0.25 by construction.
    """
    if species == 'H2':
        b11 = _B11_H2H2 * T**0.75
        mu11 = _M_H2_G / 2.0
    elif species == 'H':
        b11 = _H_COLUMN_SCALE * _B11_H2H2 * T**0.75
        mu11 = _M_H_G / 2.0
    else:
        raise ValueError(f'hydrogen route only covers H and H2, got {species!r}')
    sigma_cm2 = (3.0 * math.sqrt(math.pi) / (16.0 * b11)) * math.sqrt(2.0 * kb_cgs * T / mu11)
    return sigma_cm2 * 1e-4

def sigma_geometric(species: str) -> tuple[float, bool]:
    """Last-resort geometric hard-sphere cross section pi (2 r_vdW)^2, in m^2.

    Uses the Bondi (1964) van der Waals radius; for a composite molecule
    without a tabulated radius, the largest constituent-element radius sets
    the scale. A hard sphere has no temperature dependence, so against the
    shrinking collision integrals this fallback is roughly right at room
    temperature but overshoots by a factor of a few at 1e4 K, which biases
    the Knudsen number low and the switch toward hydrodynamic verdicts. The
    provenance class records which species sit on this fallback.

    Returns ``(sigma [m^2], tabulated)``, where ``tabulated`` is False when
    the radius came from ``FALLBACK_VDW_RADIUS_A`` rather than the published
    table, so the caller can class the two apart.
    """
    base = species.split('_')[0]
    if base in BONDI_VDW_RADIUS_A:
        return math.pi * (2.0 * BONDI_VDW_RADIUS_A[base] * 1e-10) ** 2, True
    radii = [BONDI_VDW_RADIUS_A.get(el) for el in parse_formula(base)]
    tabulated = all(r is not None for r in radii)
    r_max = max(r if r is not None else FALLBACK_VDW_RADIUS_A for r in radii)
    return math.pi * (2.0 * r_max * 1e-10) ** 2, tabulated

def sigma_species(species: str, T: float) -> tuple[float, str]:
    """Cross section for one species with its provenance class.

    The fallback order, most trusted source first: the Laricchiuta et al. (2009)
    like-pair collision integral where tabulated; the Zahnle et al. (1990)
    diffusion-inversion route for H and H2; the geometric Bondi-radius
    hard sphere as last resort. Returns ``(sigma [m^2], provenance)`` with
    provenance one of ``'laricchiuta'``, ``'zk90-scaled'``,
    ``'geometric-vdw'``, or ``'geometric-assumed-radius'`` when even the van
    der Waals radius was assumed rather than published.
    """
    base = species.split('_')[0]
    if (base, base) in LARICCHIUTA_PAIRS:
        return lar_sigma_diff((base, base), T), 'laricchiuta'
    if base in ('H', 'H2'):
        return sigma_zk90_hydrogen(base, T), 'zk90-scaled'
    sigma, tabulated = sigma_geometric(base)
    return sigma, 'geometric-vdw' if tabulated else 'geometric-assumed-radius'

def sigma_mixture(vmr: dict[str, float], T: float) -> tuple[float, dict]:
    """Density-weighted effective cross section of a mixture.

    ``sigma_C = sum_k n_k sigma_k / n``, the mixture weighting of
    Chatterjee & Pierrehumbert (2026, ApJ 998, 236, their Eq. 25).
    ``vmr`` maps species to mole fractions, renormalized internally over the
    species that are actually present: the normalization runs over the same
    positive weights the sum runs over, so a negative entry cannot shrink the
    denominator and inflate the result. A mixture with nothing present raises.
    Returns ``(sigma_C [m^2], provenance dict per species)``.
    """
    present = {sp: x for sp, x in vmr.items() if x > 0.0}
    tot = sum(present.values())
    if tot <= 0.0:
        raise ValueError('no species has a positive mole fraction')
    sig = 0.0
    prov: dict[str, str] = {}
    for sp, x in present.items():
        s, p = sigma_species(sp, T)
        sig += (x / tot) * s
        prov[sp] = p
    return sig, prov

def mean_free_path(sigma: float, n: float) -> float:
    """Maxwell mean free path 1/(sqrt(2) sigma n), in m."""
    return 1.0 / (SQRT2 * sigma * n)

def sonic_scale_height(r_sc: float, gamma: float = 1.0) -> float:
    """Analytic sonic-point density scale height, in m (CP26 Eq. 17)."""
    return (1.0 + gamma) * r_sc / (4.0 + SQRT2 * math.sqrt(5.0 - 3.0 * gamma))

def kn_sonic(
    n_sc: float, r_sc: float, vmr: dict[str, float], T_sc: float, gamma: float = 1.0
) -> tuple[float, float, dict]:
    """Sonic-point Knudsen number (CP26 Eq. 18) for a mixture.

    ``n_sc`` is the heavy-particle number density at the sonic point
    [m^-3], ``r_sc`` the sonic radius [m], ``vmr`` the (atomized) mixture
    composition, ``T_sc`` the temperature the cross sections are evaluated
    at [K], and ``gamma`` the polytropic index (1 for an isothermal wind).
    Returns ``(Kn_sc, sigma_C [m^2], provenance dict)``.

    At or below the threshold the sonic point is collisional and the wind
    hydrodynamic; above it escape is hydrostatic. The cross sections are
    neutral-neutral, the source's own onset, and ion channels would lower
    ``Kn_sc``, so a hydrostatic verdict on an ionized wind is the side
    fuller physics could overturn.
    """
    sigma, prov = sigma_mixture(vmr, T_sc)
    kn = (4.0 + SQRT2 * math.sqrt(5.0 - 3.0 * gamma)) / (
        SQRT2 * (1.0 + gamma) * sigma * n_sc * r_sc
    )
    return kn, sigma, prov

def effective_threshold(kn_crit: float, kn_hysteresis: float, prev_regime: str | None) -> float:
    """The switch threshold, with a hysteresis window for evolutionary use.

    ``prev_regime`` is the previous result's ``rate_branch``, the label of
    the branch whose rate won before the Roche screen renamed it, since a
    wind renamed ``roche_overflow`` is still a wind. With none the threshold
    is ``kn_crit`` (sharp switch). The threshold moves to
    ``kn_crit * kn_hysteresis`` while leaving a hydrodynamic wind and to
    ``kn_crit / kn_hysteresis`` while leaving a hydrostatic state, so a
    time-stepping track cannot chatter between branches on numerical noise.
    After a boil-off or an L1 transfer the switch is sharp, since neither
    sits on one side of it.
    """
    if prev_regime is None:
        return kn_crit
    if prev_regime.startswith('hydrodynamic'):
        return kn_crit * kn_hysteresis
    if prev_regime == 'hydrostatic':
        return kn_crit / kn_hysteresis
    return kn_crit