def bounce_action(s, p_parallel, *, closure_tolerance=1e-6):
    r"""Unnormalized parallel action on one explicitly closed bounce cycle.

    $$J_\parallel=\oint p_\parallel\,\mathrm ds$$

    Parameters
    ----------
    s : array_like
        Signed field-line arc coordinate in time order, returning to its
        initial value after exactly one full bounce [m].
    p_parallel : array_like
        Parallel momentum at the same samples, returning to its initial value
        and changing sign on both legs [kg m/s].
    closure_tolerance : float, optional
        Allowed endpoint mismatch relative to the sampled excursion,
        a numerical closure tolerance in [0, 0.01] [-].

    Returns
    -------
    float
        Parallel action without a factor of 1/(2 pi) [J s].

    Convention
    ----------
    s increases along B, not accumulated unsigned distance. Endpoint closure,
    both momentum signs and positive oriented action are required. The caller
    must identify one physical cycle; these checks cannot establish periodicity.

    Limitations
    -----------
    Trapezoidal quadrature needs a resolved complete periodic bounce. Partial
    histories, passing transits and separatrix orbits have no action here.

    References
    ----------
    .. [1] J. R. Cary and A. J. Brizard, Rev. Mod. Phys. 81 (2009) 693.
    """
    s, p = np.asarray(s, float), np.asarray(p_parallel, float)
    if s.ndim != 1 or p.shape != s.shape or len(s) < 5 or not np.all(np.isfinite(s + p)):
        raise ValueError("a finite resolved one-dimensional bounce cycle is required")
    if np.ptp(s) == 0 or np.ptp(p) == 0:
        raise ValueError("degenerate bounce cycle")
    if not np.isfinite(closure_tolerance) or not 0 <= closure_tolerance <= .01:
        raise ValueError("closure_tolerance must be finite and in [0,0.01]")
    if abs(s[-1] - s[0]) > closure_tolerance * np.ptp(s) or abs(p[-1] - p[0]) > closure_tolerance * np.ptp(p):
        raise ValueError("bounce cycle is not closed")
    if not (np.min(p) < 0 < np.max(p)):
        raise ValueError("bounce cycle must contain both parallel directions")
    result = float(np.sum(0.5 * (p[1:] + p[:-1]) * np.diff(s)))
    if result <= 0:
        raise ValueError("s must be oriented along the signed parallel momentum")
    return result

def bounce_harmonic_detuning(omega_b, omega_exb, omega_magnetic=0.0, *, l=1, n=1):
    r"""Frequency mismatch of a trapped orbit with a toroidal-$n$ perturbation at bounce harmonic $\ell$.

    $$\Delta\omega_\mathrm{BH} = \ell\,\omega_b - n\,\omega_E - n\,\omega_B$$

    Parameters
    ----------
    omega_b : float or np.ndarray
        Bounce frequency of the trapped orbit [rad/s].
    omega_exb : float or np.ndarray
        Toroidal $E\times B$ precession frequency [rad/s].
    omega_magnetic : float or np.ndarray
        Toroidal magnetic (grad-B and curvature) precession frequency [rad/s].
    l : int
        Bounce harmonic, any integer including 0 [-].
    n : int
        Toroidal mode number of the perturbation [-].

    Returns
    -------
    float or np.ndarray
        $\Delta\omega_\mathrm{BH}$; zero at the resonance [rad/s].

    Raises
    ------
    ValueError
        ``l`` or ``n`` is not an integer.

    Convention
    ----------
    The temporal resonance of Park, Boozer and Menard: a static perturbation
    $\propto e^{in\phi}$ seen by an orbit precessing toroidally at
    $\omega_E + \omega_B$ and bouncing at $\omega_b$. $\ell = 0$ is the
    precession resonance. It is not the *spatial* resonance $m = nq$ of a
    field line; the two are separate conditions.

    Physical interpretation
    -----------------------
    Where it vanishes the orbit sees a stationary phase of the perturbation
    each bounce, so the change of $P_\phi$ adds up secularly instead of
    averaging out -- the orbit-level origin of non-resonant transport and NTV,
    which themselves need the kinetic response of the whole distribution.

    References
    ----------
    .. [1] J.-K. Park, A. H. Boozer and J. E. Menard, Phys. Rev. Lett. 102
           (2009) 065002.
    """
    for name, value in (("l", l), ("n", n)):
        if isinstance(value, bool) or int(value) != value:
            raise ValueError(f"{name} must be an integer, not {value!r}")
    result = (int(l) * np.asarray(omega_b, dtype=float) - int(n) * np.asarray(omega_exb, dtype=float)
              - int(n) * np.asarray(omega_magnetic, dtype=float))
    return float(result) if np.ndim(result) == 0 else result

