Operations and statistical analysis

vaft.formula is the physics layer of VAFT: closed-form expressions, empirical scalings and Green’s functions, written as pure functions of NumPy arrays and scalars. How it sits beside vaft.process and vaft.code is set out on Computational layers. Nothing in it touches an ODS, reads a file, or plots. To evaluate the same physics directly on a VEST shot, use the ODS-aware wrappers in vaft.omas.formula_wrapper (the vaft.omas.compute_* functions), which pull the inputs out of the IDSs and hand them to these functions.

flowchart LR
    A["ODS (IMAS IDSs)"] --> B["vaft.omas.formula_wrapper<br/>compute_tau_E_scaling, compute_power_balance, ..."]
    B --> C["vaft.formula<br/>equilibrium / stability / green / utils / constants"]
    C --> B
    B --> D["scalars and time series,<br/>written back into the ODS"]
    E["your own NumPy arrays"] --> C

vaft/formula/__init__.py resolves its submodules lazily (PEP 562): importing one of them costs only that module, and from vaft.formula import * still exposes every public symbol. Each is reachable three ways:

import vaft                                             # lazy: vaft.formula loads on first access
from vaft.formula import greenwald_density              # flat namespace
from vaft.formula.equilibrium import q_from_phi         # explicit submodule

Every public function documents its definition, units, convention, validity, limitations and references in one standardized layout; vaft.formula.describe("greenwald_density") prints it, and the Formula reference renders the same text per submodule.

Function names encode their inputs — confinement_time_from_engineering_parameters, beta_N_from_beta_a_B0_Ip, rho_star_from_M_T_B_R_epsilon. Where two conventions exist for the same quantity, they are two functions, not a flag.

Read the units before you call anything. This module deliberately mixes SI functions with published engineering-unit regressions, because those regressions are only valid in the units their authors used. Every unit trap on this page is collected again in Unit traps at a glance at the end.

Constants

from vaft.formula.constants import MU0, EPS0, QE, ME, MI_P, E_ALPHA

The same module also holds _SCALING_COEFS, the private confinement-scaling coefficient table read by confinement_time_from_engineering_parameters — see Confinement time and scalings below.

Flux, safety factor, shear

import numpy as np
import vaft

psi  = vaft.formula.psi_from_RBtheta(R, B_theta, l, psi_axis=0.0)   # psi = int R B_theta dl + psi_a
psiN = vaft.formula.psi_normalised(psi, psi_axis, psi_boundary)     # psi_N = (psi-psi_a)/(psi_b-psi_a)

phi  = vaft.formula.phi_from_Bphi(B_phi, dA)                        # Phi = int B_phi dA
rhoN = vaft.formula.rhoN_from_phi(phi, phi_boundary)                # rho_N = sqrt(Phi/Phi_b)

q    = vaft.formula.q_from_phi(psi, phi)                            # q = dPhi/dpsi
rhoN = vaft.formula.rhoN_from_qpsiN(psiN, q)                        # rho_N from the q integral
s    = vaft.formula.shear_from_r_q(r, q)                            # s = (r/q) dq/dr

q_from_rhoN(psiN, rhoN, C=1.0) goes the other way, with C a user-supplied prefactor. magnetic_shear is a live alias of shear_from_r_q; normalize_psi is a deprecated alias of psi_normalised and emits a DeprecationWarning.

Fields and current density from the flux map:

B_r = vaft.formula.radial_magnetic_field_from_psi(psi, R, Z)      # B_r = -(1/R) dpsi/dZ
B_z = vaft.formula.vertical_magnetic_field_from_psi(psi, R, Z)    # B_z = +(1/R) dpsi/dR
# current_density_from_psi (removed in 0.8.0, #355) returned -B_Z/mu0 [A/m],
# not a current density.  For that quantity, note the minus sign:
j   = -vaft.formula.vertical_magnetic_field_from_psi(psi, R, Z) / MU0   # [A/m]
# For a real toroidal current density use the Grad-Shafranov operator:
#   vaft.process.equilibrium.grad_shafranov_operator(psi, R, Z), or the
#   ODS's own profiles_2d.j_tor.

