# plasma_plots.theory.orbits

*module*

Charged-particle orbits: gyromotion, guiding-center drifts and trapped particles in a tokamak.

The formulas hold in any consistent units (SI, or Struphy's normalized units with charge and
mass in units of e and m_p) unless a function says otherwise. Vectors are arrays with a trailing
axis of length 3, (x, y, z) or any right-handed Cartesian basis; everything broadcasts.

> **Tokamak model**
>
> The trapped-particle functions use a large-aspect-ratio tokamak with circular flux surfaces,
> B(θ) = B₀ / (1 + ε cos θ) with ε = r/R₀ the inverse aspect ratio of the flux surface and θ the
> poloidal angle (θ = 0 on the outboard midplane, where B is smallest). A particle of speed v and
> magnetic moment μ has the pitch parameter λ = μB₀/E = (v⊥²/v²) B₀/B, and
> 
>     v∥² = v² (1 + ε cos θ − λ) / (1 + ε cos θ) = 2ε v² (κ² − sin²(θ/2)) / (1 + ε cos θ)
> 
> with the trapping parameter κ² = (1 + ε − λ)/(2ε): particles with κ² < 1 are trapped (they bounce
> at sin²(θ/2) = κ²), κ² = 0 is deeply trapped at θ = 0, κ² > 1 passes. The bounce and transit
> frequencies and banana widths are to leading order in ε (the factor 1 + ε cos θ is dropped,
> and the field line length is q R₀ dθ).

> **References**
>
> P. Helander and D. J. Sigmar, Collisional Transport in Magnetized Plasmas (Cambridge, 2002),
> chapter 7.
> 
> J. Wesson, Tokamaks, 4th ed. (Oxford, 2011), section 3.12.
> 
> F. F. Chen, Introduction to Plasma Physics and Controlled Fusion, 3rd ed. (Springer, 2016),
> chapter 2.

**Examples**

The bounce frequency of a barely trapped (κ² = 0.9) particle with v = 1 on the ε = 0.1 surface of
a tokamak with q = 2 and R₀ = 3:

```pycon
>>> print(f"{bounce_frequency(1.0, 0.9, 0.1, 2.0, 3.0):.5f}")
0.02271
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/orbits.py#L1-L1)

## plasma_plots.theory.orbits.banana_width

*function*

```python
def banana_width(gyroradius, kappa2, epsilon, safety_factor)
```

Compute the full radial width of a banana orbit, Δr = 2√2 q κ ρ / √ε.

From the conservation of the canonical toroidal momentum, the radial excursion is
δr = q δv∥ / (ε Ω); v∥ swings between ±κ v √(2ε) at the outboard midplane. The width is
largest near the trapped–passing boundary, 2√2 q ρ / √ε, the "q ρ / √ε" of scaling estimates,
and vanishes for deeply trapped particles.

**Parameters**

- `gyroradius` (`float or array_like`) — Gyroradius ρ = m v / (|q| B₀) with the total speed v (for trapped particles v⊥ ≈ v).
- `kappa2` (`float or array_like`) — Trapping parameter κ², 0 ≤ κ² ≤ 1.
- `epsilon` (`float or array_like`) — Inverse aspect ratio ε of the flux surface.
- `safety_factor` (`float or array_like`) — Safety factor q of the flux surface.

**Returns**

- (`float or numpy.ndarray`) — Δr, in the units of ``gyroradius``; ``nan`` for passing particles (κ² > 1).

> **References**
>
> Helander and Sigmar, Collisional Transport in Magnetized Plasmas, section 7.3.
> 
> Wesson, Tokamaks, section 3.12.

**Examples**

```pycon
>>> print(f"{banana_width(0.01, 1.0, 0.1, 2.0):.4f}")
0.1789
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/orbits.py#L544-L584)

## plasma_plots.theory.orbits.bounce_frequency

*function*

```python
def bounce_frequency(speed, kappa2, epsilon, safety_factor, major_radius)
```

Compute the bounce frequency of a trapped particle, ω_b = π v √(2ε) / (4 q R₀ K(κ²)).

The period of the banana orbit is τ_b = 2π/ω_b = 8 q R₀ K(κ²) / (v √(2ε)), to leading order
in ε. Deeply trapped particles (κ² = 0) bounce at ω_b = v √(ε/2) / (q R₀); ω_b → 0
logarithmically at the trapped–passing boundary κ² → 1.

**Parameters**

- `speed` (`float or array_like`) — The particle speed v.
- `kappa2` (`float or array_like`) — Trapping parameter κ² (see [`trapping_parameter()`][plasma_plots.theory.orbits.trapping_parameter]), 0 ≤ κ² < 1.
- `epsilon` (`float or array_like`) — Inverse aspect ratio ε of the flux surface.
- `safety_factor` (`float or array_like`) — Safety factor q of the flux surface.
- `major_radius` (`float or array_like`) — Major radius R₀.

**Returns**

- (`float or numpy.ndarray`) — ω_b in rad per unit time; ``nan`` for passing particles (κ² > 1).

> **References**
>
> Helander and Sigmar, Collisional Transport in Magnetized Plasmas, section 7.2.
> 
> Wesson, Tokamaks, section 3.12 (the deeply trapped limit).

**Examples**

```pycon
>>> # √(ε/2) deeply trapped
>>> print(f"{bounce_frequency(1.0, 0.0, 0.02, 1.0, 1.0):.4f}")
0.1000
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/orbits.py#L442-L490)

## plasma_plots.theory.orbits.exb_drift

*function*

```python
def exb_drift(E, B)
```

Compute the E×B drift velocity v_E = E × B / B².

**Parameters**

- `E` (`array_like`) — Electric field, shape (..., 3).
- `B` (`array_like`) — Magnetic field, shape (..., 3).

**Returns**

- (`numpy.ndarray`) — v_E, shape (..., 3); independent of charge and mass, |v_E| = E⊥/B.

**Raises**

- `ValueError` — If a vector has no trailing axis of length 3.

> **References**
>
> Chen, Introduction to Plasma Physics, section 2.2.2.

**Examples**

```pycon
>>> print(exb_drift([1.0, 0.0, 0.0], [0.0, 0.0, 2.0]))
[ 0.  -0.5  0. ]
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/orbits.py#L168-L198)

## plasma_plots.theory.orbits.grad_b_drift

*function*

```python
def grad_b_drift(perpendicular_speed, B, grad_B, charge=1.0, mass=1.0)
```

Compute the grad-B drift velocity v_∇B = (m v⊥² / (2q)) B × ∇B / B³.

**Parameters**

- `perpendicular_speed` (`float or array_like`) — The speed v⊥ perpendicular to B.
- `B` (`array_like`) — Magnetic field, shape (..., 3).
- `grad_B` (`array_like`) — Gradient of the field strength ∇|B|, shape (..., 3).
- `charge` (`float or array_like`) (default: `1.0`) — Charge q; ions and electrons drift in opposite directions. Default: ``1``.
- `mass` (`float or array_like`) (default: `1.0`) — Mass m. Default: ``1``.

**Returns**

- (`numpy.ndarray`) — v_∇B, shape (..., 3).

**Raises**

- `ValueError` — If a vector has no trailing axis of length 3.

> **References**
>
> Chen, Introduction to Plasma Physics, section 2.3.1.

**Examples**

```pycon
>>> # B along z, ∇B along x
>>> print(grad_b_drift(1.0, [0.0, 0.0, 1.0], [0.1, 0.0, 0.0]))
[0.   0.05 0.  ]
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/orbits.py#L201-L240)

## plasma_plots.theory.orbits.gyrofrequency

*function*

```python
def gyrofrequency(field, charge=1.0, mass=1.0)
```

Compute the gyrofrequency Ω = q|B|/m.

**Parameters**

- `field` (`float or array_like`) — Magnetic field strength |B| (its sign doesn't matter).
- `charge` (`float or array_like`) (default: `1.0`) — Charge q; Ω has its sign (negative for electrons). Default: ``1``.
- `mass` (`float or array_like`) (default: `1.0`) — Mass m. Default: ``1``.

**Returns**

- (`float or numpy.ndarray`) — Ω, in rad per unit time (a particle with q > 0 gyrates clockwise when viewed along B).

> **References**
>
> Chen, Introduction to Plasma Physics, section 2.2.

**Examples**

```pycon
>>> print(gyrofrequency(2.0, charge=-1.0, mass=0.5))
-4.0
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/orbits.py#L68-L96)

## plasma_plots.theory.orbits.gyroradius

*function*

```python
def gyroradius(perpendicular_speed, field, charge=1.0, mass=1.0)
```

Compute the gyroradius ρ = m v⊥ / (|q| B).