def canonical_toroidal_momentum(q, m, R, v_phi, A_phi):
    r"""Canonical toroidal angular momentum of a particle in an axisymmetric field.

    $$P_\phi = m R v_\phi + q R A_\phi$$

    Parameters
    ----------
    q : float
        Signed particle charge [C].
    m : float
        Particle mass [kg].
    R : float or np.ndarray
        Major radius of the particle [m].
    v_phi : float or np.ndarray
        Physical toroidal velocity component $\hat{\boldsymbol\phi}\cdot\mathbf v$ [m/s].
    A_phi : float or np.ndarray
        Physical toroidal component of the vector potential [T m].

    Returns
    -------
    float or np.ndarray
        $P_\phi = \partial\mathcal L/\partial\dot\phi$ [kg m^2/s].

    Raises
    ------
    ValueError
        ``m`` is not positive or ``R`` is negative.

    Convention
    ----------
    Cylindrical $(R, \phi, Z)$ right-handed, $\phi$ counter-clockwise seen
    from above (the IMAS $\phi$). ``v_phi`` and ``A_phi`` are *physical*
    components; $RA_\phi$ is the covariant one, which equals the poloidal flux
    per radian $\psi$ (see ``guiding_center_toroidal_momentum``). The full
    particle form: it includes the gyration in $v_\phi$.

    Physical interpretation
    -----------------------
    The momentum conjugate to $\phi$. Where $\partial\mathcal L/\partial\phi = 0$
    (axisymmetry) it is exactly conserved, which ties the particle's toroidal
    velocity to its flux-surface position.

    Assumptions
    -----------
    The Lagrangian $\mathcal L = \tfrac12 mv^2 + q\mathbf A\cdot\mathbf v - q\Phi$;
    conservation needs axisymmetry, $\partial\mathcal L/\partial\phi = 0$.

    References
    ----------
    .. [1] R. B. White, *The Theory of Toroidally Confined Plasmas*, 3rd ed.,
           Imperial College Press (2014), Ch. 3.
    """
    m = float(m)
    R = np.asarray(R, dtype=float)
    if m <= 0.0 or np.any(R < 0.0):
        raise ValueError("m must be positive and R non-negative")
    result = m * R * np.asarray(v_phi, dtype=float) + float(q) * R * np.asarray(A_phi, dtype=float)
    return float(result) if np.ndim(result) == 0 else result

def curvature_drift_velocity(q, m, v_par, B, R_c):
    r"""Curvature drift velocity of the guiding centre.

    $$\mathbf{v}_R = \frac{m v_\parallel^{2}}{qB^{2}}\,
    \frac{\mathbf{R}_c\times\mathbf{B}}{R_c^{2}}$$

    Parameters
    ----------
    q : float
        Particle charge, signed [C].
    m : float
        Particle mass [kg].
    v_par : float or np.ndarray
        Speed along the magnetic field [m/s].
    B : array_like
        Magnetic field, last axis the three Cartesian components [T].
    R_c : array_like
        Radius-of-curvature vector, from the centre of curvature to the
        field line [m].

    Returns
    -------
    np.ndarray
        Drift velocity [m/s].

    Raises
    ------
    ValueError
        An input is not a 3-vector, ``q`` is zero, or ``B`` or ``R_c``
        vanishes.

    Convention
    ----------
    $\mathbf{R}_c$ points *outward*, from the centre of curvature to the
    line, as in Chen. Charge-signed like the grad-B drift. In a vacuum field
    $\nabla\times\mathbf{B} = 0$ the two drifts add in the same direction.

    Physical interpretation
    -----------------------
    The centrifugal force of streaming along a curved line, $m v_\parallel^2
    / R_c$, acts like any other perpendicular force $\mathbf{F}$ and drives
    $\mathbf{F}\times\mathbf{B}/(qB^2)$.

    Assumptions
    -----------
    $\rho \ll R_c$.

    References
    ----------
    .. [1] F. F. Chen, *Introduction to Plasma Physics and Controlled
           Fusion*, 3rd ed., Springer (2016), Sec. 2.3.2.
    """
    B = _vector(B, "B")
    R_c = _vector(R_c, "R_c")
    q = float(q)
    if q == 0.0:
        raise ValueError("q must be non-zero")
    B2 = np.sum(B * B, axis=-1, keepdims=True)
    R2 = np.sum(R_c * R_c, axis=-1, keepdims=True)
    if np.any(B2 == 0.0) or np.any(R2 == 0.0):
        raise ValueError("B and R_c must be non-zero")
    v_par = np.asarray(v_par, dtype=float)[..., None] if np.ndim(v_par) else float(v_par)
    return float(m) * v_par ** 2 / (q * B2) * np.cross(R_c, B) / R2