All of these differentiate with np.gradient along a single axis, so they expect 1-D slices, not a 2-D $(R, Z)$ map.

From an ODS

import vaft

ods = vaft.omas.sample_ods()                       # packaged VEST shot
p1d = ods['equilibrium.time_slice.0.profiles_1d']
s   = vaft.formula.shear_from_r_q(p1d['rho_tor_norm'], p1d['q'])   # s = (rho/q) dq/drho

shear_from_r_q accepts any monotonic radius label, so the rho_tor_norm every equilibrium carries is enough. vaft.omas.compute_magnetic_shear(ods, time_slice) is the same call with the abscissa fixed to equilibrium.time_slice[i].profiles_1d.r; that node is not part of the IMAS profiles_1d and the packaged sample does not carry it, so the wrapper raises a KeyError naming the missing path unless you have stored a minor-radius array there yourself.

Geometry and derived scalars

V     = vaft.formula.volume_from_RZ_boundary(R_bdry, Z_bdry)              # 2 pi A_poly R_bar
kappa = vaft.formula.elongation_from_RZ_boundary(R_bdry, Z_bdry)          # kappa = (Zmax-Zmin)/(2a)
delta = vaft.formula.triangularity_from_RZ_boundary(R_bdry, Z_bdry, R0)   # mean of the two below
d_u   = vaft.formula.triangularity_upper_from_RZ_boundary(R_bdry, Z_bdry, R0)  # (R0 - R at Zmax)/a
d_l   = vaft.formula.triangularity_lower_from_RZ_boundary(R_bdry, Z_bdry, R0)  # (R0 - R at Zmin)/a
eK    = vaft.formula.eK_from_K(kappa)                                     # eK = (k^2-1)/(k^2+1)
eps   = vaft.formula.inverse_aspect_ratio_from_a_R(a, R)                  # epsilon = a/R
A     = vaft.formula.aspect_ratio_from_a_R(a, R)                          # A = R/a
PF    = vaft.formula.peaking_factor(x_central, x_volume_avg)              # PF = X(0)/<X>

volume_from_RZ_boundary multiplies the shoelace polygon area by $2\pi \bar{R}$, with $\bar{R}$ the arithmetic mean of the boundary points — an approximation, not the exact Pappus centroid. calc_inverse_aspect_ratio(a_m, R_geo_m) is the validated variant (it rejects non-positive input) used by the dimensionless-parameter chain.

Energy and resistivity:

W   = vaft.formula.thermal_energy_from_p_V(p, V)               # W_th = 3/2 <p> V (IMAS energy_thermal)
PV  = vaft.formula.stored_energy_from_p_V(p, V)                # int p dV = <p> V: NOT W_th, despite the name
PV  = vaft.formula.stored_energy_from_beta_V(beta, B0, V)      # <p> V again, from beta = 2 mu0 <p> / B0^2
eta = vaft.formula.spitzer_resistivity_from_T_e_Z_eff_ln_Lambda(T_e, Z_eff=2.0, ln_Lambda=17.0)

For a self-consistent $\ln\Lambda$ instead of the 17.0 default:

ln_lambda = vaft.formula.coulomb_logarithm_from_n_T(n_m3=4e19, T_eV=200.0)   # 30.9 - ln(sqrt(n)/T)

Beta, current and power limits

beta_N = vaft.formula.beta_N_from_beta_a_B0_Ip(beta_percent, a, B0, I_p_MA)   # beta [%], I_p [MA]
beta_p = vaft.formula.beta_pol_from_beta_tor(beta_tor, q_95, epsilon)   # beta_p = beta_t (q95/eps)^2
beta_t = vaft.formula.beta_tor_from_beta_pol(beta_pol, q_95, epsilon)   # beta_t = beta_p (eps/q95)^2

