zephyrus.dispatcher

def dispatch(inputs: EscapeInputs) -> EscapeResult:
    """Dispatch one atmospheric state to its escape regime and rate.

    Runs the fixed evaluation order documented in the module notes and
    returns an :class:`EscapeResult`. Raises ``ValueError`` on malformed
    input, and on a non-finite candidate rate, which is a defect rather
    than a physical state. Every physically posed state returns a labeled,
    finite, non-negative result whose per-species rates sum to the bulk rate.
    """
    inputs.validate()
    st = inputs.settings
    flags: dict = {}
    diag: dict = {
        'documentation': {
            'murray_clay_exponents': dict(dg.MURRAY_CLAY_EXPONENTS),
            'dayside_factors': dict(dg.DAYSIDE_FACTORS),
            'kn_band': kn.KN_BAND,
            'lambda_crit_band': bl.LAMBDA_BAND,
        }
    }
    if inputs.atm_converged is False:
        flags['stale_input'] = True

    r_hill = hy.hill_radius_periapsis(inputs.M_p, inputs.M_star, inputs.a, inputs.e)
    photo, f = photospheric_level(inputs.profile, st.P_photo)
    flags.update(f)
    r_xuv = photo['r']

    # The tidal reduction of the escape barrier, needed by both rate
    # branches that measure one: the energy-limited rate and the
    # luminosity cap on the bolometric residual.
    xi_ktide = r_hill / inputs.R_p
    # The factor has a double root at xi = 1 and the rates divide by it. At
    # and below the root the barrier is gone, so the rates are computed
    # without it (the smaller reading) and the Roche screen relabels the state.
    if st.tidal and xi_ktide <= 1.0:
        k_factor = 1.0
        flags['k_tide_undefined'] = True
    elif st.tidal:
        k_factor = hy.k_tide(xi_ktide)
    else:
        k_factor = 1.0

    # Step 1: bolometric candidate, computed at every point.
    # At the launch level, the photospheric surface the threshold is
    # calibrated on; R_p is the interior radius in a coupled call.
    lam_gate = bl.lambda_restricted(inputs.M_p, photo['r'], inputs.T_eq, photo['mmw'])
    bolo_rate, bolo = bl.bolometric_candidate(
        inputs.M_p,
        inputs.R_p,
        inputs.T_eq,
        inputs.kappa_photo,
        photo,
        inputs.F_int,
        lam_gate,
        st.lambda_crit,
        k_tide=k_factor,
        residual_mode=st.residual_mode,
    )
    # Candidate warnings are merged only if that candidate wins, so the flags
    # describe the dispatched rate; input-hygiene flags describe the state
    # and merge as they are found.
    bolo_flags = dict(bolo['flags'])
    diag['lambda_gate'] = lam_gate
    diag['bolometric'] = {k: v for k, v in bolo.items() if k != 'flags'}

    # Step 2: hydrodynamic candidate (always computed; it is cheap).
    channels = dict(
        cool_atomic=st.cool_atomic,
        cool_co2_band=st.cool_co2_band,
        cool_o_finestructure=st.cool_o_finestructure,
        cool_recombination=st.cool_recombination,
    )
    # Resolved once so the hydrostatic upper structure and the extended wind
    # base are built at the same exobase temperature.
    t_exo = _resolve_t_exo(inputs, channels)
    # The wind base sets only the XUV wind and a nozzle launched on its
    # column, so its flags travel with those candidates.
    base, base_flags = _resolve_wind_base(inputs, t_exo)
    elements = atomize(base['vmr'])
    f_mean = _wind_flux(inputs.F_xuv)
    t_wind, thermo = th.solve_wind_temperature(
        inputs.T_eq,
        base,
        elements,
        f_mean,
        col_co2=th.overlying_co2_column(base, inputs.M_p),
        **channels,
    )
    # Held aside and merged only if a hydrodynamic candidate wins; a losing
    # candidate's record stays in diag['hydrodynamic'].
    hydro_flags: dict = dict(base_flags)
    if thermo.get('clamped'):
        hydro_flags['thermostat_clamped'] = thermo['clamped']
    rr = hy.rr_chain(inputs.M_p, f_mean, base['r'], t_wind, elements)
    if rr['subcritical']:
        hydro_flags['subcritical_sonic'] = True

    eps = st.efficiency
    if st.efficiency_mode == 'caldiroli':
        eta_eff, cf = hy.caldiroli_efficiency(inputs.F_xuv, inputs.M_p, inputs.R_p, k_factor)
        hydro_flags.update(cf)
        if eta_eff is not None:
            # Their efficiency is defined against an R_p^3 rate geometry.
            eps = eta_eff * (inputs.R_p / r_xuv) ** 2
        else:
            hydro_flags['efficiency_fallback_fixed'] = True
    mdot_el = hy.el_rate(eps, inputs.F_xuv, inputs.R_p, r_xuv, inputs.M_p, k_factor)
    # Off, the chain is still evaluated: its sonic point feeds the Knudsen switch.
    mdot_rr = rr['mdot_rr'] if st.recombination_limit else math.inf
    # One particle per intercepted front photon, on the EL disk so the two
    # share one photon budget (Owen & Alvarez 2016, Eq. 10).
    e_ion = rr['hnu0_eV'] * ev2joule
    m_ion = rr['mu_plus_wind'] * m_p
    mdot_pl = hy.pl_rate(inputs.F_xuv, r_xuv, e_ion, m_ion) if st.photon_limit else math.inf
    candidates = {'EL': mdot_el, 'RR': mdot_rr, 'PL': mdot_pl}
    winner = min(candidates, key=candidates.get)
    mdot_hydro = candidates[winner]
    hydro_label = f'hydrodynamic:{winner}'
    diag['hydrodynamic'] = dict(
        mdot_el=mdot_el,
        mdot_rr=mdot_rr,
        mdot_pl=mdot_pl,
        photon_limit=st.photon_limit,
        recombination_limit=st.recombination_limit,
        efficiency_photon_limit=G * inputs.M_p * k_factor * m_ion / (e_ion * inputs.R_p),
        efficiency=eps,
        K_tide=k_factor,
        T_wind=t_wind,
        selection_mechanism=hy.selection_mechanism(rr, winner),
        rr_chain={
            k: rr[k]
            for k in (
                'c_s',
                'R_s',
                'R_s_calc',
                'lambda_b',
                'rho_s',
                'f_plus_base',
                'barometric_factor',
                'mu_wind',
                'mu_plus_wind',
                'hnu0_eV',
            )
        },
    )
    diag['thermostat'] = thermo

    # The L1 nozzle candidate (Jackson et al. 2017 Eq. 3). Either setting
    # launches from one level with one temperature and one mean mass, which
    # the Bernoulli invariance of the rate along an isothermal column needs.
    if st.nozzle_temperature == 'wind':
        t_nozzle = t_wind
        mu_nozzle = rr['mu_wind'] * m_p
        # Ideal-gas density at the base pressure for the wind's state, since
        # pressure is continuous across the temperature transition and
        # density is not.
        r_nozzle = base['r']
        rho_nozzle = base['p'] * mu_nozzle / (kb * t_nozzle)
        vmr_nozzle = elements
    else:
        t_nozzle, mu_nozzle = photo['T'], photo['mmw']
        r_nozzle, rho_nozzle = photo['r'], photo['rho']
        vmr_nozzle = photo['vmr']
    sigma_nozzle, _ = kn.sigma_mixture(vmr_nozzle, t_nozzle)
    nozzle_args = dict(
        rho_ph=rho_nozzle,
        r_ph=r_nozzle,
        T=t_nozzle,
        mu_kg=mu_nozzle,
        sigma_m2=sigma_nozzle,
    )
    nozzle_full, noz = nz.nozzle_candidate(
        inputs.M_p, inputs.M_star, inputs.a, inputs.e, kn_crit=st.kn_crit, **nozzle_args
    )
    # The competed rate is duty-cycled over the arc where overflow applies;
    # the unguarded average is kept so the closed form stays comparable
    # with the primary's published rates.
    nozzle_rate = noz['rate_applicable_kg_s']
    nozzle_applicable = noz['applicable']
    noz['rate_kg_s'] = nozzle_rate
    noz['rate_full_orbit_kg_s'] = nozzle_full
    noz['temperature_mode'] = st.nozzle_temperature
    noz['r_launch'] = r_nozzle
    noz['rho_launch'] = rho_nozzle
    # The uncapped flow's lift power is reported against what the planet
    # has, to show where the isothermal assumption is strained.
    noz['L_int_W'] = 4.0 * math.pi * inputs.R_p**2 * inputs.F_int
    noz['L_bol_intercepted_W'] = math.pi * inputs.R_p**2 * inputs.F_bol
    diag['nozzle'] = dict(noz)

    # Step 3: the sonic-point Knudsen switch.
    n_sc = rr['rho_s'] / (rr['mu_plus_wind'] * m_p)  # heavy-particle density
    if n_sc > _TINY:
        kn_sc, sigma_c, sigma_prov = kn.kn_sonic(
            n_sc, rr['R_s'], elements, t_wind, gamma=st.gamma_wind
        )
    else:
        kn_sc, sigma_c, sigma_prov = math.inf, math.nan, {}
    threshold = kn.effective_threshold(st.kn_crit, st.kn_hysteresis, inputs.prev_regime)
    diag['knudsen'] = dict(
        kn_sc=kn_sc,
        threshold_applied=threshold,
        sigma_c=sigma_c,
        provenance=sigma_prov,
    )
    if threshold != st.kn_crit:
        flags['hysteresis_active'] = True

    # Step 4: hydrostatic branch (always evaluated: its exobase quantities
    # feed the diagnostics at every dispatch).
    hs_per_element, hsd = hs.hydrostatic_rates_refined(
        inputs.profile,
        inputs.M_p,
        t_exo,
        gamma_bates=st.gamma_bates,
        kzz_default=st.kzz,
        n_levels_min=st.hydrostatic_levels_min,
        n_levels_max=st.hydrostatic_levels_max,
        rtol=st.hydrostatic_rtol,
    )
    hs_flags = hsd.pop('flags')
    mdot_hs = sum(hs_per_element.values())
    hnu, _e_ion, sigma_front = th.front_constants(elements)
    f_plus_exo = th.ionization_fraction(
        hsd['n_exo'] * 1e-6,
        hsd['T_exo'],
        f_mean / ergcm2stoWm2,
        hnu,
        sigma_front,
        th.recombination_alpha(elements, hsd['T_exo']),
    )
    unstable, contested = hs.gate_unstable(hsd['T_exo'], hsd, st.gate, f_plus_exo)
    diag['hydrostatic'] = dict(
        rate_kg_s=mdot_hs,
        convergence=hsd['convergence'],
        T_exo=hsd['T_exo'],
        T_exo_mode=st.T_exo_mode,
        r_exo=hsd['r_exo'],
        f_plus_exo=f_plus_exo,
        T_esc_neutral=hsd['T_esc_neutral'],
        T_esc_plasma=hsd['T_esc_plasma'],
        gate=st.gate,
        gate_unstable=unstable,
        detail=hsd,
    )
    if contested:
        flags['contested_ion'] = True
        diag['contested_ion'] = dict(
            hydrostatic_rate=mdot_hs,
            hydrodynamic_rate=mdot_hydro,
            note=(
                'the neutral and plasma escape-temperature conventions disagree '
                'here, so the branch assignment depends on ion physics this '
                'version does not model; both rates are recorded'
            ),
        )

    # Route. ``branch`` names the physics that produced the rate and decides
    # the split; ``label`` is what the caller reads, and the Roche screen
    # may overwrite it.
    gate_open = bolo['active']
    collisional = kn_sc <= threshold
    if not gate_open and not collisional and unstable:
        # The escape-temperature gate admits the wind beside the Jeans rate,
        # outside its own domain there (Volkov et al. 2011).
        flags['gate_rerouted'] = True
    admit = admitted(gate_open, collisional, unstable, bolo['competes'], nozzle_applicable)
    rates = dict(
        hydrostatic=mdot_hs,
        boiloff=bolo_rate,
        hydrodynamic=mdot_hydro,
        residual=bolo_rate,
        nozzle=nozzle_rate,
    )
    winner = select(admit, rates)
    # A nozzle win keeps the bound winner's flow radius: the lobe radius
    # would pin xi_flow at the fixed lobe-to-Hill ratio.
    bound = select(tuple(c for c in admit if c != 'nozzle'), rates)
    flow_radii = {
        'hydrostatic': hsd['r_exo'],
        'boiloff': bolo['R_sonic'],
        'hydrodynamic': max(r_xuv, rr['R_s']),
        'residual': bolo['R_sonic'],
    }
    flow_radius = flow_radii[bound]
    rate = rates[winner]
    per_species = dict(hs_per_element) if winner == 'hydrostatic' else None
    branch_of = {
        'hydrostatic': 'hydrostatic',
        'boiloff': 'boiloff',
        'hydrodynamic': hydro_label,
        'residual': 'boiloff',
        'nozzle': 'roche_overflow',
    }
    branch = branch_of[winner]
    if winner == 'nozzle':
        winner_flags = _nozzle_flags(noz, inputs.e)
        if st.nozzle_temperature == 'wind':
            # Launched on the wind's column at the thermostat's temperature.
            winner_flags.update(base_flags)
            if thermo.get('clamped'):
                winner_flags['thermostat_clamped'] = thermo['clamped']
    else:
        winner_flags = {
            'hydrostatic': hs_flags,
            'boiloff': bolo_flags,
            'hydrodynamic': hydro_flags,
            'residual': dict(bolo_flags, bolometric_residual=True),
        }[winner]
    # Only the winner's warnings reach the result.
    flags.update(winner_flags)
    label = branch

    # Step 6: the Roche screen renames the state and never touches the rate,
    # which keeps the dispatched rate continuous across the label boundary,
    # since the tested flow radius belongs to the branch that won the rate.
    xi_flow = r_hill / flow_radius if flow_radius > 0 else math.inf
    # The atmosphere's outer extent separates the two overflow geometries
    # but does not trigger the screen, which asks whether the escaping flow
    # stays bound; triggering on it would move the label boundary.
    r_atm = max(float(inputs.profile.r[-1]), hsd['r_exo'])
    diag['roche'] = dict(
        R_hill_periapsis=r_hill,
        flow_radius=flow_radius,
        xi_flow=xi_flow,
        xi_ktide=xi_ktide,
        k_tide=k_factor,
        tidal_inflation=1.0 / k_factor,
        r_atmosphere=r_atm,
        rate_branch=branch,
    )
    # Two routes into the label, geometric or a nozzle win, with a geometric
    # subflag under either; ``near_roche`` warns about a bound rate's tidal
    # inflation, so it is not raised on a labeled state.
    reaches_lobe = r_atm >= noz['r_lobe']
    renamed, subflag, near = roche_label(xi_flow, xi_ktide, winner == 'nozzle', reaches_lobe)
    if renamed:
        label = 'roche_overflow'
        flags['roche_overflow'] = True
        flags['roche_subflag'] = subflag
    elif near:
        flags['near_roche'] = True

    # The label at either edge of the collisionality band, through the same
    # admission, selection, and screen as the verdict itself. The nozzle is
    # re-evaluated there, since its fluid condition at L1 reads the threshold.
    counterfactual = {}
    for edge in kn.KN_BAND:
        _, noz_e = nz.nozzle_candidate(
            inputs.M_p, inputs.M_star, inputs.a, inputs.e, kn_crit=edge, **nozzle_args
        )
        rates_e = dict(rates, nozzle=noz_e['rate_applicable_kg_s'])
        adm = admitted(
            gate_open, kn_sc <= edge, unstable, bolo['competes'], noz_e['applicable']
        )
        bound_e = select(tuple(c for c in adm if c != 'nozzle'), rates_e)
        radius = flow_radii[bound_e]
        xi_e = r_hill / radius if radius > 0 else math.inf
        won = select(adm, rates_e)
        renamed_e = roche_label(xi_e, xi_ktide, won == 'nozzle', reaches_lobe)[0]
        counterfactual[edge] = 'roche_overflow' if renamed_e else branch_of[won]
    diag['knudsen']['counterfactual_labels'] = counterfactual

    # Per-species split, by the branch that produced the rate.
    if per_species is None:
        if branch.startswith('hydrodynamic') and st.fractionate and rate > 0.0:
            per_species, cdiag, cflags = closure_per_species(
                rate, elements, t_wind, inputs.M_p, base['r']
            )
            flags.update(cflags)
            diag['closure'] = cdiag
        else:
            per_species, sflags = unfractionated_split(rate, inputs.reservoirs, elements)
            flags.update(sflags)

    # Boxed diagnostics (nothing above this line reads them back).
    m_bar = mean_particle_mass(elements)
    ratio, q_net, q_c = dg.q_net_over_qc(
        eps,
        inputs.F_xuv,
        r_xuv,
        rr['R_s'],
        base['r'],
        inputs.M_p,
        m_bar,
        sigma_c if sigma_c == sigma_c else 1e-19,
    )
    diag['johnson_q'] = dict(q_net_over_qc=ratio, q_net_W=q_net, q_c_W=q_c)
    diag['guo_triple'] = dg.guo_triple(
        inputs.M_p,
        inputs.R_p,
        inputs.T_eq,
        photo['mmw'],
        inputs.M_star,
        inputs.a,
        inputs.e,
        hsd['lambda_exo_bulk'],
    )
    diag['potential_screens'] = dg.potential_screens(inputs.M_p, inputs.R_p)
    diag['erkaev_tc_K'] = dg.erkaev_tc(inputs.M_p, inputs.R_p, hsd['r_exo'], r_hill)
    diag['fluid_check'] = dg.along_profile_fluid_check(
        inputs.profile, inputs.M_p, rr['R_s'], st.kn_crit
    )
    # On the closed-form wind before the luminosity cap: that cap is L/(g R_p K)
    # and the check's cooling time is set by the same L, so testing it against
    # itself would return K rather than a verdict on the wind.
    wind_rate = min(bolo['mdot_parker'], bolo['mdot_bondi'])
    diag['tang_timescale'] = bl.tang_timescale_check(
        inputs.M_p, inputs.R_p, inputs.F_int, wind_rate, inputs.reservoirs
    )
    diag['self_consistency'] = dg.self_consistency_screen(inputs.reservoirs, rate, inputs.age)
    diag['rate_floor'] = dg.rate_floor_screen(rate)
    diag['base_level'] = dict(
        p_Pa=base['p'],
        p_physical_Pa=base.get('p_physical'),
        r_m=base['r'],
        T_K=base['T'],
        clamp_decades=base_flags.get('base_clamp_decades'),
        flags=dict(base_flags),
    )

    # Output contract: non-negative bulk rate, per-species rates summing to
    # it exactly.
    rate = max(rate, 0.0)
    per_species = {el: max(v, 0.0) for el, v in per_species.items()}
    tot = sum(per_species.values())
    if tot > 0.0 and rate > 0.0 and abs(tot - rate) / rate > 1e-9:
        per_species = {el: v * rate / tot for el, v in per_species.items()}
    return EscapeResult(
        regime=label,
        rate_branch=branch,
        mdot=rate,
        per_species=per_species,
        flags=flags,
        diagnostics=diag,
    )