def effective_parallel_potential(mu, B, q=0.0, potential=0.0):
    r"""Effective parallel potential of a magnetized particle in static fields.

    $$U_\mathrm{eff}=\mu B+q\Phi$$

    Parameters
    ----------
    mu : float or np.ndarray
        Nonnegative magnetic moment [J/T].
    B : float or np.ndarray
        Positive magnetic field magnitude [T].
    q : float, optional
        Signed particle charge; zero omits the electrostatic term [C].
    potential : float or np.ndarray, optional
        Electrostatic potential in a stated gauge [V].

    Returns
    -------
    float or np.ndarray
        Effective potential energy [J].

    Assumptions
    -----------
    Static electrostatic field and adiabatic magnetic moment.

    References
    ----------
    .. [1] F. F. Chen, Introduction to Plasma Physics and Controlled Fusion,
       3rd ed. (2016), Sec. 2.3.
    """
    mu, B, q, potential = np.broadcast_arrays(*[np.asarray(v, float) for v in (mu, B, q, potential)])
    if np.any(mu < 0) or np.any(B <= 0) or not np.all(np.isfinite(mu + B + q + potential)):
        raise ValueError("finite inputs, mu >= 0 and B > 0 required")
    return mu * B + q * potential

def exb_drift_velocity(E, B):
    r"""E-cross-B drift velocity of the guiding centre.

    $$\mathbf{v}_E = \frac{\mathbf{E}\times\mathbf{B}}{B^{2}}$$

    Parameters
    ----------
    E : array_like
        Electric field, last axis the three Cartesian components [V/m].
    B : array_like
        Magnetic field, last axis the three Cartesian components [T].

    Returns
    -------
    np.ndarray
        Drift velocity, same shape as the broadcast inputs [m/s].

    Raises
    ------
    ValueError
        An input is not a 3-vector, or ``B`` vanishes.

    Convention
    ----------
    Right-handed Cartesian components. The drift is independent of charge
    and mass, so ions and electrons move together and no current flows.

    Physical interpretation
    -----------------------
    Over each gyration the particle gains energy on one side of its orbit
    and loses it on the other, so its radius changes around the orbit and
    the guiding centre steps sideways -- perpendicular to both fields.

    Assumptions
    -----------
    $E_\perp \ll cB$ (non-relativistic drift); uniform fields over a Larmor
    radius and slow in time compared with the gyroperiod.

    References
    ----------
    .. [1] F. F. Chen, *Introduction to Plasma Physics and Controlled
           Fusion*, 3rd ed., Springer (2016), Sec. 2.2.2.
    """
    E = _vector(E, "E")
    B = _vector(B, "B")
    B2 = np.sum(B * B, axis=-1, keepdims=True)
    if np.any(B2 == 0.0):
        raise ValueError("B must be non-zero")
    return np.cross(E, B) / B2

def grad_b_drift_velocity(q, m, v_perp, B, grad_B):
    r"""Grad-B drift velocity of the guiding centre.

    $$\mathbf{v}_{\nabla B} = \frac{m v_\perp^{2}}{2qB}\,
    \frac{\mathbf{B}\times\nabla B}{B^{2}}$$

    Parameters
    ----------
    q : float
        Particle charge, signed [C].
    m : float
        Particle mass [kg].
    v_perp : float or np.ndarray
        Speed perpendicular to the magnetic field [m/s].
    B : array_like
        Magnetic field, last axis the three Cartesian components [T].
    grad_B : array_like
        Gradient of the field magnitude $|\mathbf{B}|$ [T/m].

    Returns
    -------
    np.ndarray
        Drift velocity [m/s].

    Raises
    ------
    ValueError
        An input is not a 3-vector, ``q`` is zero, or ``B`` vanishes.

    Convention
    ----------
    Charge-signed: ions and electrons drift in opposite directions, which is
    what makes this drift carry current.

    Physical interpretation
    -----------------------
    The orbit is tighter where the field is stronger, so the guiding centre
    walks perpendicular to both $\mathbf{B}$ and its gradient.

    Assumptions
    -----------
    $\rho \ll B/|\nabla B|$; the drift is the first-order average over a
    gyration.

    References
    ----------
    .. [1] F. F. Chen, *Introduction to Plasma Physics and Controlled
           Fusion*, 3rd ed., Springer (2016), Sec. 2.3.1.
    """
    B = _vector(B, "B")
    grad_B = _vector(grad_B, "grad_B")
    q = float(q)
    if q == 0.0:
        raise ValueError("q must be non-zero")
    Bmag = np.linalg.norm(B, axis=-1, keepdims=True)
    if np.any(Bmag == 0.0):
        raise ValueError("B must be non-zero")
    v_perp = np.asarray(v_perp, dtype=float)[..., None] if np.ndim(v_perp) else float(v_perp)
    return float(m) * v_perp ** 2 / (2.0 * q * Bmag) * np.cross(B, grad_B) / Bmag ** 2

def guiding_center_rhs(q, m, p_parallel, mu, B, grad_B, curl_b, E):
    r"""Hamiltonian guiding-center velocity and parallel force in static fields.

    $$\mathbf B^*=\mathbf B+\frac{p_\parallel}{q}\nabla\times\mathbf b,\quad
    \dot{\mathbf X}=\frac{(p_\parallel/m)\mathbf B^*
    +\mathbf b\times(\mu\nabla B-q\mathbf E)/q}{\mathbf b\cdot\mathbf B^*},\quad
    \dot p_\parallel=\frac{\mathbf B^*\cdot(q\mathbf E-\mu\nabla B)}
    {\mathbf b\cdot\mathbf B^*}$$

    Parameters
    ----------
    q : float
        Signed nonzero charge [C].
    m : float
        Positive particle mass [kg].
    p_parallel : float
        Signed parallel kinetic momentum [kg m/s].
    mu : float
        Nonnegative magnetic moment, held fixed by the reduced model [J/T].
    B : array_like
        Cartesian magnetic field, shape (3,) [T].
    grad_B : array_like
        Gradient of field magnitude, shape (3,) [T/m].
    curl_b : array_like
        Curl of unit field direction, shape (3,) [1/m].
    E : array_like
        Static electric field, shape (3,) [V/m].

    Returns
    -------
    velocity : np.ndarray
        Guiding-center Cartesian velocity [m/s].
    force : float
        Parallel kinetic momentum derivative [N].

    Convention
    ----------
    SI units, right-handed Cartesian axes; b = B/abs(B). The low-frequency,
    weak-electric-drift Hamiltonian is p_parallel^2/(2m) + mu B + q Phi.
    No separate drift velocities are added to this expression.

    Validity
    --------
    Nonrelativistic magnetized particles in smooth static fields with small
    gyroradius compared with field scales, and weak E-cross-B flow compared
    with characteristic particle speed. This is the lowest-order Hamiltonian
    GC model, not a high-flow or higher-order finite-Larmor-radius theory.

    Limitations
    -----------
    The B-star parallel singularity must terminate integration. Holding mu
    fixed is a model assumption, not a numerical proof of adiabaticity.

    References
    ----------
    .. [1] J. R. Cary and A. J. Brizard, Rev. Mod. Phys. 81 (2009) 693,
       doi:10.1103/RevModPhys.81.693.
    """
    if not all(np.isfinite([q, m, p_parallel, mu])) or q == 0 or m <= 0 or mu < 0:
        raise ValueError("finite q != 0, m > 0 and mu >= 0 required")
    vectors = [np.asarray(a, dtype=float) for a in (B, grad_B, curl_b, E)]
    if any(a.shape != (3,) or not np.all(np.isfinite(a)) for a in vectors):
        raise ValueError("B, grad_B, curl_b and E must be finite 3-vectors")
    B, grad_B, curl_b, E = vectors
    strength = np.linalg.norm(B)
    if strength == 0:
        raise ValueError("magnetic field must be nonzero")
    b = B / strength
    star = B + p_parallel / q * curl_b
    denominator = float(b @ star)
    if abs(denominator) <= 1e-10 * strength:
        raise ValueError("guiding-center B-star parallel singularity")
    force = q * E - mu * grad_B
    return ((p_parallel / m * star - np.cross(b, force) / q) / denominator,
            float(star @ force / denominator))

def guiding_center_toroidal_momentum(q, m, v_par, R, b_phi, psi_per_radian):
    r"""Canonical toroidal momentum of a guiding centre, flux plus mechanical parts.

    $$P_{\phi,\mathrm{gc}} = q\psi + m v_\parallel R\, b_\phi$$

    Parameters
    ----------
    q : float
        Signed particle charge [C].
    m : float
        Particle mass [kg].
    v_par : float or np.ndarray
        Signed parallel velocity, positive along $\mathbf B$ [m/s].
    R : float or np.ndarray
        Major radius of the guiding centre [m].
    b_phi : float or np.ndarray
        Toroidal component of the unit vector $\mathbf b = \mathbf B/B$ [-].
    psi_per_radian : float or np.ndarray
        Poloidal flux per radian at the guiding centre, $\psi = RA_\phi$ [Wb/rad].

    Returns
    -------
    float or np.ndarray
        $P_{\phi,\mathrm{gc}}$ [kg m^2/s].

    Raises
    ------
    ValueError
        ``m`` is not positive or ``R`` is negative.

    Convention
    ----------
    $\psi$ is **per radian** and is the covariant $RA_\phi$ with the IMAS
    $\phi$: its sign follows the plasma current's direction in $\phi$, largest
    on the axis for a current along $+\phi$. A stored COCOS flux is converted
    by ``psi_per_radian_from_cocos`` -- $-\psi/(2\pi)$ for COCOS 11 (IMAS DD3),
    $+\psi/(2\pi)$ for COCOS 17; used as stored it gives the flux term the
    wrong sign or a $2\pi$ error. $b_\phi v_\parallel$ is
    the toroidal velocity of the guiding centre; $mRb_\phi v_\parallel = mIv_\parallel/B$
    with $I = RB_\phi$.

    Physical interpretation
    -----------------------
    Conserved in axisymmetry, so a change of $v_\parallel$ along the orbit is paid
    for by a change of $\psi$: $\Delta\psi = -(m/q)\Delta(v_\parallel Rb_\phi)$.
    That is the finite orbit width -- the global view of what the grad-B and
    curvature drifts do locally.

    Assumptions
    -----------
    Lowest-order guiding-centre theory: the gyration averaged out, drift
    velocity small against $v_\parallel$ in the mechanical term.

    References
    ----------
    .. [1] R. B. White, *The Theory of Toroidally Confined Plasmas*, 3rd ed.,
           Imperial College Press (2014), Ch. 3.
    .. [2] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 3.12.
    """
    m = float(m)
    R = np.asarray(R, dtype=float)
    if m <= 0.0 or np.any(R < 0.0):
        raise ValueError("m must be positive and R non-negative")
    result = (float(q) * np.asarray(psi_per_radian, dtype=float)
              + m * np.asarray(v_par, dtype=float) * R * np.asarray(b_phi, dtype=float))
    return float(result) if np.ndim(result) == 0 else result

def gyration_offset(q, m, B, v):
    r"""Position of a gyrating particle relative to its guiding centre.

    $$\boldsymbol{\rho} = \frac{m}{qB^{2}}\,\mathbf{B}\times\mathbf{v}_\perp$$

    Parameters
    ----------
    q : float
        Particle charge, signed [C].
    m : float
        Particle mass [kg].
    B : array_like
        Magnetic field, last axis the three Cartesian components [T].
    v : array_like
        Velocity in the guiding-centre frame; only its part perpendicular to
        ``B`` enters [m/s].

    Returns
    -------
    np.ndarray
        Vector from the guiding centre to the particle [m].

    Raises
    ------
    ValueError
        An input is not a 3-vector, ``q`` is zero, ``m`` is not positive, or
        ``B`` vanishes.

    Convention
    ----------
    Guiding centre = position minus this vector. Its magnitude is
    ``larmor_radius`` and its sense follows the charge: for $\mathbf{B} =
    \hat z$ and $\mathbf{v} = \hat x$, an ion sits a Larmor radius along
    $+\hat y$ from its guiding centre, an electron along $-\hat y$.

    Physical interpretation
    -----------------------
    The Lorentz force $q\mathbf{v}\times\mathbf{B}$ points from the particle
    to the centre of its orbit; this is minus that direction times the
    Larmor radius.

    Assumptions
    -----------
    Uniform field over the orbit; ``v`` measured in the frame moving with any
    drift.

    References
    ----------
    .. [1] F. F. Chen, *Introduction to Plasma Physics and Controlled
           Fusion*, 3rd ed., Springer (2016), Sec. 2.2.
    """
    B = _vector(B, "B")
    v = _vector(v, "v")
    q = float(q)
    if q == 0.0 or not float(m) > 0.0:
        raise ValueError("q must be non-zero and m positive")
    B2 = np.sum(B * B, axis=-1, keepdims=True)
    if np.any(B2 == 0.0):
        raise ValueError("B must be non-zero")
    return float(m) / (q * B2) * np.cross(B, v)

def gyrofrequency(q, m, B):
    r"""Signed gyrofrequency of a charged particle.

    $$\omega_c = \frac{qB}{m}$$

    Parameters
    ----------
    q : float or np.ndarray
        Particle charge, signed [C].
    m : float or np.ndarray
        Particle mass [kg].
    B : float or np.ndarray
        Magnetic-field magnitude [T].

    Returns
    -------
    float or np.ndarray
        Gyrofrequency, positive for positive charge [rad/s].

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

    Convention
    ----------
    Signed: $\omega_c > 0$ for an ion, $< 0$ for an electron. An ion gyrates
    in the left-handed sense about $\mathbf{B}$ (clockwise seen with
    $\mathbf{B}$ pointing at the viewer), an electron in the right-handed
    sense; the sign carries that. Use ``abs`` for the frequency alone.

    Physical interpretation
    -----------------------
    The rate at which the Lorentz force turns the perpendicular velocity:
    the fastest time scale of magnetised motion, and the reason drifts are
    slow by comparison.

    References
    ----------
    .. [1] F. F. Chen, *Introduction to Plasma Physics and Controlled
           Fusion*, 3rd ed., Springer (2016), Sec. 2.2.
    .. [2] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 2.2.
    """
    m = np.asarray(m, dtype=float)
    if np.any(m <= 0.0):
        raise ValueError("m must be positive")
    return np.asarray(q, dtype=float) * np.asarray(B, dtype=float) / m

def larmor_radius(q, m, v_perp, B):
    r"""Larmor (gyro) radius of a charged particle.

    $$\rho = \frac{m\,v_\perp}{|q|\,B}$$

    Parameters
    ----------
    q : float or np.ndarray
        Particle charge; only its magnitude enters [C].
    m : float or np.ndarray
        Particle mass [kg].
    v_perp : float or np.ndarray
        Speed perpendicular to the magnetic field [m/s].
    B : float or np.ndarray
        Magnetic-field magnitude [T].

    Returns
    -------
    float or np.ndarray
        Radius of the gyro-orbit [m].

    Raises
    ------
    ValueError
        ``q`` or ``B`` is zero.

    Convention
    ----------
    Always positive: the gyration sense is in ``gyrofrequency``. $v_\perp$ is
    the speed in the guiding-centre frame, i.e. after removing any drift.

    Physical interpretation
    -----------------------
    The scale below which a particle does not follow field lines; drift
    approximations need it small against the field's gradient length.

    References
    ----------
    .. [1] F. F. Chen, *Introduction to Plasma Physics and Controlled
           Fusion*, 3rd ed., Springer (2016), Sec. 2.2.
    """
    q = np.abs(np.asarray(q, dtype=float))
    B = np.abs(np.asarray(B, dtype=float))
    if np.any(q == 0.0) or np.any(B == 0.0):
        raise ValueError("q and B must be non-zero")
    return np.asarray(m, dtype=float) * np.abs(np.asarray(v_perp, dtype=float)) / (q * B)

def magnetic_moment(m, v_perp, B):
    r"""Magnetic moment of a gyrating particle, the first adiabatic invariant.

    $$\mu = \frac{m v_\perp^2}{2B}$$

    Parameters
    ----------
    m : float
        Particle mass [kg].
    v_perp : float or np.ndarray
        Speed perpendicular to $\mathbf B$ [m/s].
    B : float or np.ndarray
        Field strength at the guiding centre [T].

    Returns
    -------
    float or np.ndarray
        $\mu$, non-negative [J/T].

    Raises
    ------
    ValueError
        ``m`` or ``B`` is not positive.

    Convention
    ----------
    Independent of the charge sign; ``v_perp`` enters squared, so its sign is
    irrelevant. With $E = \tfrac12 mv_\parallel^2 + \mu B$ in a static field it
    gives the mirror relation of ``parallel_speed_from_mu``.

    Physical interpretation
    -----------------------
    The gyro-orbit's magnetic dipole moment, equivalently its gyro-action
    $J_\perp = (2m/|q|)\mu$ up to a constant: the invariant of the fastest
    motion, which is why $v_\perp^2$ tracks $B$ along a field line.

    Assumptions
    -----------
    Conserved only under the adiabatic ordering: the field changes little over
    a gyro-period and a Larmor radius. It is defined for any inputs; its
    conservation is not implied.

    References
    ----------
    .. [1] F. F. Chen, *Introduction to Plasma Physics and Controlled Fusion*,
           3rd ed., Springer (2016), Sec. 2.3.3.
    """
    m = float(m)
    B = np.asarray(B, dtype=float)
    if m <= 0.0 or np.any(B <= 0.0):
        raise ValueError("m and B must be positive")
    result = 0.5 * m * np.asarray(v_perp, dtype=float) ** 2 / B
    return float(result) if np.ndim(result) == 0 else result

def mirror_force(mu, dB_ds):
    r"""Parallel magnetic mirror force along field-line arc length.

    $$F_\parallel=-\mu\,\frac{\mathrm dB}{\mathrm ds}$$

    Parameters
    ----------
    mu : float or np.ndarray
        Nonnegative magnetic moment [J/T].
    dB_ds : float or np.ndarray
        Signed derivative along the magnetic field [T/m].

    Returns
    -------
    float or np.ndarray
        Parallel force [N].

    Assumptions
    -----------
    Adiabatic magnetized motion, with s increasing along B.

    References
    ----------
    .. [1] F. F. Chen, Introduction to Plasma Physics and Controlled Fusion,
       3rd ed. (2016), Sec. 2.3.
    """
    mu, gradient = np.broadcast_arrays(np.asarray(mu, float), np.asarray(dB_ds, float))
    if np.any(mu < 0) or not np.all(np.isfinite(mu + gradient)):
        raise ValueError("finite gradient and nonnegative finite mu required")
    return -mu * gradient

def mirror_loss_cone(B_min, B_max):
    r"""Loss-cone pitch angle at the minimum of a magnetic mirror.

    $$\alpha_\mathrm{lc}=\arcsin\sqrt{B_\min/B_\max}$$

    Parameters
    ----------
    B_min : float or np.ndarray
        Positive minimum field [T].
    B_max : float or np.ndarray
        Maximum field, at least the minimum [T].

    Returns
    -------
    float or np.ndarray
        Loss-cone half angle; equality is the marginal separatrix [rad].

    Convention
    ----------
    Pitch is measured from either parallel direction at B_min. Particles
    inside either cone pass; outside they reflect in the adiabatic model.

    Assumptions
    -----------
    Static mirror without electric potential and conserved magnetic moment.

    References
    ----------
    .. [1] F. F. Chen, Introduction to Plasma Physics and Controlled Fusion,
       3rd ed. (2016), Sec. 2.3.
    """
    lo, hi = np.broadcast_arrays(np.asarray(B_min, float), np.asarray(B_max, float))
    if np.any(lo <= 0) or np.any(hi < lo) or not np.all(np.isfinite(lo + hi)):
        raise ValueError("finite 0 < B_min <= B_max required")
    return np.arcsin(np.sqrt(lo / hi))

