Skip to main content

Documentation Index

Fetch the complete documentation index at: https://mintlify.com/stourgai/WPIT/llms.txt

Use this file to discover all available pages before exploring further.

The EMIC_ion_mod sub-module implements the relativistic equations of motion for a single ion (proton, helium, or oxygen) interacting with an obliquely propagating Electromagnetic Ion Cyclotron (EMIC) wave in a dipole magnetic field. All coupling is expressed through right-hand and left-hand polarised wave components decomposed via Bessel functions of the first kind. Compared with whistler_electron_mod, this module adds two unique evolution equations — dkpardt and dwcdt — that track the variation of the parallel wave number and the cyclotron frequency along the particle trajectory, which are needed for accurate nonlinear trapping analysis in an inhomogeneous field.
import WPIT.WPI_mod.EMIC_ion_mod as emic

# Or import individual functions
from WPIT.WPI_mod.EMIC_ion_mod import (
    wpi_params,
    dzdt, dlamdadt,
    dppardt, dpperdt,
    dEkdt, dgammadt,
    dalphadt, daeqdt,
    detadt,
    dwcdt,
    nonlinear_C0, nonlinear_C1m, nonlinear_C1p,
    nonlinear_H, nonlinear_S, nonlinear_theta,
)
# dkpardt is defined in dkpardt.py but is not re-exported via __init__.py;
# import it directly if needed:
from WPIT.WPI_mod.EMIC_ion_mod.dkpardt import dkpardt

wpi_params

Computes derived wave–particle coupling parameters required by all time-derivative functions. Call this once per integration step before evaluating any d*/dt function. The function decomposes the transverse wave fields into right- and left-hand circularly polarised components and forms the corresponding reference momenta and frequencies.
beta, BwR, BwL, EwR, EwL, pwR, pwL, wR, wL = emic.wpi_params(
    pper_arg, kper_arg, qi_arg, mi_arg,
    Bmag_arg, Exw_arg, Eyw_arg,
    Bxw_arg, Byw_arg, gamma_arg
)
pper_arg
float
Ion momentum component perpendicular to the background magnetic field B₀, in kg m s⁻¹.
kper_arg
float
Perpendicular component of the wave number vector, in rad m⁻¹.
qi_arg
float
Ion charge, in Coulombs. Use WPIT.Environment_mod.const.qi for protons.
mi_arg
float
Ion rest mass, in kg. Use const.mH, const.mHe, or const.mO for proton, helium, or oxygen respectively.
Bmag_arg
float
Magnitude of the background geomagnetic field, in T.
Exw_arg
float
x-component of the wave electric field amplitude, in V m⁻¹.
Eyw_arg
float
y-component of the wave electric field amplitude, in V m⁻¹.
Bxw_arg
float
x-component of the wave magnetic field amplitude, in T.
Byw_arg
float
y-component of the wave magnetic field amplitude, in T.
gamma_arg
float
Relativistic Lorentz factor of the ion (dimensionless).
Returns
beta
float
Bessel function argument beta = -(kper * pper) / (qi * Bmag) (dimensionless).
BwR
float
Right-hand circularly polarised wave magnetic field amplitude, BwR = 0.5*(Bxw + Byw), in T.
BwL
float
Left-hand circularly polarised wave magnetic field amplitude, BwL = 0.5*(Bxw - Byw), in T.
EwR
float
Right-hand circularly polarised wave electric field amplitude, in V m⁻¹.
EwL
float
Left-hand circularly polarised wave electric field amplitude, in V m⁻¹.
pwR
float
Right-hand reference momentum pwR = gamma * mi * (EwR / BwR), in kg m s⁻¹.
pwL
float
Left-hand reference momentum pwL = gamma * mi * (EwL / BwL), in kg m s⁻¹.
wR
float
Right-hand effective cyclotron-like frequency wR = (qi * BwR) / (gamma * mi), in rad s⁻¹.
wL
float
Left-hand effective cyclotron-like frequency wL = (qi * BwL) / (gamma * mi), in rad s⁻¹.

dzdt

Returns the time rate of change of the distance along the magnetic field line.
dz_dt = emic.dzdt(ppar_arg, gamma_arg, mi_arg)
ppar_arg
float
Ion momentum component parallel to B₀, in kg m s⁻¹.
gamma_arg
float
Lorentz factor (dimensionless).
mi_arg
float
Ion rest mass, in kg.
dz_dt
float
Time derivative of the field-line distance dz/dt = ppar / (mi * gamma), in m s⁻¹.