I_max  = vaft.formula.current_limit_from_q(q_95, a, B0)           # I_p = 2 pi a^2 B0 / (mu0 q95)
P_max  = vaft.formula.power_limit_from_beta(beta_N, B0, V)
P_max  = vaft.formula.power_limit_from_q(q_95, I_p, R0)

beta_N_from_beta_a_B0_Ip works in the Troyon convention, %·m·T/MA: pass $\beta$ in percent and $I_p$ in megaamperes (the parameters are named beta_percent and I_p_MA for that reason). It does not rescale, so a fraction and amperes give $10^{-8}$ times the number you would compare with 2.8.

The two $\beta_p \leftrightarrow \beta_t$ conversions are the cylindrical relation $\beta_p = \beta_t\,(q_{95}/\varepsilon)^2$ with $\varepsilon = a/R_0$ the inverse aspect ratio, which is a required, strictly positive argument (ValueError otherwise).

Kink safety factor, with three geometry models:

q_kink, q_min, beta_max, beta_crit, ip_max = vaft.formula.kink_safety_factor(
    R, a, kappa, Ip, Bt, 'ST')       # 'circular' | 'conventional' | 'ST'

Only the 'conventional' and 'ST' branches return beta_max / beta_crit; 'circular' returns None for both. Any other string raises ValueError. The ST branch uses

Normalised plasma current, the Troyon normalisation $I_N = I_p/(a B_t)$ (Plasma Phys. Control. Fusion 26, 209) and the current axis of the Menard et al. low-aspect-ratio stability maps (Nucl. Fusion 37, 595):

Ip_star = vaft.formula.normalized_plasma_current(Ip, R, a, Bt)   # I_p [A] in, MA/(m T) out

Stability limits

from vaft.formula import (
    greenwald_density, greenwald_fraction,
    ballooning_alpha_from_p_B_R, ballooning_stability_criterion,
    kink_stability_criterion, sawtooth_stability_criterion,
    beta_stability_boundary, plasma_stability_margins,
)

n_G = greenwald_density(I_p=0.1, a=0.4)        # I_p in MA  ->  n_G in 1e19 m^-3
f_G = greenwald_fraction(n_e=1.5, n_G=n_G)     # same units and definition on both arguments

The implementation returns $n_G$ in units of $10^{19}\ \mathrm{m^{-3}}$ (that is, $10 I_p / \pi a^2$), and the Greenwald fraction is conventionally evaluated with the line-averaged density — pair them accordingly.

The local and global criteria all return a (margin, critical_value) pair, so the sign of the first element is the answer:

alpha        = ballooning_alpha_from_p_B_R(p, B, R, q)       # alpha = -2 mu0 R q^2 (dp/dR) / B^2
d_alpha, a_c = ballooning_stability_criterion(alpha, s)      # alpha_crit = 0.6 s
d_beta, b_c  = kink_stability_criterion(q_95, beta_N)        # deprecated (#350): beta_N_crit = 2.8 q95
d_bp,   bp_c = sawtooth_stability_criterion(q_0, beta_pol)   # deprecated (#350): beta_p_crit = 0.3 (1 - q0)
d_bN,   bN_c = beta_stability_boundary(beta_N, q_95)         # deprecated (#350): beta_N_crit = 0.028 q95

beta_margin, q_margin, density_margin = plasma_stability_margins(beta_N, q_95, n_e, n_G)  # deprecated (#350)