**Parameters**

- `perpendicular_speed` (`float or array_like`) — The speed v⊥ perpendicular to B.
- `field` (`float or array_like`) — Magnetic field strength |B|.
- `charge` (`float or array_like`) (default: `1.0`) — Charge q (its sign doesn't matter). Default: ``1``.
- `mass` (`float or array_like`) (default: `1.0`) — Mass m. Default: ``1``.

**Returns**

- (`float or numpy.ndarray`) — ρ (positive).

> **References**
>
> Chen, Introduction to Plasma Physics, section 2.2.

**Examples**

```pycon
>>> print(gyroradius(3.0, 2.0))
1.5
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/orbits.py#L99-L132)

## plasma_plots.theory.orbits.magnetic_moment

*function*

```python
def magnetic_moment(perpendicular_speed, field, mass=1.0)
```

Compute the magnetic moment μ = m v⊥² / (2B), the adiabatic invariant of the gyromotion.

**Parameters**

- `perpendicular_speed` (`float or array_like`) — The speed v⊥ perpendicular to B.
- `field` (`float or array_like`) — Magnetic field strength |B|.
- `mass` (`float or array_like`) (default: `1.0`) — Mass m. Default: ``1``.

**Returns**

- (`float or numpy.ndarray`) — μ.

> **References**
>
> Chen, Introduction to Plasma Physics, section 2.3.

**Examples**

```pycon
>>> print(magnetic_moment(2.0, 4.0, mass=2.0))
1.0
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/orbits.py#L135-L162)

## plasma_plots.theory.orbits.pitch_parameter

*function*

```python
def pitch_parameter(pitch, epsilon, theta=0.0)
```

Compute the pitch parameter λ = μB₀/E from the pitch v∥/v at the poloidal angle θ.

λ = (1 − (v∥/v)²) B₀/B(θ) = (1 − (v∥/v)²)(1 + ε cos θ).

**Parameters**

- `pitch` (`float or array_like`) — v∥/v, between −1 and 1, at the angle ``theta``.
- `epsilon` (`float or array_like`) — Inverse aspect ratio ε of the flux surface.
- `theta` (`float or array_like`) (default: `0.0`) — Poloidal angle θ where the pitch is given (0: outboard midplane). Default: ``0``.

**Returns**

- (`float or numpy.ndarray`) — λ.

> **References**
>
> Helander and Sigmar, Collisional Transport in Magnetized Plasmas, section 7.1.

**Examples**

```pycon
>>> print(f"{pitch_parameter(0.3, 0.1):.4f}")
1.0010
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/orbits.py#L329-L359)

## plasma_plots.theory.orbits.transit_frequency

*function*

```python
def transit_frequency(speed, kappa2, epsilon, safety_factor, major_radius)
```

Compute the poloidal transit frequency of a passing particle.

ω_t = π κ v √(2ε) / (2 q R₀ K(1/κ²)), to leading order in ε: 2π/ω_t is the time for one
poloidal turn. Far from the boundary (κ² ≫ 1) it tends to |v∥| / (q R₀); ω_t → 0
logarithmically at κ² → 1.

**Parameters**

- `speed` (`float or array_like`) — The particle speed v.
- `kappa2` (`float or array_like`) — Trapping parameter κ² (see [`trapping_parameter()`][plasma_plots.theory.orbits.trapping_parameter]), κ² > 1.
- `epsilon` (`float or array_like`) — Inverse aspect ratio ε of the flux surface.
- `safety_factor` (`float or array_like`) — Safety factor q of the flux surface.
- `major_radius` (`float or array_like`) — Major radius R₀.

**Returns**

- (`float or numpy.ndarray`) — ω_t in rad per unit time; ``nan`` for trapped particles (κ² < 1).

> **References**
>
> Helander and Sigmar, Collisional Transport in Magnetized Plasmas, section 7.2.

**Examples**

```pycon
>>> # ≈ v∥/(qR₀) = 1
>>> print(f"{transit_frequency(1.0, 50.0, 0.01, 1.0, 1.0):.4f}")
0.9950
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/orbits.py#L493-L541)

## plasma_plots.theory.orbits.trapped_fraction

*function*

```python
def trapped_fraction(epsilon, approximation='exact')
```

Compute the effective fraction of trapped particles on a flux surface of a circular tokamak.

f_t = 1 − (3/4) ⟨B²⟩ ∫₀^{1/B_max} λ dλ / ⟨√(1 − λB)⟩ with B = B₀/(1 + ε cos θ) and the
flux-surface average ⟨A⟩ = ∮ A (1 + ε cos θ) dθ / 2π, the fraction that enters neoclassical
transport (bootstrap current, neoclassical resistivity).

**Parameters**

- `epsilon` (`float or array_like`) — Inverse aspect ratio ε = r/R₀ of the flux surface, 0 ≤ ε < 1.
- `approximation` (`('exact', 'lin-liu', 'sqrt')`) (default: `"exact"`) — ``"exact"``: the integral above, by quadrature (relative accuracy 1e-8 or better for ε ≥ 1e-8); ``"lin-liu"``: Lin-Liu and Miller's fit 1 − (1 − ε)² / (√(1 − ε²) (1 + 1.46 √ε)), within a few per cent at any ε; ``"sqrt"``: the small-ε limit 1.46 √ε. Default: ``"exact"``.

**Returns**

- (`float or numpy.ndarray`) — f_t, 0 at ε = 0 and 1 at ε = 1.

**Raises**

- `ValueError` — If ``approximation`` is unknown.

> **References**
>
> Helander and Sigmar, Collisional Transport in Magnetized Plasmas, section 11.2 (f_t ≈ 1.46 √ε).
> 
> Y. R. Lin-Liu and R. L. Miller, "Upper and lower bounds of the effective trapped particle
> fraction in general tokamak equilibria", Phys. Plasmas 2, 1666 (1995).

**Examples**

```pycon
>>> fractions = [
...     trapped_fraction(0.1, a) for a in ("exact", "lin-liu", "sqrt")
... ]
>>> print(", ".join(f"{f:.4f}" for f in fractions))
0.4492, 0.4431, 0.4617
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/orbits.py#L270-L326)

## plasma_plots.theory.orbits.trapping_boundary

*function*

```python
def trapping_boundary(epsilon, quantity='lambda', theta=0.0)
```

Compute the trapped–passing boundary on a flux surface.

Particles with λ = μB₀/E > λ_c = B₀/B_max = 1 − ε are trapped, equivalently those with a
pitch |v∥/v| < √(1 − (1 − ε)/(1 + ε cos θ)) at the angle θ (√(2ε/(1 + ε)) on the outboard
midplane).

**Parameters**

- `epsilon` (`float or array_like`) — Inverse aspect ratio ε of the flux surface.
- `quantity` (`('lambda', 'pitch')`) (default: `"lambda"`) — ``"lambda"``: the critical λ_c; ``"pitch"``: the critical |v∥/v| at ``theta``. Default: ``"lambda"``.
- `theta` (`float or array_like`) (default: `0.0`) — Poloidal angle for ``quantity="pitch"``. Default: ``0``.

**Returns**

- (`float or numpy.ndarray`) — λ_c or the critical pitch.

**Raises**

- `ValueError` — If ``quantity`` is unknown.

> **References**
>
> Wesson, Tokamaks, section 3.12.

**Examples**

```pycon
>>> print(
...     f"{trapping_boundary(0.1):.2f}, {trapping_boundary(0.1, 'pitch'):.4f}"
... )
0.90, 0.4264
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/orbits.py#L362-L406)

## plasma_plots.theory.orbits.trapping_parameter

*function*

```python
def trapping_parameter(lam, epsilon)
```

Compute the trapping parameter κ² = (1 + ε − λ) / (2ε).

κ² < 1 for trapped particles, which bounce at sin²(θ/2) = κ²; κ² = 0 for particles deeply
trapped at θ = 0 (λ = 1 + ε), κ² = 1 on the trapped–passing boundary (λ = 1 − ε) and κ² > 1
for passing ones.

**Parameters**

- `lam` (`float or array_like`) — Pitch parameter λ = μB₀/E (see [`pitch_parameter()`][plasma_plots.theory.orbits.pitch_parameter]).
- `epsilon` (`float or array_like`) — Inverse aspect ratio ε > 0 of the flux surface.

**Returns**

- (`float or numpy.ndarray`) — κ².

> **References**
>
> Helander and Sigmar, Collisional Transport in Magnetized Plasmas, section 7.2.

**Examples**

```pycon
>>> print(f"{trapping_parameter(1.0, 0.1):.2f}")
0.50
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/orbits.py#L409-L439)