dlamdadt

Returns the time rate of change of the magnetic latitude lambda along the field line, accounting for dipole geometry.
dlamda_dt = emic.dlamdadt(ppar_arg, lamda_arg, gamma_arg, mi_arg, L_arg)
ppar_arg
float
Parallel ion momentum, in kg m s⁻¹.
lamda_arg
float
Magnetic latitude (lambda), in radians.
gamma_arg
float
Lorentz factor (dimensionless).
mi_arg
float
Ion rest mass, in kg.
L_arg
float
Magnetic L-shell value (dimensionless).
dlamda_dt
float
Time derivative of the magnetic latitude, in rad s⁻¹.

dppardt

Returns the time rate of change of the parallel momentum, combining the wave-driven force (expressed via Bessel functions) with the adiabatic mirror force.
dppar_dt = emic.dppardt(
    pper_arg, eta_arg, gamma_arg, m_res_arg,
    qi_arg, mi_arg, Ewz_arg, beta_arg,
    wR_arg, wL_arg, Bmag_arg, dBdz_arg
)
pper_arg
float
Perpendicular momentum, in kg m s⁻¹.
eta_arg
float
Wave–particle phase angle, in radians.
gamma_arg
float
Lorentz factor (dimensionless).
m_res_arg
int
Resonance harmonic order m (integer).
qi_arg
float
Ion charge, in Coulombs.
mi_arg
float
Ion rest mass, in kg.
Ewz_arg
float
Field-aligned (z) component of the wave electric field amplitude, in V m⁻¹.
beta_arg
float
Bessel argument beta from wpi_params (dimensionless).
wR_arg
float
Right-hand effective frequency from wpi_params, in rad s⁻¹.
wL_arg
float
Left-hand effective frequency from wpi_params, in rad s⁻¹.
Bmag_arg
float
Background magnetic field magnitude, in T.
dBdz_arg
float
Spatial gradient of the magnetic field magnitude along the field line, in T m⁻¹.
dppar_dt
float
Time derivative of the parallel momentum dppar/dt, in kg m s⁻².

dpperdt

Returns the time rate of change of the perpendicular momentum. The wave term involves the difference between parallel momentum and the reference momenta pwR and pwL.
dpper_dt = emic.dpperdt(
    pper_arg, ppar_arg, eta_arg, gamma_arg, m_res_arg,
    qi_arg, mi_arg, pwR_arg, pwL_arg,
    beta_arg, wR_arg, wL_arg, Bmag_arg, dBdz_arg
)
pper_arg
float
Perpendicular momentum, in kg m s⁻¹.
ppar_arg
float
Parallel momentum, in kg m s⁻¹.
eta_arg
float
Wave–particle phase angle, in radians.
gamma_arg
float
Lorentz factor (dimensionless).
m_res_arg
int
Resonance harmonic order m.
qi_arg
float
Ion charge, in Coulombs.
mi_arg
float
Ion rest mass, in kg.
pwR_arg
float
Right-hand reference momentum from wpi_params, in kg m s⁻¹.
pwL_arg
float
Left-hand reference momentum from wpi_params, in kg m s⁻¹.
beta_arg
float
Bessel argument beta (dimensionless).
wR_arg
float
Right-hand effective frequency, in rad s⁻¹.
wL_arg
float
Left-hand effective frequency, in rad s⁻¹.
Bmag_arg
float
Background magnetic field magnitude, in T.
dBdz_arg
float
Spatial gradient of the magnetic field along the field line, in T m⁻¹.
dpper_dt
float
Time derivative of the perpendicular momentum dpper/dt, in kg m s⁻².

dEkdt

Returns the time rate of change of the relativistic kinetic energy of the ion due to wave–particle interaction.
dEk_dt = emic.dEkdt(
    pper_arg, ppar_arg, eta_arg, gamma_arg, m_res_arg,
    qi_arg, mi_arg, Ewz_arg, beta_arg,
    EwR_arg, EwL_arg, wmega_arg, kappa_arg
)
pper_arg
float
Perpendicular momentum, in kg m s⁻¹.
ppar_arg
float
Parallel momentum, in kg m s⁻¹.
eta_arg
float
Wave–particle phase angle, in radians.
gamma_arg
float
Lorentz factor (dimensionless).
m_res_arg
int
Resonance harmonic order m.
qi_arg
float
Ion charge, in Coulombs.
mi_arg
float
Ion rest mass, in kg.
Ewz_arg
float
Field-aligned wave electric field component, in V m⁻¹.
beta_arg
float
Bessel argument beta (dimensionless).
EwR_arg
float
Right-hand wave electric field amplitude, in V m⁻¹.
EwL_arg
float
Left-hand wave electric field amplitude, in V m⁻¹.
wmega_arg
float
Wave angular frequency, in rad s⁻¹.
kappa_arg
float
Wave number magnitude, in rad m⁻¹.
dEk_dt
float
Time derivative of the kinetic energy, in J s⁻¹.