The four heuristics above are deprecated (#350): no source multiplies the $\beta_N$ limit by $q_{95}$, and Porcelli’s sawtooth trigger is a set of conditions, not $0.3\,(1 - q_0)$. The sourced limits live in vaft.formula.boundaries ("troyon", "low_q", "greenwald", "murakami"), evaluated with evaluate_boundary.

kink_stability_criterion and beta_stability_boundary use different critical-$\beta_N$ coefficients (2.8 versus 0.028) — the same Troyon-type relation written in two unit conventions (%·m·T/MA versus the dimensionless fraction). plasma_stability_margins is built on beta_stability_boundary, so it lives in the 0.028 convention; its q_margin is simply $q_{95} - 2$.

Legacy vertices of the JET empirical $l_i$–$q_\psi$ stability boundary (Wesson et al., Nucl. Fusion 29, 1989, Fig. 6, lower boundary). The audited boundaries, with provenance and quantity identity, are the registered wesson_1989_jet_li_qpsi_* entries of vaft.formula.boundaries (see Stability, “Internal inductance against edge q”):

qa_ref, li_ref = vaft.formula.empirical_li_qa()          # 18 surveyed points
li             = vaft.formula.li_from_qa_empirical(qa)   # piecewise-linear interpolation

Characteristic speeds and collisionality:

v_A     = vaft.formula.v_alfven_from_B_n_mi(B, n)            # n in m^-3; m_i defaults to the proton mass
c_s     = vaft.formula.c_s_from_Te_Ti_mi(T_e_keV, T_i_keV)   # keV in, m/s out
nu_star = vaft.formula.collisionality_from_n_T_B_R(n_e, T_e_keV, B_t, R0)

collisionality_from_n_T_B_R takes $n_e$ in $10^{19}\ \mathrm{m^{-3}}$, $T_e$ in keV, $B_t$ in T and $R_0$ in m — note that v_alfven_from_B_n_mi on the line above wants $n$ in $\mathrm{m^{-3}}$, so the two density arguments are not interchangeable.

Confinement time and scalings

The measured confinement time is a one-liner; the fitted one goes through the scaling evaluator:

tau_exp = vaft.formula.confinement_time_from_P_loss_W_th(P_loss, W_th)   # tau = W_th / P_loss

tau_iter89p = vaft.formula.confinement_time_from_engineering_parameters(
    I_p=1.0e5,        # [A]    -> MA internally
    B_t=0.18,         # [T]
    P_loss=3.0e5,     # [W]    -> MW internally
    n_e=2.0e19,       # [m^-3] -> 1e19 m^-3 internally
    M=1.0,            # [amu]
    R=0.4,            # [m]
    epsilon=0.75,     # a/R [-]
    kappa=1.6,        # [-]
    scaling="ITER89P",             # default
    input_density_definition="line_avg",
    line_to_volume_factor=None,    # required only when definitions differ
)

H = vaft.formula.confinement_factor_ITER89P(tau_exp, tau_iter89p)   # H = tau_exp / tau_ITER89P

with the product running only over the variables the selected scaling actually declares. Eight scaling names are accepted; anything else raises ValueError listing the available keys.

The ITPA20 fits add a triangularity factor $(1+\delta)^{\alpha}$: pass delta= (the average of the upper and lower LCFS triangularity) or they raise ValueError; the other scalings ignore it. Their kappa is the volume elongation $\kappa_a = V/(2\pi R\,\pi a^2)$ and their power the thermal loss power through the LCFS with radiation not subtracted. The paper’s ITER point (15 MA, 5.3 T, $10.3\times10^{19}$ m⁻³, 87 MW, 6.2 m, $\delta = 0.48$, $\kappa_a = 1.7$, $\epsilon = 0.32$, $M = 2.5$) gives 3.08 s (ITPA20) and 2.89 s (ITPA20-IL) with the printed coefficients, against the paper’s 3.07 s and 2.90 s from the unrounded fits. ITPA20-IL was fitted on the paper’s ITER-like subset (Sec. 2.1: eight conventional tokamaks, $q_{95} > 2.8$, $1.3 < \kappa < 2.2$, $\epsilon < 0.5$, $Z_{eff} < 5$), so it is an ITER extrapolation tool, not a fit for a spherical tokamak.

On the ODS path, vaft.omas.formula_wrapper.compute_tau_E_scaling takes $\delta$ from equilibrium.time_slice[i].boundary: the mean of triangularity_upper and triangularity_lower when both exist, else triangularity. Without either, an ITPA20 scaling does not raise there: the slice gets NaN and a warning that names the missing triangularity, as the wrapper does for every input it cannot map. That wrapper passes boundary.elongation (b/a) as kappa for every scaling, not $\kappa_a$.

The two NSTX entries omit R, epsilon, kappa and Mi on purpose: only the dependencies the paper publishes are implemented, and unspecified ones are not assumed. The signature is the same for every scaling, so you still have to supply those arguments — they are simply not raised to any power, and not even range-checked, when the selected scaling ignores them.

Things the evaluator does for you, and things it does not:

• It converts SI to engineering units. I_p A → MA, P_loss W → MW, n_e m⁻³ → $10^{19}$ m⁻³. B_t, R, epsilon, kappa, M pass through unchanged. Do not pre-scale.
• The prefactor $C$ is unit-dependent, not universal. The NSTX coefficients were converted from the papers’ SI-like form into this common convention, so the stored $C$ is only correct with the inputs above.
• Every variable used by the selected scaling must be finite and strictly positive, or you get a ValueError naming the offending argument.
• Density definitions are explicit. Each table entry declares a density_definition; all eight currently declare "line_avg". If your input_density_definition differs from the entry’s target, you must also pass line_to_volume_factor — the evaluator refuses to guess and raises ValueError.

From an ODS

These wrappers need more than the packaged sample holds — equilibrium...boundary.elongation, a core_profiles density profile and the power balance — so run them on a database shot whose equilibrium and profile stages are filled:

ods  = vaft.database.load(39915, source="public")

eng  = vaft.omas.compute_tau_E_engineering_parameters(ods, 0)   # dict: I_p, B_t, P_loss,
                                                                # n_e_line_avg, n_e_vol_avg,
                                                                # R, epsilon, kappa, M
tau_fit = vaft.omas.compute_tau_E_scaling(ods, 0, scaling="H98y2")
tau_exp = vaft.omas.compute_tau_E_exp(ods, 0)

(tau_ITER89P, tau_H98y2, tau_NSTX2006H, tau_NSTX2006L,
 tau_Kurskiev2022, H_factor, tau_exp) = vaft.omas.compute_confiment_time_paramters(ods, 0)

compute_tau_E_scaling picks the line- or volume-averaged density out of eng to match whatever the requested scaling declares, so you never pass input_density_definition at the ODS level.

scaling= is required. The wrappers are skip-tolerant: an unknown key or a missing input logs a warning and returns NaN for that time slice instead of raising.

Two experimental confinement times, kept apart (#1905). compute_tau_E_exp and the P_loss above are the legacy modelled budget: Spitzer $\int\eta J^2 dV$ heating ($Z_\mathrm{eff}=2$, $\ln\Lambda=17$), modelled radiation (line + bremsstrahlung + synchrotron) always subtracted, and $W = \tfrac32\langle p\rangle V$ with $p = 2n_eT_e$. They keep that definition, and calling compute_tau_E_exp says so with a FutureWarning. The experimental Ohmic chain is vaft.omas.confinement.ohmic_confinement_ods:

• $P_\mathrm{OH} = I_p V_\mathrm{res}$, from the measured inboard-midplane loop voltage with the internal-inductive voltage removed (Romero 2010 eq. 24). Both terms, $L_i\dot I_p$ and $\tfrac14\mu_0R_0I_p\dot l_{i3}$, are kept.
• $P_\mathrm{net} = P_\mathrm{OH} - dW/dt$, with unmeasured radiation not subtracted.
• An explicit stored-energy source (w_source) and recorded windows.
• A NaN with a reason where a rate is missing, never $dW/dt = 0$.

$V_\mathrm{res}$ is inferred, not measured, and exact only at the plasma boundary. The core_profiles summary sheet carries both definitions in separate columns (CORE_PROFILES_DEFINITIONS).

Unit traps at a glance

Everything on this page that will silently give you a wrong number if you feed it the wrong convention: