def alpha_heating_power_from_n_D_n_T_T_keV_V(
    n_D_1e19: float, n_T_1e19: float, T_keV: float, V_m3: float
) -> float:
    r"""D-T alpha heating power with the rough $\langle\sigma v\rangle \propto T^2$ fit.

    $$P_\alpha = n_D\,n_T\,\langle\sigma v\rangle\,E_\alpha\,V, \qquad
      \langle\sigma v\rangle \approx 1.1\times10^{-24}\,T_{\mathrm{keV}}^2\ \mathrm{m^3/s}$$

    Parameters
    ----------
    n_D_1e19 : float
        Deuterium density [1e19 m^-3].
    n_T_1e19 : float
        Tritium density [1e19 m^-3].
    T_keV : float
        Ion temperature [keV].
    V_m3 : float
        Plasma volume [m^3].

    Returns
    -------
    float
        Alpha heating power [W].

    Assumptions
    -----------
    Flat profiles (the product of densities is taken at one temperature over
    the whole volume); all alpha energy deposited in the plasma.

    Validity
    --------
    Empirical fit.  The quadratic $\langle\sigma v\rangle$ is Wesson's
    interpolation of the D-T reactivity, accurate to ~10 % for $10 < T <
    20$ keV [1]_; outside that window use the Bosch-Hale parametrisation [2]_.

    Limitations
    -----------
    Irrelevant for a hydrogen device such as VEST, which runs hydrogen and
    occasionally helium and so sustains no D-T reaction; kept for power-balance
    completeness.  Tracked in #360 (Bosch-Hale replacement).

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 1.3.
    .. [2] H.-S. Bosch and G. M. Hale, Nucl. Fusion 32 (1992) 611, Table VII.
    """
    n_D = n_D_1e19 * 1e19  # m^-3
    n_T = n_T_1e19 * 1e19  # m^-3
    sigma_v = SIGMA_V_COEF * T_keV**2  # m^3/s (rough fit)
    return n_D * n_T * sigma_v * E_ALPHA * V_m3  # W

def aspect_ratio_from_a_R(a: float, R: float) -> float:
    r"""Aspect ratio $A = R/a = 1/\varepsilon$.

    $$A = \frac{R}{a}$$

    Parameters
    ----------
    a : float
        Minor radius [m].
    R : float
        Major radius [m].

    Returns
    -------
    float
        Aspect ratio [-].

    Notes
    -----
    Not to be confused with the elongation $\kappa$.
    """
    return R / a

def auxiliary_heating_power(P_aux: float,
                          eta_CD: float) -> Tuple[float, float]:
    r"""Split auxiliary power into heating and current-drive parts.

    $$P_{CD} = \frac{P_{aux}}{1 + \eta_{CD}}, \qquad P_{heat} = P_{aux} - P_{CD}$$

    Parameters
    ----------
    P_aux : float
        Total auxiliary power [W].
    eta_CD : float
        Current-drive efficiency figure from :func:`current_drive_efficiency` [-].

    Returns
    -------
    P_heat : float
        Part of the power counted as heating [W].
    P_CD : float
        Part of the power counted as current drive [W].

    Limitations
    -----------
    A bookkeeping split with the same unsourced normalisation as
    ``eta_CD``; the two parts always sum to ``P_aux``.
    """
    P_CD = P_aux / (1 + eta_CD)
    P_heat = P_aux - P_CD
    return P_heat, P_CD

def ballooning_alpha_from_volume(V, dV_dpsi, dp_dpsi, R0):
    r"""Normalised pressure gradient $\alpha$ of local ballooning theory on a general equilibrium.

    $$\alpha = -\frac{2\mu_0}{(2\pi)^2}\,\frac{dV}{d\psi}\,\frac{dp}{d\psi}\,
      \sqrt{\frac{V}{2\pi^2R_0}}$$

    Parameters
    ----------
    V : float or np.ndarray
        Volume enclosed by the surface, positive [m^3].
    dV_dpsi : float or np.ndarray
        $dV/d\psi$ with $\psi$ the poloidal flux **per radian** [m^3 rad/Wb].
    dp_dpsi : float or np.ndarray
        $dp/d\psi$ against the same per-radian flux [Pa rad/Wb].
    R0 : float
        Reference major radius naming the volume radius, positive; the
        magnetic axis, as GPEC.jl's local ballooning uses [m].

    Returns
    -------
    float or np.ndarray
        $\alpha$, positive where the pressure falls outward [-].

    Raises
    ------
    ValueError
        ``V`` or ``R0`` is not positive, or an input is not finite.

    Convention
    ----------
    Unlike $\hat s_V$, this depends on the flux label: $\psi$ must be the
    poloidal flux per radian, $|\nabla\psi| = RB_p$. With the full flux in Wb
    each derivative carries a $2\pi$ and $\alpha$ comes out $(2\pi)^2$ too
    small. The sign of $\psi$ cancels (both derivatives flip). In the
    large-aspect-ratio circular limit, $d\psi/dr = rB_0/q$ and
    $V = 2\pi^2R_0r^2$, it is exactly the Connor-Hastie-Taylor
    $\alpha = -2\mu_0R_0q^2p'(r)/B_0^2$ with $r$ the **minor** radius
    (``ballooning_alpha_from_p_B_R`` differentiates against the major radius
    instead). $\alpha \propto \sqrt{R_0}$ through $r_V$: Miller et al. use each
    surface's geometric centre instead of the axis, which at a spherical
    tokamak's aspect ratio moves $\alpha$ by several per cent, so state which
    $R_0$ when comparing.

    Physical interpretation
    -----------------------
    The pressure-gradient drive of high-$n$ ballooning, measured against the
    field-line bending it must overcome, without choosing a minor radius or a
    field strength: both enter through $dV/d\psi$.

    Assumptions
    -----------
    Nested surfaces; local (high-$n$) ballooning ordering. It is a
    normalisation, not a stability criterion: the boundary it is compared
    with depends on the shaping the reduced $s$-$\alpha$ model drops.

    References
    ----------
    .. [1] R. L. Miller, M. S. Chu, J. M. Greene, Y. R. Lin-Liu and
           R. E. Waltz, Phys. Plasmas 5 (1998) 973.
    .. [2] J. W. Connor, R. J. Hastie and J. B. Taylor, Phys. Rev. Lett. 40
           (1978) 396.
    """
    V = np.asarray(V, dtype=float)
    dV_dpsi = np.asarray(dV_dpsi, dtype=float)
    dp_dpsi = np.asarray(dp_dpsi, dtype=float)
    for name, value in (("V", V), ("dV_dpsi", dV_dpsi), ("dp_dpsi", dp_dpsi), ("R0", np.asarray(R0, dtype=float))):
        if not np.all(np.isfinite(value)):
            raise ValueError(f"{name} must be finite")
    if np.any(V <= 0.0) or not float(R0) > 0.0:
        raise ValueError("V and R0 must be positive")
    result = -2.0 * MU0 / (2.0 * np.pi) ** 2 * dV_dpsi * dp_dpsi * np.sqrt(V / (2.0 * np.pi ** 2 * float(R0)))
    return float(result) if result.ndim == 0 else result

def beta_normal_from_beta_tor(beta_tor: float,
                              a: float,
                              B0: float,
                              Ip: float) -> float:
    r"""Normalized beta (Troyon) from toroidal beta, minor radius, field and current.

    $$\beta_N = 100\,\beta_t \frac{a |B_0|}{|I_p[\mathrm{MA}]|}$$

    Parameters
    ----------
    beta_tor : float
        Toroidal beta [-].
    a : float
        Minor radius [m].
    B0 : float
        Vacuum toroidal field at ``r0`` [T].
    Ip : float
        Plasma current; converted to MA internally [A].

    Returns
    -------
    float
        Normalized beta [% m T/MA].

    Reduction
    ---------
    input: scalar_0d
    output: scalar_0d
    kind: normalization
    locality: global
    role: stability_coordinate

    References
    ----------
    .. [1] F. Troyon et al., Plasma Phys. Control. Fusion 26 (1984) 209.
    .. [2] IMAS Data Dictionary, ``equilibrium.time_slice[:].global_quantities.beta_normal``.
    """
    return 100 * float(beta_tor) * float(a) * abs(float(B0)) / abs(float(Ip) / 1e6)

def beta_normal_from_beta_volume(beta_volume: float,
                                 a: float,
                                 B0: float,
                                 Ip: float) -> float:
    r"""Normalized volume beta: the volume beta normalized like the Troyon beta_N.

    $$\langle \beta_N \rangle = 100\,\beta_B \frac{a |B_0|}{|I_p[\mathrm{MA}]|}$$

    Parameters
    ----------
    beta_volume : float
        Volume beta, $2\mu_0\langle p\rangle_V/\langle B^2\rangle_V$
        (:func:`beta_volume_from_p_B2`) [-].
    a : float
        Minor radius [m].
    B0 : float
        Vacuum toroidal field at ``r0`` [T].
    Ip : float
        Plasma current; converted to MA internally [A].

    Returns
    -------
    float
        Normalized volume beta [% m T/MA].

    Convention
    ----------
    Menard's $\langle\beta_N\rangle$: the normalization $a B_0 / I_p$ is the
    conventional one, only the beta differs. Which $B_0$ is the caller's: Menard
    et al. take the vacuum field at the plasma's geometric centre, VAFT's
    ``beta_normal`` the one at ``vacuum_toroidal_field.r0``; where the two radii
    differ, so does the value, by their ratio. It is not the conventional
    :func:`beta_normal_from_beta_tor`; at low aspect ratio the conventional
    $\beta_N$ of an optimized no-wall sequence nearly doubles (3.15 at A = 10 to
    5.85 at A = 1.25) while this one stays at 3.2 within 3 % [2]_, [3]_.

    Semantics
    ---------
    consumes: minor_radius, b_t, plasma_current

    References
    ----------
    .. [1] F. Troyon et al., Plasma Phys. Control. Fusion 26 (1984) 209.
    .. [2] J. E. Menard et al., Phys. Plasmas 11 (2004) 639, doi:10.1063/1.1640623.
    .. [3] J. E. Menard et al., "Unified ideal stability limits for advanced tokamak
           and spherical torus plasmas", PPPL-3779 (2003), Fig. 3.
    """
    return 100 * float(beta_volume) * float(a) * abs(float(B0)) / abs(float(Ip) / 1e6)