dgammadt

Returns the time rate of change of the Lorentz factor due to the wave–particle interaction.
dgamma_dt = emic.dgammadt(
    pper_arg, ppar_arg, eta_arg, gamma_arg, m_res_arg,
    qi_arg, mi_arg, Ewz_arg, beta_arg, EwR_arg, EwL_arg
)
pper_arg
float
Perpendicular momentum, in kg m s⁻¹.
ppar_arg
float
Parallel momentum, in kg m s⁻¹.
eta_arg
float
Wave–particle phase angle, in radians.
gamma_arg
float
Lorentz factor (dimensionless).
m_res_arg
int
Resonance harmonic order m.
qi_arg
float
Ion charge, in Coulombs.
mi_arg
float
Ion rest mass, in kg.
Ewz_arg
float
Field-aligned wave electric field component, in V m⁻¹.
beta_arg
float
Bessel argument beta (dimensionless).
EwR_arg
float
Right-hand wave electric field amplitude, in V m⁻¹.
EwL_arg
float
Left-hand wave electric field amplitude, in V m⁻¹.
dgamma_dt
float
Time derivative of the Lorentz factor, in s⁻¹.

dalphadt

Returns the time rate of change of the local pitch angle.
dalpha_dt = emic.dalphadt(
    pper_arg, ppar_arg, eta_arg, Ewz_arg, m_res_arg,
    qi_arg, pwR_arg, pwL_arg, beta_arg, wR_arg, wL_arg
)
pper_arg
float
Perpendicular momentum, in kg m s⁻¹.
ppar_arg
float
Parallel momentum, in kg m s⁻¹.
eta_arg
float
Wave–particle phase angle, in radians.
Ewz_arg
float
Field-aligned wave electric field component, in V m⁻¹.
m_res_arg
int
Resonance harmonic order m.
qi_arg
float
Ion charge, in Coulombs.
pwR_arg
float
Right-hand reference momentum from wpi_params, in kg m s⁻¹.
pwL_arg
float
Left-hand reference momentum from wpi_params, in kg m s⁻¹.
beta_arg
float
Bessel argument beta (dimensionless).
wR_arg
float
Right-hand effective frequency, in rad s⁻¹.
wL_arg
float
Left-hand effective frequency, in rad s⁻¹.
dalpha_dt
float
Time derivative of the local pitch angle, in rad s⁻¹.

daeqdt

Returns the time rate of change of the equatorial pitch angle, combining wave forcing and the tan(aeq) geometric factor that maps local pitch angle changes to their equatorial equivalent.
daeq_dt = emic.daeqdt(
    pper_arg, ppar_arg, eta_arg, aeq_arg,
    Ewz_arg, gamma_arg, m_res_arg,
    qi_arg, mi_arg, pwR_arg, pwL_arg,
    beta_arg, wR_arg, wL_arg
)
pper_arg
float
Perpendicular momentum, in kg m s⁻¹.
ppar_arg
float
Parallel momentum, in kg m s⁻¹.
eta_arg
float
Wave–particle phase angle, in radians.
aeq_arg
float
Equatorial pitch angle, in radians.
Ewz_arg
float
Field-aligned wave electric field component, in V m⁻¹.
gamma_arg
float
Lorentz factor (dimensionless).
m_res_arg
int
Resonance harmonic order m.
qi_arg
float
Ion charge, in Coulombs.
mi_arg
float
Ion rest mass, in kg.
pwR_arg
float
Right-hand reference momentum, in kg m s⁻¹.
pwL_arg
float
Left-hand reference momentum, in kg m s⁻¹.
beta_arg
float
Bessel argument beta (dimensionless).
wR_arg
float
Right-hand effective frequency, in rad s⁻¹.
wL_arg
float
Left-hand effective frequency, in rad s⁻¹.
daeq_dt
float
Time derivative of the equatorial pitch angle, in rad s⁻¹.

detadt

Returns the time rate of change of the wave–particle phase angle eta, which governs whether the particle is trapped or passing relative to the wave potential.
deta_dt = emic.detadt(ppar_arg, mres_arg, wc_arg, gamma_arg, kpar_arg, mi_arg, w_arg)
ppar_arg
float
Parallel momentum, in kg m s⁻¹.
mres_arg
int
Resonance harmonic order m.
wc_arg
float
Ion cyclotron frequency magnitude, in rad s⁻¹.
gamma_arg
float
Lorentz factor (dimensionless).
kpar_arg
float
Parallel wave number component, in rad m⁻¹.
mi_arg
float
Ion rest mass, in kg.
w_arg
float
Wave angular frequency, in rad s⁻¹.
deta_dt
float
Time derivative of the wave–particle phase: deta/dt = (m*wc/gamma) + (kpar*ppar)/(gamma*mi) - w, in rad s⁻¹.

dkpardt

Unique to this module. Returns the rate of change of the parallel wave number experienced by the ion as it moves along the magnetic field line. This is needed to compute the nonlinear parameter H accurately in an inhomogeneous field.
from WPIT.WPI_mod.EMIC_ion_mod.dkpardt import dkpardt
dkpar_dt = dkpardt(ppar, m, gamma, psi, dkdz)
ppar
float
Parallel ion momentum, in kg m s⁻¹.
m
float
Ion rest mass, in kg.
gamma
float
Lorentz factor (dimensionless).
psi
float
Wave–particle phase angle (used as the cosine projection factor), in radians.
dkdz
float
Spatial gradient of the parallel wave number along the field line, in rad m⁻².
dkpar_dt
float
Time derivative of the parallel wave number: dkpar/dt = (ppar / (m*gamma)) * dkdz * cos(psi), in rad m⁻¹ s⁻¹.
dkpardt is defined in EMIC_ion_mod/dkpardt.py but is not re-exported through __init__.py. Import it directly with from WPIT.WPI_mod.EMIC_ion_mod.dkpardt import dkpardt. There is no equivalent in whistler_electron_mod; the EMIC ion interaction requires explicit tracking of parallel wave-number evolution because the EMIC dispersion surface varies significantly with field-line position and ion composition.

dwcdt

Unique to this module. Returns the rate of change of the ion cyclotron frequency along the particle trajectory as it moves in the inhomogeneous dipole field.
dwc_dt = emic.dwcdt(ppar, m, gamma, dwcdz)
ppar
float
Parallel ion momentum, in kg m s⁻¹.
m
float
Ion rest mass, in kg.
gamma
float
Lorentz factor (dimensionless).
dwcdz
float
Spatial gradient of the ion cyclotron frequency along the field line, in rad m⁻¹ s⁻¹.
dwc_dt
float
Time derivative of the cyclotron frequency: dwc/dt = (ppar / (m*gamma)) * dwcdz, in rad s⁻².

Nonlinear trapping diagnostics

The four nonlinear functions compute the trapping-island Hamiltonian terms following the formalism of Omura et al. (2010) adapted for EMIC waves with Bessel-function corrections for oblique propagation.

nonlinear_C0

Computes the inhomogeneity coefficient C0, the field-aligned component of the nonlinear force.
C0 = emic.nonlinear_C0(
    ppar_arg, kpar_arg, mres_arg, gamma_arg,
    qi_arg, mi_arg, wce_arg, Ewz_arg
)
ppar_arg
float
Parallel momentum, in kg m s⁻¹.
kpar_arg
float
Parallel wave number, in rad m⁻¹.
mres_arg
int
Resonance harmonic order m.
gamma_arg
float
Lorentz factor (dimensionless).
qi_arg
float
Ion charge, in Coulombs.
mi_arg
float
Ion rest mass, in kg.
wce_arg
float
Ion cyclotron frequency, in rad s⁻¹.
Ewz_arg
float
Field-aligned wave electric field amplitude, in V m⁻¹.
C0
float
The C0 nonlinear coefficient (SI units).

nonlinear_C1m

Computes the left-hand coupling coefficient C1m for the nonlinear trapping force.
C1m = emic.nonlinear_C1m(
    pper_arg, ppar_arg, kpar_arg, mres_arg,
    qi_arg, mi_arg, gamma_arg,
    wL_arg, EwL_arg, wce_arg
)
pper_arg
float
Perpendicular momentum, in kg m s⁻¹.
ppar_arg
float
Parallel momentum, in kg m s⁻¹.
kpar_arg
float
Parallel wave number, in rad m⁻¹.
mres_arg
int
Resonance harmonic order m.
qi_arg
float
Ion charge, in Coulombs.
mi_arg
float
Ion rest mass, in kg.
gamma_arg
float
Lorentz factor (dimensionless).
wL_arg
float
Left-hand effective frequency from wpi_params, in rad s⁻¹.
EwL_arg
float
Left-hand wave electric field amplitude, in V m⁻¹.
wce_arg
float
Ion cyclotron frequency, in rad s⁻¹.
C1m
float
The C1m left-hand nonlinear coupling coefficient (SI units).