def parallel_speed_from_mu(v, pitch_ref, B_ref, B):
    r"""Parallel speed where the field is $B$, from energy and magnetic-moment conservation.

    $$\frac{v_\parallel}{v} = \pm\sqrt{1 - \left(1 - \xi_\mathrm{ref}^2\right)\frac{B}{B_\mathrm{ref}}}$$

    Parameters
    ----------
    v : float
        Speed, conserved [m/s].
    pitch_ref : float
        Pitch $\xi = v_\parallel/v$ where the field is ``B_ref``, in $[-1, 1]$ [-].
    B_ref : float
        Field strength at the reference point [T].
    B : float or np.ndarray
        Field strength along the orbit [T].

    Returns
    -------
    float or np.ndarray
        $|v_\parallel|$, NaN where the particle cannot reach (it is reflected
        before $B$) [m/s].

    Raises
    ------
    ValueError
        ``pitch_ref`` lies outside $[-1, 1]$ or a field is not positive.

    Convention
    ----------
    Magnitude only: the sign is the direction of travel, which the caller
    tracks. $\mu = mv_\perp^2/(2B)$ and $E = mv^2/2$ are the invariants;
    a static field, no electric potential.

    Physical interpretation
    -----------------------
    The magnetic mirror: $v_\parallel$ falls as the particle moves into
    stronger field and vanishes where $B = B_\mathrm{ref}/(1 - \xi_\mathrm{ref}^2)$,
    the bounce point of a trapped particle.

    Assumptions
    -----------
    Adiabatic motion: $\mu$ conserved, field varying slowly over a gyro-orbit.

    See Also
    --------
    vaft.diagram.trapped_and_passing_orbits : the canonical diagram of this relation.

    References
    ----------
    .. [1] F. F. Chen, *Introduction to Plasma Physics and Controlled Fusion*,
           3rd ed., Springer (2016), Sec. 2.3.3.
    """
    xi = float(pitch_ref)
    B_ref = float(B_ref)
    B = np.asarray(B, dtype=float)
    if not -1.0 <= xi <= 1.0:
        raise ValueError(f"pitch_ref must lie in [-1, 1], not {pitch_ref!r}")
    if B_ref <= 0.0 or np.any(B <= 0.0):
        raise ValueError("fields must be positive")
    arg = 1.0 - (1.0 - xi * xi) * B / B_ref
    # a bounce point computed from its own mirror ratio lands within round-off of zero
    arg = np.where(np.abs(arg) < 1e-12, 0.0, arg)
    result = float(v) * np.sqrt(np.where(arg >= 0.0, arg, np.nan))
    return float(result) if np.ndim(result) == 0 else result

def psi_per_radian_from_cocos(psi, cocos):
    r"""The guiding-centre $\psi = RA_\phi$ from a flux stored in a COCOS convention.

    $$\psi_{RA_\phi} = -\,\sigma_{Bp}\,\sigma_{R\phi Z}\,\frac{\psi_\mathrm{COCOS}}{(2\pi)^{e_{Bp}}}$$

    Parameters
    ----------
    psi : float or np.ndarray
        Poloidal flux as the COCOS convention stores it [Wb or Wb/rad].
    cocos : int
        The COCOS index, 1-8 or 11-18 [-].

    Returns
    -------
    float or np.ndarray
        $RA_\phi$ with the IMAS $\phi$, the flux per radian that
        ``guiding_center_toroidal_momentum`` takes [Wb/rad].

    Raises
    ------
    ValueError
        ``cocos`` is not one of the sixteen indices.

    Convention
    ----------
    $\mathbf B_\mathrm{pol} = \nabla(RA_\phi)\times\nabla\phi$ with $\phi$
    counter-clockwise from above, so $RA_\phi$ is largest on the axis for a
    current along $+\phi$. Sauter's COCOS flux has $d\psi/d\rho$ of sign
    $\sigma_{Ip}\sigma_{Bp}$ and its own $\phi$ orientation $\sigma_{R\phi Z}$,
    hence the factor: $-1/(2\pi)$ for COCOS 11 (IMAS DD3), $+1/(2\pi)$ for
    COCOS 17 (DD4), $-1$ for COCOS 1. $e_{Bp} = 1$ for 11-18 (whole-turn flux in Wb).

    Physical interpretation
    -----------------------
    Only this signed, per-radian flux makes $q\psi + mv_\parallel Rb_\phi$ a
    conserved quantity; a COCOS flux used as stored gives the flux term the
    wrong sign or a $2\pi$ error.

    References
    ----------
    .. [1] O. Sauter and S. Yu. Medvedev, Comput. Phys. Commun. 184 (2013)
           293, Eq. 8 and Table I.
    """
    if isinstance(cocos, bool) or int(cocos) != cocos or int(cocos) not in (*range(1, 9), *range(11, 19)):
        raise ValueError(f"cocos must be one of 1-8 or 11-18, not {cocos!r}")
    c = int(cocos)
    base = c % 10
    sigma_bp = 1 if base in (1, 2, 5, 6) else -1
    sigma_rphiz = 1 if base % 2 == 1 else -1
    exponent = 1 if c >= 11 else 0
    result = -sigma_bp * sigma_rphiz * np.asarray(psi, dtype=float) / (2.0 * np.pi) ** exponent
    return float(result) if np.ndim(result) == 0 else result