def beta_poloidal_from_circumference(p_average: float,
                                     Ip: float,
                                     length_pol: float) -> float:
    r"""Poloidal beta in the EFIT/OMFIT convention, normalized by the LCFS circumference.

    $$\beta_{p,\mathrm{circ}} = \frac{2\mu_0 \langle p \rangle_V}{B_{pa}^2}, \qquad
      B_{pa} = \frac{\mu_0 I_p}{L_{pol}}$$

    Parameters
    ----------
    p_average : float
        Volume-averaged plasma pressure [Pa].
    Ip : float
        Plasma current [A].
    length_pol : float
        Poloidal circumference of the last closed flux surface [m].

    Returns
    -------
    float
        Poloidal beta, circumference convention [-].

    Convention
    ----------
    **Not** the IMAS DD's ``beta_pol``, and not an estimate of it: this
    normalizes the volume-averaged pressure by the poloidal field implied by
    the LCFS *circumference* rather than by ``R_0 mu_0 Ip^2``.  The two differ
    by the geometric factor ``R_0 L_pol^2 / (2 V)`` -- 26% on the packaged
    kineticEfit reference, where this form reproduces the stored value to
    0.1% and the DD form does not.  It exists so the database summary can keep
    reporting what OMFIT reported; ``global_quantities.beta_pol`` stays
    DD-normative (see :func:`beta_poloidal_from_pressure_integral` and issue
    #318, which owns the sensitivity study of the two).

    References
    ----------
    .. [1] L. L. Lao, H. St. John, R. D. Stambaugh, A. G. Kellman and
           W. Pfeiffer, Nucl. Fusion 25 (1985) 1611 (EFIT definitions).
    """
    b_pa = MU0 * abs(float(Ip)) / float(length_pol)
    return 2 * MU0 * float(p_average) / b_pa**2

def beta_poloidal_from_pressure_integral(pressure_integral: float,
                                         R0: float,
                                         Ip: float) -> float:
    r"""Poloidal beta in the IMAS definition, from the pressure volume integral.

    $$\beta_p = \frac{4 \int p \, dV}{R_0 \mu_0 I_p^2}$$

    Parameters
    ----------
    pressure_integral : float
        Plasma pressure integrated over the plasma volume [Pa m^3].
    R0 : float
        Reference major radius the DD normalizes by [m].
    Ip : float
        Plasma current [A].

    Returns
    -------
    float
        Poloidal beta [-].

    Convention
    ----------
    This is the DD-normative ``beta_pol``, normalized by ``R_0 mu_0 Ip^2``.
    The EFIT/OMFIT circumference form is a different definition; see
    :func:`beta_poloidal_from_circumference`.

    Reduction
    ---------
    input: scalar_0d
    output: scalar_0d
    kind: dimensionless_normalization
    locality: global
    role: global_descriptor

    References
    ----------
    .. [1] IMAS Data Dictionary, ``equilibrium.time_slice[:].global_quantities.beta_pol``.
    .. [2] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Ch. 3, Equilibrium (definition of beta).
    """
    return 4 * float(pressure_integral) / (float(R0) * MU0 * float(Ip) ** 2)

def beta_t_from_n_T_B(
    n_m3: Union[float, np.ndarray],
    T_eV: Union[float, np.ndarray],
    B_t_T: Union[float, np.ndarray],
    output: str = "percent",
) -> Union[float, np.ndarray]:
    r"""Toroidal beta from density, temperature and field.

    $$\beta_t\,[\%] = 8.05\times10^{-23}\,\frac{n\,[\mathrm{m^{-3}}]\;T\,[\mathrm{eV}]}{B_t^2\,[\mathrm{T^2}]}$$

    Parameters
    ----------
    n_m3 : float or np.ndarray
        Electron density, strictly positive [m^-3].
    T_eV : float or np.ndarray
        Temperature, strictly positive [eV].
    B_t_T : float or np.ndarray
        Toroidal field, strictly positive [T].
    output : str, optional
        ``"percent"`` (default) or ``"fraction"`` [str].

    Returns
    -------
    float or np.ndarray
        Toroidal beta in percent, or as a fraction [%].

    Raises
    ------
    ValueError
        For non-finite or non-positive input, or an unknown ``output``.

    Convention
    ----------
    $8.05\times10^{-23} = 100\times2\mu_0\times2e$: the total pressure is
    taken as $2n_eT$ (equal electron and ion temperatures, $n_i = n_e$), and
    the default output is a percentage.  Verdoolaege's ITER example is
    reproduced by construction.

    References
    ----------
    .. [1] G. Verdoolaege et al., Nucl. Fusion 61 (2021) 076006, Sec. 2.
    .. [2] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011), Sec. 3.5.
    """
    n_arr = _validate_positive("n_m3", n_m3)
    t_arr = _validate_positive("T_eV", T_eV)
    b_arr = _validate_positive("B_t_T", B_t_T)
    beta_percent = 8.05e-23 * n_arr * t_arr / (b_arr**2)

    normalized = output.strip().lower()
    if normalized == "percent":
        return _maybe_scalar(beta_percent)
    if normalized == "fraction":
        return _maybe_scalar(beta_percent / 100.0)
    raise ValueError("output must be either 'percent' or 'fraction'")

def beta_toroidal_from_p_B0(p_average: float,
                            B0: float) -> float:
    r"""Toroidal beta from the volume-averaged pressure and the vacuum field.

    $$\beta_t = \frac{2\mu_0 \langle p \rangle_V}{B_0^2}$$

    Parameters
    ----------
    p_average : float
        Volume-averaged plasma pressure [Pa].
    B0 : float
        Vacuum toroidal field at ``r0`` (``equilibrium.vacuum_toroidal_field.b0``),
        as the DD requires [T].

    Returns
    -------
    float
        Toroidal beta [-].

    Reduction
    ---------
    input: scalar_0d
    output: scalar_0d
    kind: dimensionless_normalization
    locality: global
    role: global_descriptor

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Ch. 3, Equilibrium (definition of beta).
    .. [2] IMAS Data Dictionary, ``equilibrium.time_slice[:].global_quantities.beta_tor``.
    """
    return 2 * MU0 * float(p_average) / float(B0) ** 2

def beta_volume_from_p_B2(p_average: float,
                          B2_average: float) -> float:
    r"""Volume beta: the volume-averaged pressure over the volume-averaged total magnetic energy density.

    $$\beta_B = \frac{2\mu_0 \langle p \rangle_V}{\langle B^2 \rangle_V}
               = \frac{2\mu_0 \int p \, dV}{\int B^2 \, dV}$$

    Parameters
    ----------
    p_average : float
        Volume-averaged plasma pressure inside the last closed flux surface [Pa].
    B2_average : float
        Volume average of the total field squared, $B_R^2 + B_Z^2 + B_\phi^2$, over the
        same volume [T^2].

    Returns
    -------
    float
        Volume beta [-].

    Convention
    ----------
    A ratio of volume averages -- the beta Menard et al. attribute to Troyon -- not the average of the local ratio
    $\langle 2\mu_0 p / B^2 \rangle_V$, and not the toroidal beta
    (:func:`beta_toroidal_from_p_B0`), which divides by the vacuum field at one radius.
    The two agree at large aspect ratio and low beta; at low aspect ratio the $1/R$
    variation of $B_\phi$ and the poloidal field make $\langle B^2 \rangle_V$ differ from
    $B_0^2$, which is why Menard et al. normalize by it to compare aspect ratios.

    References
    ----------
    .. [1] F. Troyon et al., Plasma Phys. Control. Fusion 26 (1984) 209
           (beta as twice the pressure over the magnetic energy integrals).
    .. [2] J. E. Menard et al., Phys. Plasmas 11 (2004) 639, doi:10.1063/1.1640623
           (PPPL-3908): definition of the volume-averaged total-field beta.
    """
    B2_average = float(B2_average)
    if not B2_average > 0:
        raise ValueError(f"B2_average must be positive, got {B2_average!r}")
    return 2 * MU0 * float(p_average) / B2_average

def bootstrap_current_fraction(n_e: float,
                             T_e_keV: float,
                             R0: float,
                             a: float,
                             q_95: float) -> float:
    r"""Heuristic bootstrap-current fraction $f_{BS}$.

    $$f_{BS} = 0.3\sqrt{\beta_p}, \qquad
      \beta_p = \frac{0.4\,n_e\,T_e\,a}{R_0\,q_{95}^2}$$

    Parameters
    ----------
    n_e : float
        Electron density in the unit the coefficient was fitted for [any].
        The original source does not record which.
    T_e_keV : float
        Electron temperature [keV].
    R0 : float
        Major radius [m].
    a : float
        Minor radius [m].
    q_95 : float
        Safety factor at the 95% flux surface [-].

    Returns
    -------
    float
        Bootstrap fraction of the plasma current [-].

    Physical interpretation
    -----------------------
    Neoclassical bootstrap current scales as $\epsilon^{1/2}\beta_p$; the inner
    expression is a crude $\beta_p$ estimate from an $n T$ pressure and the
    cylindrical $q$-$I_p$ relation, and the outer square root is a fit.

    Validity
    --------
    Empirical fit.  Labelled "ITER scaling for bootstrap current" in the original
    VAFT source without a publication or unit system for the coefficients 0.3
    and 0.4.  The physics it approximates is the $\sqrt{\epsilon}\,\beta_p$
    scaling of Peeters [1]_ and Wesson [2]_.

    Limitations
    -----------
    Unsourced coefficients and unstated density unit; the $\beta_p$ estimate
    ignores profile shape and the ion pressure.  Tracked in #361.

    References
    ----------
    .. [1] A. G. Peeters, Plasma Phys. Control. Fusion 42 (2000) B231, Sec. 2.
    .. [2] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 4.9 (bootstrap current).
    """
    beta_p = 0.4 * n_e * T_e_keV * a / (R0 * q_95**2)
    return 0.3 * np.sqrt(beta_p)

def bremsstrahlung_power_density_from_T_e_p_Z_eff(
    T_e: float,
    p: float,
    Z_eff: float = 2.0
) -> float:
    r"""Bremsstrahlung power density in the pressure form.

    $$S_B = Z_{\mathrm{eff}}\,K_B\,\frac{p_{\mathrm{bar}}^2}{T_{\mathrm{keV}}^{3/2}}
      \ [\mathrm{MW/m^3}], \qquad K_B = 0.052$$

    which is the NRL $P_{br} = 1.69\times10^{-38}\,Z_{\mathrm{eff}}n_e^2\sqrt{T_e}$
    rewritten with $n_e = p/(2T)$ (equal electron and ion pressure).

    Parameters
    ----------
    T_e : float
        Electron temperature [eV].
    p : float
        Total plasma pressure [Pa].
    Z_eff : float, optional
        Effective charge; default 2 [-].

    Returns
    -------
    float
        Radiated power density [W/m^3].

    Convention
    ----------
    Pressure is converted to bar ($10^5$ Pa) and temperature to keV internally;
    the prefactor 0.052 MW/m^3 follows exactly from the NRL coefficient and
    $p = 2n_eT$, so a pressure that already includes fast ions or unequal
    $T_i$ overestimates $n_e$.

    Assumptions
    -----------
    Maxwellian electrons, Gaunt factor 1, $n_i = n_e$ and $T_i = T_e$; no
    recombination or line radiation.

    References
    ----------
    .. [1] NRL Plasma Formulary (2019), p. 58 (bremsstrahlung).
    .. [2] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Ch. 4 (radiation losses).
    """
    # convert pressure to 1e5 Pa units
    p_bar = p / 1e5
    # convert temperature to keV
    T_k_keV = T_e * 1e-3

    S_B_MW_m3 = K_B_COEF * (p_bar ** 2) / (T_k_keV ** 1.5)

    return Z_eff * S_B_MW_m3 * 1e6  # W/m^3

def bremsstrahlung_power_density_from_n_e_T_e_Z_eff(
    n_e_m3: float,
    T_e_eV: float,
    *,
    Z_eff: float = 2.0
) -> float:
    r"""Maxwellian free-free (bremsstrahlung) power density from first principles.

    $$S_B = \frac{2^{1/2}}{3\pi^{5/2}}\,\frac{e^6}{\varepsilon_0^3c^3h\,m_e^{3/2}}\;
      Z_{\mathrm{eff}}\,n_e^2\,\sqrt{k_BT_e}$$

    Parameters
    ----------
    n_e_m3 : float
        Electron density [m^-3].
    T_e_eV : float
        Electron temperature, converted to joules internally [eV].
    Z_eff : float, optional
        Effective charge, keyword-only; default 2 [-].

    Returns
    -------
    float
        Radiated power density [W/m^3].

    Convention
    ----------
    The prefactor evaluates to $4.2\times10^{-29}$ W m^3 J^-1/2, equal to the NRL
    $1.69\times10^{-38}\,n_e^2\sqrt{T_{\mathrm{eV}}}$ used by
    :func:`bremsstrahlung_radiation_power_from_z_eff_n_e_t_e`; the two agree to
    0.2 %, which is the rounding in NRL's published coefficient.

    ``Z_eff`` is keyword-only, unlike its two siblings.  The NRL form takes
    $(Z_{\mathrm{eff}}, n_e, T_e)$ and this one takes $(n_e, T_e)$, so a caller
    moving between them by name used to swap the arguments silently and get a
    number wrong by 27 orders of magnitude; there is no argument order that can
    now do that without raising.  Tracked in #760.

    Assumptions
    -----------
    Maxwellian, non-relativistic electrons; Gaunt factor 1; $Z_{\mathrm{eff}}$
    absorbs $\sum_Z Z^2 n_Z / n_e$.

    Validity
    --------
    $T_e \ll m_ec^2$; relativistic corrections exceed 10 % above ~50 keV.

    See Also
    --------
    bremsstrahlung_radiation_power_from_z_eff_n_e_t_e
    bremsstrahlung_power_density_from_T_e_p_Z_eff

    References
    ----------
    .. [1] G. B. Rybicki and A. P. Lightman, *Radiative Processes in
           Astrophysics*, Wiley (1979), Eq. (5.15b).
    .. [2] I. H. Hutchinson, *Principles of Plasma Diagnostics*, 2nd ed.,
           Cambridge University Press (2002), Sec. 5.3.
    """

    T_J = T_e_eV * QE
    return C_B * Z_eff * n_e_m3**2 * np.sqrt(T_J)

def bremsstrahlung_radiation_power_from_z_eff_n_e_t_e(Z_eff: float,
                                        n_e: float,
                                        T_e_eV: float) -> float:
    r"""Bremsstrahlung power density, NRL engineering form.

    $$p_{br} = 1.69\times10^{-38}\,Z_{\mathrm{eff}}\,n_e^2\,\sqrt{T_e\,[\mathrm{eV}]}
      \ [\mathrm{W/m^3}]$$

    Parameters
    ----------
    Z_eff : float
        Effective charge [-].
    n_e : float
        Electron density [m^-3].
    T_e_eV : float
        Electron temperature [eV].

    Returns
    -------
    float
        Radiated power density [W/m^3].

    Convention
    ----------
    The NRL Formulary $1.69\times10^{-32}\,n_e T_e^{1/2}\sum Z^2n_Z$ W/cm^3 with
    cm^-3 densities, converted to SI, and $\sum Z^2 n_Z = Z_{\mathrm{eff}}n_e$.

    Assumptions
    -----------
    Maxwellian electrons, Gaunt factor 1.

    References
    ----------
    .. [1] NRL Plasma Formulary (2019), p. 58.
    """
    return 1.69e-38 * Z_eff * n_e**2 * np.sqrt(T_e_eV)

def calc_inverse_aspect_ratio(
    a_m: Union[float, np.ndarray],
    R_geo_m: Union[float, np.ndarray],
) -> Union[float, np.ndarray]:
    r"""Inverse aspect ratio $\varepsilon = a/R_{geo}$ with input validation.

    $$\varepsilon = \frac{a}{R_{geo}}$$

    Parameters
    ----------
    a_m : float or np.ndarray
        Minor radius, strictly positive [m].
    R_geo_m : float or np.ndarray
        Geometric major radius, strictly positive [m].

    Returns
    -------
    float or np.ndarray
        Inverse aspect ratio [-].

    Raises
    ------
    ValueError
        For non-finite or non-positive input.

    See Also
    --------
    inverse_aspect_ratio_from_a_R : the unvalidated form.
    """
    a_arr = _validate_positive("a_m", a_m)
    r_arr = _validate_positive("R_geo_m", R_geo_m)
    epsilon = inverse_aspect_ratio_from_a_R(a_arr, r_arr)
    return _maybe_scalar(epsilon)

def check_kadomtsev_constraint(
    alpha_I: Union[float, np.ndarray],
    alpha_B: Union[float, np.ndarray],
    alpha_P: Union[float, np.ndarray],
    alpha_n: Union[float, np.ndarray],
    alpha_R: Union[float, np.ndarray],
    tol: float = 1e-6,
) -> Union[bool, np.ndarray]:
    r"""Whether engineering exponents satisfy the Kadomtsev constraint within a tolerance.

    $$|\alpha_K| \le \mathrm{tol}$$

    with $\alpha_K$ from :func:`kadomtsev_constraint_from_engineering_exponents`.

    Parameters
    ----------
    alpha_I : float or np.ndarray
        Exponent of the plasma current [-].
    alpha_B : float or np.ndarray
        Exponent of the toroidal field [-].
    alpha_P : float or np.ndarray
        Exponent of the heating power [-].
    alpha_n : float or np.ndarray
        Exponent of the density [-].
    alpha_R : float or np.ndarray
        Exponent of the major radius [-].
    tol : float, optional
        Absolute tolerance on $|\alpha_K|$; default 1e-6 [-].

    Returns
    -------
    bool or np.ndarray
        ``True`` where the constraint holds [bool].

    Limitations
    -----------
    The default tolerance is far tighter than any published fit satisfies
    (ITER89P misses by 0.15); pass a physically motivated ``tol``.
    :func:`verify_kadomtsev_constraint` returns the same residual from the
    dimensionless indices (#351).

    References
    ----------
    .. [1] B. B. Kadomtsev, Sov. J. Plasma Phys. 1 (1975) 295.
    """
    tol_arr = _validate_positive("tol", tol)
    alpha_k = np.asarray(
        kadomtsev_constraint_from_engineering_exponents(
            alpha_I=alpha_I,
            alpha_B=alpha_B,
            alpha_P=alpha_P,
            alpha_n=alpha_n,
            alpha_R=alpha_R,
        ),
        dtype=float,
    )
    satisfied = np.abs(alpha_k) <= tol_arr
    if np.asarray(satisfied).ndim == 0:
        return bool(satisfied)
    return satisfied

def confinement_factor_ITER89P(tau_E_exp: float, tau_E_ITER89P: float) -> float:

    r"""Confinement enhancement factor $H_{89}$ relative to ITER89P.

    $$H_{89} = \frac{\tau_{E,\mathrm{exp}}}{\tau_{E,\mathrm{ITER89P}}}$$

    Parameters
    ----------
    tau_E_exp : float
        Measured energy confinement time [s].
    tau_E_ITER89P : float
        ITER89P prediction for the same discharge [s].

    Returns
    -------
    float
        $H_{89}$; values above 1 beat the L-mode scaling [-].

    References
    ----------
    .. [1] P. N. Yushmanov et al., Nucl. Fusion 30 (1990) 1999.
    """

    return tau_E_exp / tau_E_ITER89P

def confinement_scaling_basis(scaling: str) -> ConfinementScalingBasis:
    r"""Energy and power basis a published confinement scaling was fitted on.

    $$\tau_{E,th} = W_{th}/P,\qquad \tau_{E,global} = W/P,\qquad W = W_{th} + W_{fast}$$

    Parameters
    ----------
    scaling : str
        Scaling name: a key of ``_SCALING_COEFS`` (``"ITER89P"``, ``"ITER97L"``,
        ``"H98y2"``, ``"NSTX2006H"``, ``"NSTX2006L"``, ``"Kurskiev2022"``,
        ``"ITPA20"``, ``"ITPA20_IL"``) or
        one of the ohmic/L-mode forms ``"NeoAlcator"``, ``"Goldston84L"``,
        ``"Goldston84OhmicL"`` [str].

    Returns
    -------
    ConfinementScalingBasis
        ``energy_basis`` (``"thermal"``, ``"global"`` or ``"unaudited"``),
        ``power_basis`` (``"p_loss"``, ``"p_abs"``, ``"p_heat"``, ``"none"`` or
        ``"unaudited"``) and the source of each assignment [-].

    Raises
    ------
    KeyError
        A scaling with no declared basis.

    Convention
    ----------
    An H factor is the conventional one only when the observed confinement
    time has the scaling's energy basis: a thermal scaling against
    $\tau_{E,th}$, a global one against $\tau_{E,global}$.  ``"unaudited"``
    means the original paper has not been checked for that definition; the
    value is never guessed, and a caller must treat it as unknown.

    References
    ----------
    .. [1] ITER Physics Expert Groups, Nucl. Fusion 39 (1999) 2175, Ch. 2, Sec. 6.
    .. [2] S. M. Kaye et al., Nucl. Fusion 46 (2006) 848.
    """
    if scaling not in _SCALING_BASES:
        raise KeyError(f"no energy/power basis declared for scaling {scaling!r}; known: {sorted(_SCALING_BASES)}")
    entry = _SCALING_BASES[scaling]
    return ConfinementScalingBasis(scaling, entry["energy_basis"], entry["power_basis"],
                                   entry["energy_source"], entry["power_source"])

@_jacobian(_confinement_time_jacobian, wrt=("P_loss", "W_th"),
           domain=lambda P_loss, W_th, **_: P_loss > 0 and W_th >= 0)
def confinement_time_from_P_loss_W_th(P_loss: float, W_th: float) -> float:
    r"""Energy confinement time as stored energy over loss power.

    $$\tau_E = \frac{W_{th}}{P_{loss}}$$

    Parameters
    ----------
    P_loss : float
        Loss power [W].
    W_th : float
        Thermal stored energy [J].

    Returns
    -------
    float
        Energy confinement time [s].

    Convention
    ----------
    Thermal energy over *loss* power ($P_{heat} - dW/dt$, with or without
    radiation subtracted depending on the database); the ITER definition of
    $\tau_{E,th}$ needs $P_{loss}$ from :func:`loss_power_from_p_heat_dWdt_p_rad`.

    Semantics
    ---------
    produces: energy_confinement_time

    Uncertainty propagation
    -----------------------
    Inputs in order $(P_{loss}, W_{th})$ [W, J]; analytic Jacobian
    $J = (-W_{th}/P_{loss}^2,\ 1/P_{loss})$ [s/W, s/J], so
    $(\sigma_\tau/\tau)^2 = (\sigma_P/P)^2 + (\sigma_W/W)^2 - 2\,\mathrm{cov}(P, W)/(PW)$.
    First order is local: the ratio is nonlinear in $P_{loss}$, and the
    linearization stops describing the spread as $\sigma_P/P_{loss}$ grows
    (it is singular at $P_{loss} = 0$ and undefined for $P_{loss} \le 0$, a
    transient with $dW/dt$ larger than the heating power; the declared domain is
    $P_{loss} > 0$, $W_{th} \ge 0$, so such a point or Monte Carlo draw is refused). A positive
    correlation (both inputs from one equilibrium reconstruction, or $dW/dt$
    in $P_{loss}$ from the same $W_{th}$) *reduces* the spread. Measured-input
    covariance only; no confinement-scaling or model-form uncertainty.

    References
    ----------
    .. [1] ITER Physics Expert Groups, Nucl. Fusion 39 (1999) 2175, Ch. 2, Sec. 3.
    """
    return W_th / P_loss

def coulomb_logarithm_from_n_T(
    n_m3: Union[float, np.ndarray],
    T_eV: Union[float, np.ndarray],
) -> Union[float, np.ndarray]:
    r"""Coulomb logarithm $\ln\Lambda$ for electron collisions above 10 eV.

    $$\ln\Lambda = 30.9 - \ln\!\left(\frac{\sqrt{n_e\,[\mathrm{m^{-3}}]}}{T_e\,[\mathrm{eV}]}\right)$$

    Parameters
    ----------
    n_m3 : float or np.ndarray
        Electron density, strictly positive [m^-3].
    T_eV : float or np.ndarray
        Electron temperature, strictly positive [eV].

    Returns
    -------
    float or np.ndarray
        Coulomb logarithm [-].

    Raises
    ------
    ValueError
        For non-finite or non-positive input.

    Convention
    ----------
    The NRL Formulary electron-electron/electron-ion form $24 - \ln(n_e^{1/2}
    T_e^{-1})$ with $n_e$ in cm^-3, converted to m^-3 ($24 + \ln 10^3 = 30.9$);
    also the convention of the Verdoolaege confinement-database analysis.

    Validity
    --------
    $T_e > 10$ eV; below that the NRL low-temperature branch applies.

    References
    ----------
    .. [1] NRL Plasma Formulary (2019), p. 34 (Coulomb logarithm).
    .. [2] G. Verdoolaege et al., Nucl. Fusion 61 (2021) 076006, Sec. 2.
    """
    n_arr = _validate_positive("n_m3", n_m3)
    t_arr = _validate_positive("T_eV", T_eV)
    ln_lambda = 30.9 - np.log(np.sqrt(n_arr) / t_arr)
    return _maybe_scalar(ln_lambda)

def current_density_from_B(B: Union[float, np.ndarray],
                          R: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
    r"""Toroidal current density from the radial derivative of a poloidal field.

    $$j_\varphi \approx \frac{1}{\mu_0}\,\frac{dB}{dR}$$

    the slab form of Ampere's law $\mu_0 j_\varphi = \partial B_Z/\partial R -
    \partial B_R/\partial Z$ with the second term dropped.

    Parameters
    ----------
    B : float or np.ndarray
        Poloidal field component (normally $B_Z$) along a 1-D radial cut [T].
    R : float or np.ndarray
        Major radius of the samples, monotonic [m].

    Returns
    -------
    float or np.ndarray
        Current density along the cut [A/m^2].

    Convention
    ----------
    Sign follows the COCOS $\sigma_{B_p}$ of the supplied component: with the
    midplane $B_Z$ of a standard equilibrium the result has the sign of the
    plasma current.  The $\partial B_R/\partial Z$ contribution is neglected, so
    the answer is exact only on the midplane of an up-down symmetric equilibrium.

    Numerical notes
    ---------------
    ``numpy.gradient`` along the single supplied axis (second-order interior,
    first-order ends); pass a 1-D slice, not a 2-D map.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 3.1 (Ampere's law in the tokamak).
    """
    return gradient(R, B) / MU0

def current_drive_efficiency(n_e: float,
                           T_e_keV: float,
                           Z_eff: float = 1.0) -> float:
    r"""Heuristic lower-hybrid current-drive efficiency $\eta_{CD}$.

    $$\eta_{CD} = 0.3\,\sqrt{\frac{n_e\,T_e}{Z_{\mathrm{eff}}}}$$

    Parameters
    ----------
    n_e : float
        Electron density in the unit the coefficient was fitted for [any].
        The original source does not record which.
    T_e_keV : float
        Electron temperature [keV].
    Z_eff : float, optional
        Effective charge; default 1 [-].

    Returns
    -------
    float
        Efficiency figure in the coefficient's own normalisation [-].

    Physical interpretation
    -----------------------
    Lower-hybrid current drive becomes more efficient at higher temperature and
    lower $Z_{\mathrm{eff}}$ because the wave-driven fast electrons slow down on
    a hotter, cleaner background; the density factor here is unusual (the
    standard figure of merit $\eta = n_e I_{CD} R / P$ *divides* by density).

    Validity
    --------
    Empirical fit.  Labelled "ITER scaling for lower hybrid current drive" in the
    original VAFT source; no publication, dataset or unit system for the
    coefficient 0.3 was recorded, so the number is a placeholder.  The theory of
    the efficiency and its $T_e/Z_{\mathrm{eff}}$ dependence is Fisch [1]_.

    Limitations
    -----------
    Unsourced coefficient and unstated density unit; treat the result as
    qualitative.  Tracked in #361.

    References
    ----------
    .. [1] N. J. Fisch, Rev. Mod. Phys. 59 (1987) 175, Sec. VI (lower-hybrid
           current-drive efficiency).
    """
    return 0.3 * (n_e * T_e_keV / Z_eff)**0.5

def current_limit_from_beta(beta_N: float,
                          a: float,
                          B0: float) -> float:
    r"""Current figure obtained by substituting $\beta_N$ for $q$ in the cylindrical relation.

    $$I_p = \frac{2\pi a^2 B_0}{\mu_0\,\beta_N}$$

    Parameters
    ----------
    beta_N : float
        Normalised beta in whatever convention the caller uses [-].
    a : float
        Minor radius [m].
    B0 : float
        Toroidal field on axis [T].

    Returns
    -------
    float
        The expression above [A m].

    Limitations
    -----------
    Dimensionally the cylindrical-$q$ formula with $\beta_N$ in the place of $q$;
    no derivation or source records what this is meant to bound, and the result
    depends on the units chosen for $\beta_N$.  Kept for compatibility; prefer
    :func:`vaft.formula.stability.beta_N_from_beta_a_B0_Ip` and the Troyon limit.
    Tracked in #362.

    References
    ----------
    .. [1] F. Troyon et al., Plasma Phys. Control. Fusion 26 (1984) 209 (the
           $\beta_N$ limit this appears to invert).
    """
    return 2 * np.pi * a**2 * B0 / (MU0 * beta_N)

def current_limit_from_q(q_95: float,
                        a: float,
                        B0: float) -> float:
    r"""Plasma current at a prescribed edge safety factor, cylindrical approximation.

    $$I_p = \frac{2\pi a^2 B_0}{\mu_0\,q_{95}}$$

    the inversion of $q_{cyl} = 2\pi a^2 B_0/(\mu_0 R\,I_p)\times R/R$ for a
    circular cross-section written per unit major radius.

    Parameters
    ----------
    q_95 : float
        Target safety factor at the 95% flux surface [-].
    a : float
        Minor radius [m].
    B0 : float
        Toroidal field on axis [T].

    Returns
    -------
    float
        Plasma current reaching ``q_95`` [A m].

    Assumptions
    -----------
    Circular, large-aspect-ratio cylinder with $q_{95}$ standing in for the
    cylindrical $q_a$.

    Limitations
    -----------
    As written the expression lacks the $1/R$ of the cylindrical safety factor
    $q = 2\pi a^2 B_0/(\mu_0 R I_p)$ and therefore returns $I_p R$ [A m], not a
    current; divide by $R$ for amperes.  Elongation is ignored.  Tracked in #362.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 3.4 (cylindrical safety factor).
    """
    return 2 * np.pi * a**2 * B0 / (MU0 * q_95)

def cyclotron_synchrotron_power_density_scaling_from_n_e_B_t_T_e(
    n_e_m3: float,
    B_t_T: float,
    T_e_eV: float,
) -> float:
    r"""Classical electron-cyclotron emission power density, no reabsorption.

    $$p_{\mathrm{cyc}} = \frac{e^4}{3\pi\varepsilon_0 m_e^3c^3}\,n_e\,B_t^2\,k_BT_e
      \approx 6.2\times10^{-17}\,B_t^2\,n_e\,T_{\mathrm{keV}}\ \mathrm{W/m^3}$$

    Parameters
    ----------
    n_e_m3 : float
        Electron density [m^-3].
    B_t_T : float
        Magnetic field [T].
    T_e_eV : float
        Electron temperature, converted to joules internally [eV].

    Returns
    -------
    float
        Emitted cyclotron power density [W/m^3].

    Physical interpretation
    -----------------------
    Total single-particle cyclotron radiation of a non-relativistic Maxwellian;
    the plasma is optically thick at the low harmonics, so the net loss is a
    small fraction of this, set by wall reflectivity and $\beta$ (Trubnikov).

    Validity
    --------
    Non-relativistic electrons; an upper bound on the loss, intended as a
    start-up loss-channel estimate rather than a radiation-transport result.

    Limitations
    -----------
    Ignores reabsorption and wall reflection, which reduce the net loss by one
    to two orders of magnitude in a tokamak.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Ch. 4 (cyclotron radiation).
    .. [2] B. A. Trubnikov, in *Reviews of Plasma Physics*, Vol. 7, Consultants
           Bureau (1979), p. 345.
    """
    # QE/EPS0/ME/C_LIGHT, not the lower-case `e`/`epsilon_0`/`m_e`/`c` this
    # module used to define: #368 moved those to `constants.py` because the
    # lower-case spellings leaked into `vaft.formula.__all__`, and this function
    # was the one caller left behind, so it has raised NameError ever since
    # (#753).
    coeff = QE**4 / (3.0 * np.pi * EPS0 * ME**3 * C_LIGHT**3)
    return coeff * n_e_m3 * B_t_T**2 * (T_e_eV * QE)

def cylindrical_safety_factor_from_R_B_epsilon_I_f_kappa_delta(R: float,
                                           B: float,
                                           epsilon: float,
                                           I: float,
                                           f_kappa_delta: float) -> float:
    r"""Cylindrical safety factor with a shape function.

    $$q_{cyl} = \frac{2\pi\,\varepsilon^2 R\,B_T}{\mu_0\,I_p\,f(\kappa,\delta)}
      = \frac{2\pi a^2 B_T}{\mu_0 R I_p}\,\frac{1}{f(\kappa,\delta)}$$

    Parameters
    ----------
    R : float
        Major radius [m].
    B : float
        Toroidal field [T].
    epsilon : float
        Inverse aspect ratio $a/R$ [-].
    I : float
        Plasma current [A].
    f_kappa_delta : float
        Shape function; $1/\kappa$ reproduces the ITER $q_{cyl} = 5a^2\kappa B/(RI_{MA})$ [-].

    Returns
    -------
    float
        Cylindrical safety factor [-].

    Convention
    ----------
    $2\pi/\mu_0 = 5\times10^6$, so with SI current this is the ITER Physics
    Basis $q_{cyl}$ once $f = 1/\kappa$; the shape function is left to the
    caller because the ITER-89P and IPB98 databases used different $\kappa$
    definitions.  Unlike the flux-derivative $q$ it never changes sign.

    References
    ----------
    .. [1] ITER Physics Expert Groups, Nucl. Fusion 39 (1999) 2175, Ch. 2,
           Sec. 3 (definition of $q_{cyl}$).
    .. [2] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011), Sec. 3.4.
    """
    return (2.0 * np.pi * epsilon**2 * R * B) / (MU0 * I * f_kappa_delta)

def decay_index_from_bz(
    r: np.ndarray, b_z: np.ndarray, *, axis: int = -1
) -> np.ndarray:
    r"""Field decay index $n = -\dfrac{R}{B_Z}\dfrac{\partial B_Z}{\partial R}$.

    How fast the vertical field falls off with major radius, which decides
    whether the radial force balance holding a current ring is *stable*.

    Parameters
    ----------
    r : array_like
        Major radius, monotonic, along ``axis`` [m].
    b_z : array_like
        Vertical field sampled on ``r``; may carry extra leading axes [T].
    axis : int, optional
        Axis of ``b_z`` along which ``r`` varies [-].

    Returns
    -------
    np.ndarray
        Decay index, ``nan`` where $B_Z$ vanishes [-].

    Convention
    ----------
    $0 < n < 1.5$ is the passively stable window of a rigid current ring:
    below zero the ring is vertically unstable (the field lines curve the
    wrong way to restore a vertical displacement), above 1.5 it is radially
    unstable (the field falls off too fast to restore a radial one).

    Limitations
    -----------
    Returns ``nan`` where $B_Z$ crosses zero rather than a large number.  The
    index is genuinely undefined on that surface, and a spike there is an
    artefact of the division, not a physical instability -- which is what a
    reader would otherwise take from it.

    Validity
    --------
    Machine-independent.  The stable window quoted above is the rigid-ring
    result: it assumes a thin current ring and no conducting wall, so a real
    vessel widens it.

    References
    ----------
    .. [1] V. S. Mukhovatov and V. D. Shafranov, Nucl. Fusion 11 (1971) 605,
           Sec. 2 (equilibrium of a current ring in a vertical field).
    .. [2] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 3.7 (vertical field and positional stability).
    """
    r_arr = np.asarray(r, dtype=float)
    b_arr = np.asarray(b_z, dtype=float)
    gradient_bz = np.gradient(b_arr, r_arr, axis=axis, edge_order=2)
    shape = [1] * b_arr.ndim
    shape[axis] = r_arr.size
    with np.errstate(divide="ignore", invalid="ignore"):
        index = -r_arr.reshape(shape) * gradient_bz / b_arr
    return np.where(np.isfinite(index), index, np.nan)

def eK_from_K(K: float) -> float:
    r"""Elongation parameter $e_K$ of the virial relations.

    $$e_K = \frac{K^2 - 1}{K^2 + 1}$$

    Parameters
    ----------
    K : float
        Elongation $\kappa$ [-].

    Returns
    -------
    float
        $e_K$, 0 for a circle and $\to1$ for infinite elongation [-].

    Physical interpretation
    -----------------------
    The ellipticity measure in which the Martynov-Pustovitov virial
    approximations are linear.

    References
    ----------
    .. [1] A. A. Martynov and V. D. Pustovitov, Phys. Plasmas 31 (2024), "Virial
           relations for elongated plasmas in tokamaks", definition preceding
           Eq. (21).
    """
    return (K**2 - 1) / (K**2 + 1)

def ec_heating_power_from_I_ec_V_ec(I_ec: float,
                                    V_ec: float) -> float:
    r"""Electron-cyclotron launched power $P_{ec} = I_{ec}V_{ec}$.

    $$P_{\mathrm{ec}} = I_{\mathrm{ec}}\,V_{\mathrm{ec}}$$

    Parameters
    ----------
    I_ec : float
        Gyrotron beam current [A].
    V_ec : float
        Gyrotron beam voltage [V].

    Returns
    -------
    float
        Electrical beam power of the source [W].

    Limitations
    -----------
    Gyrotron electrical power, not RF power (efficiency ~30-50 %) and not the
    power absorbed by the plasma.
    """
    return I_ec * V_ec

def electron_pressure(n_e: Union[float, np.ndarray],
                      T_e: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
    r"""Electron thermal pressure from density and temperature, $p_e = n_e T_e$.

    $$p_e = n_e\,k_B T_e = n_e\,T_e[\mathrm{eV}]\;e$$

    Parameters
    ----------
    n_e : float or np.ndarray
        Electron density, finite and non-negative [m^-3].
    T_e : float or np.ndarray
        Electron temperature, finite and non-negative, broadcastable
        against ``n_e`` [eV].

    Returns
    -------
    float or np.ndarray
        Electron pressure; $10^{19}\,\mathrm{m^{-3}}$ at 100 eV is 160.2 Pa [Pa].

    Raises
    ------
    ValueError
        A non-finite or negative density or temperature.

    Convention
    ----------
    Temperature in electronvolts, converted with the elementary charge
    :data:`vaft.formula.constants.QE`; density in m^-3, not $10^{19}$ m^-3.
    This is one species' pressure: the total kinetic pressure is
    $p_e + \sum_i p_i$, which is **not** $2 n_e T_e$ unless $T_i = T_e$ and
    $\sum_i n_i = n_e$ (hydrogen, no impurity).

    Assumptions
    -----------
    Isotropic Maxwellian electrons; a fast or anisotropic population needs
    its own $p_\parallel$, $p_\perp$.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Ch. 2 (the ideal-gas pressure of a plasma species).

    See Also
    --------
    ion_pressure
    """
    return _thermal_pressure(n_e, T_e, n_name="n_e", T_name="T_e")

def elongation_from_RZ_boundary(R: np.ndarray,
                               Z: np.ndarray) -> float:
    r"""Boundary elongation $\kappa$ from the extremal points of a contour.

    $$\kappa = \frac{Z_{\max} - Z_{\min}}{2a}, \qquad a = \frac{R_{\max} - R_{\min}}{2}$$

    Parameters
    ----------
    R : np.ndarray
        Major radius of the boundary points [m].
    Z : np.ndarray
        Height of the boundary points [m].

    Returns
    -------
    float
        Elongation [-].

    Convention
    ----------
    The IMAS ``boundary.elongation`` definition (half-height over half-width of
    the bounding box), not an area-based or flux-surface-averaged elongation.

    Limitations
    -----------
    A single-point or degenerate contour divides by zero; no validation.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 3.1.
    .. [2] IMAS Data Dictionary, ``equilibrium.time_slice[:].boundary.elongation``.
    """
    a = (R.max() - R.min()) / 2
    return (Z.max() - Z.min()) / (2 * a)

def engineering_exponents_from_dimensionless_coeffs(
    mu_rho: float,
    mu_beta: float,
    mu_nu: float,
    a_P: float,
    mu_q: float = 0.0,
) -> Tuple[float, float, float, float, float]:
    r"""Engineering exponents $(\alpha_I, \alpha_B, \alpha_P, \alpha_n, \alpha_R)$ from dimensionless indices.

    $$\alpha_I = -D\mu_q, \quad
      \alpha_B = D(\mu_q - \mu_\rho - 2\mu_\beta - 1), \quad
      \alpha_n = D(\mu_\beta + \mu_\nu) - \alpha_P, \quad
      \alpha_R = D(\mu_q + \mu_\nu - \mu_\rho) - 3\alpha_P$$

    with $D = 1 + \alpha_P$, for
    $\Omega_i\tau_E \propto \rho_*^{\mu_\rho}\beta^{\mu_\beta}\nu_*^{\mu_\nu}q^{\mu_q}$.

    Parameters
    ----------
    mu_rho : float
        Gyroradius index [-].
    mu_beta : float
        Beta index [-].
    mu_nu : float
        Collisionality index [-].
    a_P : float
        Engineering exponent of the heating power [-].
    mu_q : float, optional
        Safety-factor index; the default 0 returns the fixed-$q$ form, in which
        the current exponent is folded into ``alpha_B`` and ``alpha_R`` [-].

    Returns
    -------
    alpha_I : float
        Exponent of the plasma current [-].
    alpha_B : float
        Exponent of the toroidal field [-].
    alpha_P : float
        Exponent of the heating power, equal to ``a_P`` [-].
    alpha_n : float
        Exponent of the density [-].
    alpha_R : float
        Exponent of the size at fixed aspect ratio [-].

    Raises
    ------
    ValueError
        For non-finite input, or when $1 + \alpha_P$ vanishes [-].

    Assumptions
    -----------
    The same closure as
    :func:`dimensionless_scaling_coeffs_from_engineering_scaling_coeffs`, run
    backwards on the $n$, $B$, $R$ and $I$ exponents.  A ``a_P`` that does not
    match the indices' own temperature exponent gives exponents that violate
    the Kadomtsev constraint by exactly that mismatch.

    References
    ----------
    .. [1] T. C. Luce, C. C. Petty and J. G. Cordey, Plasma Phys. Control.
           Fusion 50 (2008) 043001, Sec. 3.
    """
    values = np.asarray([mu_rho, mu_beta, mu_nu, a_P, mu_q], dtype=float)
    if np.any(~np.isfinite(values)):
        raise ValueError(f"indices and a_P must be finite. Got {values!r}")
    denom = 1.0 + float(a_P)
    if abs(denom) < 1e-9:
        raise ValueError(
            "a_P = -1 leaves 1 + a_P = 0, so the indices fix no finite "
            f"engineering exponent. Got a_P={a_P!r}."
        )
    alpha_i = -denom * mu_q
    alpha_b = denom * (mu_q - mu_rho - 2.0 * mu_beta - 1.0)
    alpha_n = denom * (mu_beta + mu_nu) - a_P
    alpha_r = denom * (mu_q + mu_nu - mu_rho) - 3.0 * a_P
    return alpha_i, alpha_b, float(a_P), alpha_n, alpha_r

def exact_volume_from_RZ_contour(R: np.ndarray,
                                 Z: np.ndarray) -> float:
    r"""Plasma volume from a closed $(R, Z)$ contour by Green's theorem.

    $$V = \pi\oint R^2\,dZ$$

    exact for the solid of revolution swept by a closed poloidal contour, with
    no mean-radius approximation.

    Parameters
    ----------
    R : np.ndarray
        Major radius of the contour vertices, at least 3 [m].
    Z : np.ndarray
        Height of the contour vertices [m].

    Returns
    -------
    float
        Volume of the solid of revolution [m^3].

    Convention
    ----------
    Orientation-agnostic: the sign of the contour integral follows the traversal
    direction and the magnitude is returned.  The contour is closed automatically
    when its first and last points differ.

    Limitations
    -----------
    Unlike :func:`volume_from_RZ_boundary`, which factors the integral as
    $2\pi A_{\mathrm{poly}}\bar R$ with $\bar R = \mathrm{mean}(R)$, this
    evaluates the contour integral itself.  On VEST flux surfaces the two differ
    by up to ~6 % at the plasma edge, where the $\bar R$ factorisation is
    weakest.  Use this one whenever the volume is a reported quantity rather than
    an intermediate.

    Numerical notes
    ---------------
    Trapezoidal rule on $R^2$ between successive vertices, second-order in the
    vertex spacing.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 3.1.
    """
    R = np.asarray(R, dtype=float).reshape(-1)
    Z = np.asarray(Z, dtype=float).reshape(-1)
    if R.size != Z.size:
        raise ValueError("R and Z must have the same length")
    if R.size < 3:
        raise ValueError("a contour needs at least 3 points")
    if R[0] != R[-1] or Z[0] != Z[-1]:
        R = np.append(R, R[0])
        Z = np.append(Z, Z[0])
    # Trapezoidal ∮ R² dZ; the sign follows the traversal direction, so take
    # the magnitude and let the caller stay orientation-agnostic.
    return float(abs(np.pi * np.sum(0.5 * (R[:-1] ** 2 + R[1:] ** 2) * np.diff(Z))))

def flux_perturbation_from_normal_displacement(xi_n, grad_psi):
    r"""Ideal (flux-frozen) perturbed flux of a displacement normal to the flux surfaces.

    $$\delta\psi = -\boldsymbol\xi\cdot\nabla\psi_0 = -\xi_n\,|\nabla\psi_0|$$

    Parameters
    ----------
    xi_n : float or np.ndarray
        Displacement along $\hat{\mathbf n} = \nabla\psi_0/|\nabla\psi_0|$ [m].
    grad_psi : float or np.ndarray
        $|\nabla\psi_0|$ of the equilibrium flux, non-negative, in the unit of $\psi_0$ per metre [Wb/m or Wb/(rad m)].

    Returns
    -------
    float or np.ndarray
        Eulerian flux perturbation $\delta\psi$, in the unit of $\psi_0$ [Wb or Wb/rad].

    Raises
    ------
    ValueError
        ``grad_psi`` is negative.

    Convention
    ----------
    $\hat{\mathbf n}$ points up the gradient of $\psi_0$, so a positive
    $\xi_n$ moves a surface towards larger $\psi_0$ -- outward when $\psi$
    increases from the axis (COCOS 11 with positive $I_p$, VAFT's usual ODS
    storage), inward when it decreases (COCOS 11 with negative $I_p$, or a
    per-radian COCOS-3 flux with positive current). Only the
    normal component enters; a tangential displacement moves the surface
    into itself.

    Physical interpretation
    -----------------------
    Ideal MHD freezes the flux into the fluid, so a surface displaced by
    $\xi_n$ carries its $\psi_0$ with it and the flux at a fixed point
    changes by the linearised amount above. It is what turns a displacement
    into a perturbed flux map, and back: an eigenfunction solver's
    $\xi_n$ and its $\delta\psi$ are the same information.

    Assumptions
    -----------
    Linear, $|\xi_n\nabla\ln|\nabla\psi_0|| \ll 1$; ideal (no reconnection),
    so it fails in a resistive layer at a rational surface.

    References
    ----------
    .. [1] J. P. Freidberg, *Ideal MHD*, Cambridge University Press (2014),
           Sec. 8.3.
    """
    grad_psi = np.asarray(grad_psi, dtype=float)
    if np.any(grad_psi < 0.0):
        raise ValueError("grad_psi must be non-negative (it is |grad psi|)")
    return -np.asarray(xi_n, dtype=float) * grad_psi

def generalized_parabolic_profile(psi_n, *, core_value=1.0, edge_value=0.0, alpha=1.0, beta=1.0):
    r"""Smooth core profile that falls from its axis value to its edge value.

    $$f(\psi_N) = f_\mathrm{edge} + (f_\mathrm{core} - f_\mathrm{edge})\,(1 - \psi_N^{\alpha})^{\beta}$$

    Parameters
    ----------
    psi_n : float or np.ndarray
        Normalized poloidal flux, 0 on the magnetic axis and 1 on the boundary [-].
    core_value : float
        Value on the magnetic axis [any].
    edge_value : float
        Value on the boundary [any].
    alpha : float
        Radial exponent, positive; larger values flatten the core [-].
    beta : float
        Peaking exponent, positive; larger values peak the profile on axis [-].

    Returns
    -------
    f : float or np.ndarray
        Profile value, in the unit of *core_value* and *edge_value*, with the
        shape of *psi_n* [any].

    Raises
    ------
    ValueError
        *psi_n* is not finite or lies outside ``[0, 1]``, a parameter is not
        finite, or *alpha* or *beta* is not positive.

    Convention
    ----------
    The radial coordinate is the normalized poloidal flux
    $\psi_N = (\psi - \psi_\mathrm{axis})/(\psi_\mathrm{boundary} - \psi_\mathrm{axis})$,
    so the shape does not depend on the flux sign or on whether the flux is
    stored in Wb or Wb/rad. ``alpha = beta = 1`` is linear in $\psi_N$.

    Physical interpretation
    -----------------------
    The usual smooth L-mode or core shape for pressure, density, temperature
    or a prescribed current shape. It carries no pedestal: its gradient near
    the edge is set by *beta* alone.

    Assumptions
    -----------
    A monotone profile between two prescribed end values; nothing here
    solves for them.

    Limitations
    -----------
    The derivative is singular where an exponent is below one: at the axis
    for ``alpha < 1`` and at the boundary for ``beta < 1`` (see
    :func:`generalized_parabolic_profile_derivative`). Outside ``[0, 1]`` the
    form is undefined for non-integer exponents, so it is refused rather than
    extrapolated; a rounding excess below 1e-12 is clipped.
    """
    x = _profile_psi_n(psi_n, bounded=True)
    _finite_parameters(core_value=core_value, edge_value=edge_value, alpha=alpha, beta=beta)
    if alpha <= 0.0 or beta <= 0.0:
        raise ValueError("alpha and beta must be positive")
    return edge_value + (core_value - edge_value) * np.power(1.0 - np.power(x, alpha), beta)

def generalized_parabolic_profile_derivative(psi_n, *, core_value=1.0, edge_value=0.0, alpha=1.0, beta=1.0):
    r"""Gradient against normalized flux of :func:`generalized_parabolic_profile`.

    $$\frac{df}{d\psi_N} = -\alpha\beta\,(f_\mathrm{core} - f_\mathrm{edge})\,\psi_N^{\alpha-1}(1 - \psi_N^{\alpha})^{\beta-1}$$

    Parameters
    ----------
    psi_n : float or np.ndarray
        Normalized poloidal flux, in ``[0, 1]`` [-].
    core_value : float
        Value on the magnetic axis [any].
    edge_value : float
        Value on the boundary [any].
    alpha : float
        Radial exponent, positive [-].
    beta : float
        Peaking exponent, positive [-].

    Returns
    -------
    df_dpsi_n : float or np.ndarray
        Derivative with respect to $\psi_N$, in the profile's unit; infinite
        at an endpoint where the form is singular [any].

    Raises
    ------
    ValueError
        As for :func:`generalized_parabolic_profile`.

    Convention
    ----------
    This is $df/d\psi_N$, not $df/d\psi$. The physical gradient is
    $df/d\psi = (df/d\psi_N)/(\psi_\mathrm{boundary} - \psi_\mathrm{axis})$
    in whatever flux unit and sign the equilibrium carries; that conversion
    is left to the equilibrium layer so there is one owner of it.

    Physical interpretation
    -----------------------
    With a pressure profile this is the shape of $p'$ up to the constant flux
    span, which is what a Grad-Shafranov source term needs.

    Limitations
    -----------
    At ``psi_n = 0`` the value is $-\beta(f_\mathrm{core}-f_\mathrm{edge})$
    for ``alpha = 1``, zero for ``alpha > 1`` and infinite for
    ``alpha < 1``; at ``psi_n = 1`` it is zero for ``beta > 1``,
    $-\alpha(f_\mathrm{core}-f_\mathrm{edge})$ for ``beta = 1`` and infinite
    for ``beta < 1``. A flat profile (``core_value == edge_value``) has a
    zero gradient everywhere, singular exponents included.

    """
    x = _profile_psi_n(psi_n, bounded=True)
    _finite_parameters(core_value=core_value, edge_value=edge_value, alpha=alpha, beta=beta)
    if alpha <= 0.0 or beta <= 0.0:
        raise ValueError("alpha and beta must be positive")
    amplitude = core_value - edge_value
    if amplitude == 0.0:
        # a flat profile, even where the shape is singular; a scalar psi_n gets
        # the same np.float64 the sloped path returns, not a 0-d ndarray
        return np.float64(0.0) if x.ndim == 0 else np.zeros_like(x)
    with np.errstate(divide="ignore", invalid="ignore"):
        inner = np.power(x, alpha - 1.0) if alpha != 1.0 else np.ones_like(x)
        outer = np.power(1.0 - np.power(x, alpha), beta - 1.0) if beta != 1.0 else np.ones_like(x)
        return -alpha * beta * amplitude * inner * outer

def goldston_l_mode_confinement_time_from_I_P_R_a_kappa(
    I_p: Union[float, np.ndarray],
    P: Union[float, np.ndarray],
    R: Union[float, np.ndarray],
    a: Union[float, np.ndarray],
    kappa: Union[float, np.ndarray],
) -> Union[float, np.ndarray]:
    r"""Goldston L-mode energy confinement time.

    $$\tau_{L} = 6.4\times10^{-8}\, I_p\, P^{-1/2} R^{1.75} a^{-0.37} \kappa^{1/2}$$

    in the paper's units ($I_p$ in A, $P$ in W, $R$ and $a$ in cm, $\tau$ in s).

    Parameters
    ----------
    I_p : float or np.ndarray
        Plasma current [A].
    P : float or np.ndarray
        Total heating (loss) power [W].
    R : float or np.ndarray
        Major radius, converted to cm internally [m].
    a : float or np.ndarray
        Minor radius, converted to cm internally [m].
    kappa : float or np.ndarray
        Elongation [-].

    Returns
    -------
    float or np.ndarray
        Energy confinement time [s].

    Raises
    ------
    ValueError
        A non-finite or non-positive input.

    Convention
    ----------
    Strict SI in; only $R$ and $a$ are converted (m to cm), since the paper
    already uses A and W.  This is Goldston's eq. (6) [1]_ without his isotope
    factor $(A_i/1.5)^{1/2}$, i.e. evaluated at $A_i = 1.5$.

    Validity
    --------
    Empirical fit.  Regressed on the auxiliary-heated L-mode data of 1984
    (PDX, ISX-B, ASDEX, Doublet III); no density dependence.

    Limitations
    -----------
    Pure L-mode; the ohmic phase needs the quadrature of
    :func:`ohmic_l_mode_confinement_time_from_tau_ohmic_tau_aux`.

    Semantics
    ---------
    consumes: plasma_current, major_radius, minor_radius, elongation
    produces: energy_confinement_time

    References
    ----------
    .. [1] R. J. Goldston, Plasma Phys. Control. Fusion 26 (1984) 87, Eq. (6).
    """
    I_A = _validate_positive("I_p", I_p)
    P_W = _validate_positive("P", P)
    R_cm = _validate_positive("R", R) * 100.0
    a_cm = _validate_positive("a", a) * 100.0
    kappa_arr = _validate_positive("kappa", kappa)
    tau = 6.4e-8 * I_A * P_W ** -0.5 * R_cm ** 1.75 * a_cm ** -0.37 * kappa_arr ** 0.5
    return float(tau) if np.ndim(tau) == 0 else tau

def grad_shafranov_source(R, J_phi):
    r"""Right-hand side of the poloidal-flux equation from Ampere's law, any region.

    $$\Delta^*\psi = -\mu_0 R J_\phi,\qquad
      \Delta^* = R\frac{\partial}{\partial R}\frac{1}{R}\frac{\partial}{\partial R}
      + \frac{\partial^2}{\partial Z^2}$$

    Parameters
    ----------
    R : float or np.ndarray
        Major radius, positive [m].
    J_phi : float or np.ndarray
        Toroidal current density: plasma, coil or zero [A/m^2].

    Returns
    -------
    float or np.ndarray
        $\Delta^*\psi$ for a per-radian $\psi$ [Wb/(rad m^2)].

    Raises
    ------
    ValueError
        ``R`` is not positive.

    Convention
    ----------
    Per-radian flux $\psi = RA_\phi$ with $\mathbf B_p = \nabla\psi\times\nabla\phi$
    and $(R, \phi, Z)$ right-handed, $\phi$ counter-clockwise from above
    (Freidberg; COCOS 3, $\sigma_{B_p} = -1$): a positive $J_\phi$ makes $\psi$
    a maximum on the magnetic axis. A flux stored in another COCOS is
    converted first with ``psi_per_radian_from_cocos`` -- for the full-weber
    COCOS 11 of an IMAS ODS, $\psi = -\psi_{11}/(2\pi)$. The operator itself
    on a grid is ``vaft.process.equilibrium.grad_shafranov_operator``.

    Physical interpretation
    -----------------------
    One elliptic equation for one flux function everywhere; only the source
    differs by region. In the plasma $J_\phi$ is constrained by force balance
    (``toroidal_current_density_from_p_prime_ff_prime``), in a coil it is the
    prescribed coil current density, and in vacuum it is zero, leaving the
    homogeneous equation $\Delta^*\psi = 0$ -- Laplace-type, but not the
    scalar Laplacian ($\Delta^*$ has $-R^{-1}\partial_R$ where $\nabla^2$ has
    $+R^{-1}\partial_R$).

    Assumptions
    -----------
    Axisymmetry, $\partial_\phi = 0$; magnetostatics (no displacement current).

    References
    ----------
    .. [1] J. P. Freidberg, *Ideal MHD*, Cambridge University Press (2014),
           Sec. 6.2.
    .. [2] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 3.3.
    """
    R = np.asarray(R, dtype=float)
    if np.any(R <= 0.0):
        raise ValueError("R must be positive")
    return -MU0 * R * np.asarray(J_phi, dtype=float)

def heating_power_from_p_ohm_p_aux(P_ohm: float, P_aux: float) -> float:
    r"""Total heating power $P_{heat} = P_{ohm} + P_{aux}$.

    $$P_{\mathrm{heat}} = P_{\mathrm{ohm}} + P_{\mathrm{aux}}$$

    Parameters
    ----------
    P_ohm : float
        Ohmic heating power [W].
    P_aux : float
        Absorbed auxiliary heating power [W].

    Returns
    -------
    float
        Total heating power [W].
    """
    return P_ohm + P_aux

def inductive_voltage_from_dW_magdt_I_p(dW_magdt: float, I_p: float) -> float:
    r"""Inductive voltage from the rate of change of magnetic energy.

    $$V_{\mathrm{ind}} = \frac{1}{I_p}\,\frac{dW_{\mathrm{mag}}}{dt}$$

    Parameters
    ----------
    dW_magdt : float
        Rate of change of the poloidal magnetic energy [W].
    I_p : float
        Plasma current [A].

    Returns
    -------
    float
        Inductive voltage [V].

    Convention
    ----------
    **Exact for the internal field, whether or not the inductance changes.**
    With $W_{\mathrm{mag}} = \tfrac12 L_i I_p^2$ the energy inside the
    plasma, this is Romero's $V_{\mathrm{ind}} = L_i\dot I_p +
    \tfrac12 I_p\dot L_i$ (eq. 24), and $V_B = R_p I_p + V_{\mathrm{ind}}$
    closes exactly.  The circuit form $\mathrm{d}(L_i I_p)/\mathrm{d}t$ is
    *not* the larger "full" voltage it is sometimes taken for: $L_i I_p$ is not
    a flux linked by one loop, and that form overstates the profile term by
    $\tfrac12 I_p\dot L_i$.  The external inductance is different -- there
    $\mathrm{d}(L_e I_p)/\mathrm{d}t$ is right, because $L_e I_p$ is the
    boundary flux.

    References
    ----------
    .. [1] J. A. Romero and JET-EFDA contributors, Nucl. Fusion 50 (2010)
           115002, eqs. (22)-(24).

    See Also
    --------
    vaft.formula.transformer.plasma_current_rate_from_L_i_V_B_V_C_V_R
    """
    return dW_magdt / I_p

def internal_inductance_from_W_int_Ip(W_int: float, Ip: float) -> float:
    r"""Dimensional internal inductance from the poloidal-field energy inside the plasma.

    $$L_i = \frac{2W_{p,\mathrm{int}}}{I_p^2},\qquad
      W_{p,\mathrm{int}} = \int_{V_p}\frac{B_p^2}{2\mu_0}\,dV$$

    Parameters
    ----------
    W_int : float
        Poloidal magnetic energy inside the last closed flux surface [J].
    Ip : float
        Plasma current; only its magnitude enters [A].

    Returns
    -------
    float
        Internal inductance [H].

    Raises
    ------
    ValueError
        Zero or non-finite plasma current.

    Convention
    ----------
    **Dimensional, in henry**, and only the field *inside* the plasma counts:
    the energy outside contains plasma-coil cross terms, which is why the
    external inductance is defined from the boundary flux and not from
    $2W_{\mathrm{outside}}/I_p^2$.  Its dimensionless counterpart is the IMAS
    ``li_3`` (:func:`li_3_from_internal_inductance_R0`); the same energy
    written as $\int B_p^2\,dV = 2\mu_0 W_{p,\mathrm{int}}$ gives
    :func:`li_3_from_Bp2_volume_integral` directly.

    References
    ----------
    .. [1] J. A. Romero and JET-EFDA contributors, Nucl. Fusion 50 (2010)
           115002, Sec. II, eq. (22).
    .. [2] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Ch. 3, Equilibrium (internal inductance).
    """
    current = float(Ip)
    if not np.isfinite(current) or current == 0.0:
        raise ValueError(f"Ip must be finite and non-zero; got {Ip!r}")
    return 2.0 * float(W_int) / current**2

def internal_inductance_from_li_3_R0(li_3: float, R0: float) -> float:
    r"""Dimensional internal inductance from the IMAS ``li_3``.

    $$L_i = \frac{\mu_0 R_0}{2}\,l_{i3}$$

    Parameters
    ----------
    li_3 : float
        Internal inductance in the IMAS ``li_3`` definition [-].
    R0 : float
        Major radius ``li_3`` was normalised by [m].

    Returns
    -------
    float
        Internal inductance [H].

    Convention
    ----------
    **Exact only for** ``li_3``, and only with the $R_0$ the equilibrium
    normalised by: $l_{i3} = 2\int B_p^2\,dV/(\mu_0^2 I_p^2 R_0)$ is
    $2L_i/(\mu_0 R_0)$ by definition.  ``li_1`` normalises by the edge
    poloidal field instead, $l_{i1}/l_{i3} = L_{pol}^2 R_0/(2V)$, so it is not
    a valid input here -- about 8.5 % too large at $\kappa = 1.6$.

    References
    ----------
    .. [1] IMAS Data Dictionary, ``equilibrium.time_slice[:].global_quantities.li_3``.
    .. [2] J. A. Romero and JET-EFDA contributors, Nucl. Fusion 50 (2010)
           115002, eq. (71), which normalises by the magnetic-axis radius.
    """
    return 0.5 * MU0 * float(R0) * float(li_3)

def inverse_aspect_ratio_from_a_R(a: float, R: float) -> float:
    r"""Inverse aspect ratio $\varepsilon = a/R$.

    $$\varepsilon = \frac{a}{R}$$

    Parameters
    ----------
    a : float
        Minor radius [m].
    R : float
        Major radius [m].

    Returns
    -------
    float
        Inverse aspect ratio [-].

    See Also
    --------
    calc_inverse_aspect_ratio : the same ratio with positivity validation.
    """
    return a / R

def ion_pressure(n_i: Union[float, np.ndarray],
                 T_i: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
    r"""Ion thermal pressure of one species from density and temperature, $p_i = n_i T_i$.

    $$p_i = n_i\,k_B T_i = n_i\,T_i[\mathrm{eV}]\;e$$

    Parameters
    ----------
    n_i : float or np.ndarray
        Density of one ion species, finite and non-negative [m^-3].
    T_i : float or np.ndarray
        Temperature of that species, finite and non-negative, broadcastable
        against ``n_i`` [eV].

    Returns
    -------
    float or np.ndarray
        Pressure of that species [Pa].

    Raises
    ------
    ValueError
        A non-finite or negative density or temperature.

    Convention
    ----------
    Same law and units as :func:`electron_pressure`.  The density is the
    *ion* density of that species, not $n_e$: quasi-neutrality gives
    $\sum_i Z_i n_i = n_e$.  Sum the species yourself for the total ion
    pressure.

    Assumptions
    -----------
    Isotropic Maxwellian ions of the stated temperature; beam or
    fusion-product ions are a separate, non-thermal pressure.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Ch. 2 (the ideal-gas pressure of a plasma species).

    See Also
    --------
    electron_pressure
    """
    return _thermal_pressure(n_i, T_i, n_name="n_i", T_name="T_i")

def kadomtsev_constraint_from_engineering_exponents(
    alpha_I: Union[float, np.ndarray],
    alpha_B: Union[float, np.ndarray],
    alpha_P: Union[float, np.ndarray],
    alpha_n: Union[float, np.ndarray],
    alpha_R: Union[float, np.ndarray],
) -> Union[float, np.ndarray]:
    r"""Residual of the Kadomtsev high-beta constraint on engineering exponents.

    $$\alpha_K = 4\alpha_R - 8\alpha_n - \alpha_I - 3\alpha_P - 5\alpha_B - 5$$

    for $\tau_E \propto I^{\alpha_I}B^{\alpha_B}P^{\alpha_P}n^{\alpha_n}R^{\alpha_R}$;
    $\alpha_K = 0$ when the scaling is expressible in the three dimensionless
    parameters $\rho_*$, $\beta$, $\nu_*$ alone.

    Parameters
    ----------
    alpha_I : float or np.ndarray
        Exponent of the plasma current [-].
    alpha_B : float or np.ndarray
        Exponent of the toroidal field [-].
    alpha_P : float or np.ndarray
        Exponent of the heating power [-].
    alpha_n : float or np.ndarray
        Exponent of the density [-].
    alpha_R : float or np.ndarray
        Exponent of the major radius [-].

    Returns
    -------
    float or np.ndarray
        Constraint residual, 0 for exact satisfaction [-].

    Raises
    ------
    ValueError
        For non-finite input.

    Physical interpretation
    -----------------------
    Dimensional analysis of the Vlasov-Maxwell system: only three of the
    engineering variables are independent once $\rho_*$, $\beta$ and $\nu_*$ are
    fixed at constant geometry.  ITER89P gives $\alpha_K = -0.15$ and IPB98(y,2)
    gives $-0.01$.

    References
    ----------
    .. [1] B. B. Kadomtsev, Sov. J. Plasma Phys. 1 (1975) 295.
    .. [2] ITER Physics Expert Groups, Nucl. Fusion 39 (1999) 2175, Ch. 2,
           Sec. 6.2 (Kadomtsev constraint).
    .. [3] G. Verdoolaege et al., Nucl. Fusion 61 (2021) 076006.
    """
    a_i = np.asarray(alpha_I, dtype=float)
    a_b = np.asarray(alpha_B, dtype=float)
    a_p = np.asarray(alpha_P, dtype=float)
    a_n = np.asarray(alpha_n, dtype=float)
    a_r = np.asarray(alpha_R, dtype=float)
    for name, arr in (
        ("alpha_I", a_i),
        ("alpha_B", a_b),
        ("alpha_P", a_p),
        ("alpha_n", a_n),
        ("alpha_R", a_r),
    ):
        if np.any(~np.isfinite(arr)):
            raise ValueError(f"{name} must be finite. Got {arr!r}")
    alpha_k = 4.0 * a_r - 8.0 * a_n - a_i - 3.0 * a_p - 5.0 * a_b - 5.0
    return _maybe_scalar(alpha_k)