nonlinear_C1p

Computes the right-hand coupling coefficient C1p for the nonlinear trapping force.
C1p = emic.nonlinear_C1p(
    pper_arg, ppar_arg, kpar_arg, mres_arg,
    qi_arg, mi_arg, gamma_arg,
    wR_arg, EwR_arg, wce_arg
)
pper_arg
float
Perpendicular momentum, in kg m s⁻¹.
ppar_arg
float
Parallel momentum, in kg m s⁻¹.
kpar_arg
float
Parallel wave number, in rad m⁻¹.
mres_arg
int
Resonance harmonic order m.
qi_arg
float
Ion charge, in Coulombs.
mi_arg
float
Ion rest mass, in kg.
gamma_arg
float
Lorentz factor (dimensionless).
wR_arg
float
Right-hand effective frequency from wpi_params, in rad s⁻¹.
EwR_arg
float
Right-hand wave electric field amplitude, in V m⁻¹.
wce_arg
float
Ion cyclotron frequency, in rad s⁻¹.
C1p
float
The C1p right-hand nonlinear coupling coefficient (SI units).

nonlinear_H

Computes the nonlinear trapping Hamiltonian H, which measures the distance of a particle from exact resonance including field inhomogeneity corrections. When H/wtsq (the trapping parameter S) satisfies |S| < 1, the particle is phase-trapped.
H = emic.nonlinear_H(
    pper_arg, ppar_arg, kpar_arg, gamma_arg,
    mres_arg, mi_arg, wce_arg,
    dkpar_dt_arg, dwcdz_arg, dwdt_arg
)
pper_arg
float
Perpendicular momentum, in kg m s⁻¹.
ppar_arg
float
Parallel momentum, in kg m s⁻¹.
kpar_arg
float
Parallel wave number, in rad m⁻¹.
gamma_arg
float
Lorentz factor (dimensionless).
mres_arg
int
Resonance harmonic order m.
mi_arg
float
Ion rest mass, in kg.
wce_arg
float
Ion cyclotron frequency, in rad s⁻¹.
dkpar_dt_arg
float
Time derivative of the parallel wave number (from dkpardt), in rad m⁻¹ s⁻¹.
dwcdz_arg
float
Spatial gradient of the cyclotron frequency along the field line, in rad m⁻¹ s⁻¹.
dwdt_arg
float
Time derivative of the wave frequency (typically zero for monochromatic waves), in rad s⁻².
H
float
The nonlinear Hamiltonian H. Pass to nonlinear_S together with the trapping frequency squared.

nonlinear_S

Computes the dimensionless trapping parameter S = H / omega_t_sq. Phase trapping occurs when |S| < 1.
S = emic.nonlinear_S(H_arg, wtsq_arg)
H_arg
float
The nonlinear Hamiltonian H from nonlinear_H.
wtsq_arg
float
The trapping frequency squared from nonlinear_theta (in rad² s⁻²).
S
float
Dimensionless nonlinear trapping parameter. |S| less than 1 indicates phase trapping; |S| greater than 1 indicates phase liberation.

nonlinear_theta

Computes the nonlinear coupling term theta and the trapping frequency squared, using Bessel functions evaluated at beta to account for the oblique propagation geometry.
theta, wtsq = emic.nonlinear_theta(C0_arg, Cp1_arg, Cm1_arg, m_res_arg, beta_arg)
C0_arg
float
Field-aligned coefficient from nonlinear_C0.
Cp1_arg
float
Right-hand coefficient from nonlinear_C1p.
Cm1_arg
float
Left-hand coefficient from nonlinear_C1m.
m_res_arg
int
Resonance harmonic order m.
beta_arg
float
Bessel argument beta from wpi_params (dimensionless).
theta
float
Combined nonlinear trapping coupling term (Bessel-weighted sum of C0, C1p, C1m).
wtsq
float
Trapping frequency squared |theta| (in rad² s⁻²), passed to nonlinear_S as wtsq_arg.

Complete integration step example

The following example shows how to assemble one RK4-style step of the EMIC–ion equations of motion, mimicking how the functions are chained together in practice.
import numpy as np
import WPIT.WPI_mod.EMIC_ion_mod as emic
from WPIT.WPI_mod.EMIC_ion_mod.dkpardt import dkpardt as emic_dkpardt
from WPIT.Environment_mod import const

# --- Ion species: proton ---
mi = const.mH        # proton mass [kg]
qi = const.qi        # proton charge [C]

# --- Simulation state at time t ---
ppar = 1.2e-19       # parallel momentum [kg m/s]
pper = 3.5e-19       # perpendicular momentum [kg m/s]
lamda = np.deg2rad(5.0)   # magnetic latitude [rad]
eta   = 0.3          # wave-particle phase [rad]
aeq   = np.deg2rad(60.0)  # equatorial pitch angle [rad]
L_shell = 4.5        # L-shell

# --- Background field quantities ---
Bmag  = 150e-9       # |B0| in T (illustrative)
dBdz  = 2e-12        # dB/dz [T/m]
dwcdz = 5.0          # d(wc)/dz [rad/m/s]
dkdz  = 0.01         # dk_par/dz [rad/m^2]

# --- Wave parameters ---
Exw, Eyw, Ewz = 1e-4, 1e-4, 0.5e-4   # wave E-field [V/m]
Bxw, Byw      = 5e-9, 5e-9            # wave B-field [T]
kpar  = 3.0e-5       # parallel wave number [rad/m]
kperp = 1.5e-5       # perpendicular wave number [rad/m]
w_wave = 0.3 * qi * Bmag / mi         # wave frequency [rad/s]

# --- Lorentz factor ---
p_tot = np.sqrt(ppar**2 + pper**2)
gamma = np.sqrt(1 + (p_tot / (mi * const.c_light))**2)

# --- Resonance order ---
m_res = 1

# Step 1: compute coupling parameters
beta, BwR, BwL, EwR, EwL, pwR, pwL, wR, wL = emic.wpi_params(
    pper, kperp, qi, mi, Bmag, Exw, Eyw, Bxw, Byw, gamma
)

# Step 2: evaluate equations of motion
dz    = emic.dzdt(ppar, gamma, mi)
dlamda = emic.dlamdadt(ppar, lamda, gamma, mi, L_shell)
dppar = emic.dppardt(pper, eta, gamma, m_res, qi, mi, Ewz,
                     beta, wR, wL, Bmag, dBdz)
dpper = emic.dpperdt(pper, ppar, eta, gamma, m_res, qi, mi,
                     pwR, pwL, beta, wR, wL, Bmag, dBdz)
dEk   = emic.dEkdt(pper, ppar, eta, gamma, m_res, qi, mi,
                   Ewz, beta, EwR, EwL, w_wave, kpar)
dgam  = emic.dgammadt(pper, ppar, eta, gamma, m_res, qi, mi,
                      Ewz, beta, EwR, EwL)
dalpha = emic.dalphadt(pper, ppar, eta, Ewz, m_res, qi, pwR, pwL,
                       beta, wR, wL)
wc = qi * Bmag / (gamma * mi)
daeq  = emic.daeqdt(pper, ppar, eta, aeq, Ewz, gamma, m_res, qi, mi,
                    pwR, pwL, beta, wR, wL)
deta  = emic.detadt(ppar, m_res, wc, gamma, kpar, mi, w_wave)
dkpar = emic_dkpardt(ppar, mi, gamma, eta, dkdz)
dwc   = emic.dwcdt(ppar, mi, gamma, dwcdz)

# Step 3: nonlinear trapping diagnostics
C0  = emic.nonlinear_C0(ppar, kpar, m_res, gamma, qi, mi, wc, Ewz)
C1m = emic.nonlinear_C1m(pper, ppar, kpar, m_res, qi, mi, gamma, wL, EwL, wc)
C1p = emic.nonlinear_C1p(pper, ppar, kpar, m_res, qi, mi, gamma, wR, EwR, wc)
theta_val, wtsq = emic.nonlinear_theta(C0, C1p, C1m, m_res, beta)
H = emic.nonlinear_H(pper, ppar, kpar, gamma, m_res, mi, wc,
                     dkpar, dwcdz, 0.0)
S = emic.nonlinear_S(H, wtsq)

print(f"Trapping parameter S = {S:.4f}")
print(f"  |S| < 1 → trapped" if abs(S) < 1 else "  |S| >= 1 → passing/librating")

Build docs developers (and LLMs) love