Equilibrium Module

DirectBField

Internal mutable struct to hold B-field components and their derivatives at a point. It is used as a temporary workspace to avoid allocations in tight loops.

DirectRunInput(...)

A container struct that bundles all necessary inputs for the direct_run function. It is created by one of the equilibrium file read-in functions after processing the raw equilibrium data and preparing the initial splines.

Fields

• config::EquilibriumConfig The equilibrium configuration object.

• sq_in 1D spline data versus normalized poloidal flux psin. Quantities:

  1. F = R * B_t — toroidal flux function [m·T]
  2. μ₀ * Pressure — plasma pressure (non-negative) [T²]
  3. q — safety factor profile
  4. √ψ_norm — square root of normalized flux

• psi_in 2D cubic interpolant on the (R, Z) grid [m]. The values correspond to the poloidal flux adjusted to be zero at the boundary [Wb/rad]. Definitions:

  1. ψ(R, Z) = ψ_boundary - ψ(R, Z)

  2. ψ = ψ * sign(ψ(centerR, centerZ))

     • 1D profiles are represented by CubicInterpolant or CubicSeriesInterpolant
     • 2D flux surfaces by CubicInterpolantND

• psi_in_xs::Vector{Float64} — R coordinate grid for psi_in [m]

• psi_in_ys::Vector{Float64} — Z coordinate grid for psi_in [m]

• rmin::Float64 — Minimum R-coordinate of the computational grid [m]

• rmax::Float64 — Maximum R-coordinate of the computational grid [m]

• zmin::Float64 — Minimum Z-coordinate of the computational grid [m]

• zmax::Float64 — Maximum Z-coordinate of the computational grid [m]

• psio::Float64 — Total flux difference |ψ_axis - ψ_boundary| [Wb/rad]

Outer constructor for EquilibriumConfig from a parsed TOML dictionary

EquilibriumConfig(...)

A mutable struct containing configuration parameters for equilibrium reconstruction. Bundles all necessary settings originally specified in the equil fortran namelists.

Fields

• eq_type::String - Type of equilibrium file ("efit", "solovev", "lar", etc.)
• eq_filename::String - Path to equilibrium input file
• jac_type::String - Jacobian coordinate type ("hamada", "pest", "equal_arc", "boozer", "park", "other")
• power_bp::Int - Poloidal field power exponent for Jacobian
• power_b::Int - Toroidal field power exponent for Jacobian
• power_r::Int - Major radius power exponent for Jacobian
• r0exp::Float64 - Major radius normalization for CHEASE/EQDSK [m]
• b0exp::Float64 - On-axis toroidal field normalization for CHEASE/EQDSK [T]
• grid_type::String - Grid type for flux surface discretization ("ldp", etc.)
• psilow::Float64 - Lower limit of normalized flux coordinate
• psihigh::Float64 - Upper limit of normalized flux coordinate
• mpsi::Int - Number of radial grid points
• mtheta::Int - Number of poloidal grid points
• newq0::Int - Override for on-axis safety factor (0 = use input value)
• etol::Float64 - Error tolerance for equilibrium solver
• force_termination::Bool - Terminate after equilibrium setup (skip stability calculations)
• use_galgrid::Bool - Use the same grid as galerkin method

Outer constructor for EquilibriumConfig that enables a toml file interface for specifying the configuration settings

DEPRECATED: Use [Equilibrium] section in gpec.toml instead

EquilibriumParameters

A mutable struct containing computed equilibrium parameters and diagnostic flags.

Fields

• ro::Union{Nothing,Float64} - R-coordinate of the magnetic axis [m]
• zo::Union{Nothing,Float64} - Z-coordinate of the magnetic axis [m]
• psio::Union{Nothing,Float64} - Total flux difference |ψaxis - ψboundary| [Wb/rad]
• rsep::Union{Nothing,Vector{Float64}} - R-coordinates of the plasma boundary [m]
• zsep::Union{Nothing,Vector{Float64}} - Z-coordinates of the plasma boundary [m]
• rext::Union{Nothing,Vector{Float64}} - R-coordinates of the plasma edge [m]
• zext::Union{Nothing,Vector{Float64}} - Z-coordinates of the plasma edge [m]
• psi0::Union{Nothing,Float64} - Normalized poloidal flux at reference location
• b0::Union{Nothing,Float64} - Total magnetic field strength at the axis [T]
• q0::Union{Nothing,Float64} - Safety factor at the axis
• qmin::Union{Nothing,Float64} - Minimum safety factor in the plasma
• qmax::Union{Nothing,Float64} - Maximum safety factor in the plasma
• qa::Union{Nothing,Float64} - Safety factor at the plasma edge
• q95::Union{Nothing,Float64} - Safety factor at 95% flux surface
• qextrema_psi::Union{Nothing,Vector{Float64}} - Normalized flux values at q extrema
• qextrema_q::Union{Nothing,Vector{Float64}} - Safety factor values at extrema
• mextrema::Union{Nothing,Int} - Number of extrema in q-profile
• psi_norm::Union{Nothing,Float64} - Normalized poloidal flux
• b_norm::Union{Nothing,Float64} - Normalized magnetic field strength
• psi_axis::Union{Nothing,Float64} - Poloidal flux at the axis
• psi_boundary::Union{Nothing,Float64} - Poloidal flux at the boundary
• psi_boundary_norm::Union{Nothing,Float64} - Normalized boundary flux
• psi_axis_norm::Union{Nothing,Float64} - Normalized axis flux
• psi_boundary_offset::Union{Nothing,Float64} - Boundary flux offset
• psi_axis_offset::Union{Nothing,Float64} - Axis flux offset
• psi_boundary_sign::Union{Nothing,Int} - Sign of boundary flux
• psi_axis_sign::Union{Nothing,Int} - Sign of axis flux
• psi_boundary_zero::Union{Nothing,Bool} - Whether boundary flux is zero
• rmean::Union{Nothing,Float64} - Mean major radius [m]
• amean::Union{Nothing,Float64} - Mean minor radius [m]
• aratio::Union{Nothing,Float64} - Aspect ratio (R0/a)
• kappa::Union{Nothing,Float64} - Plasma elongation
• delta1::Union{Nothing,Float64} - Upper triangularity
• delta2::Union{Nothing,Float64} - Lower triangularity
• bt0::Union{Nothing,Float64} - Toroidal field at axis [T]
• crnt::Union{Nothing,Float64} - Plasma current [A]
• bwall::Union{Nothing,Float64} - Toroidal field at wall [T]
• verbose::Bool - Enable verbose output
• diagnose_src::Bool - Enable source data diagnostics
• diagnose_maxima::Bool - Enable extrema diagnostics
• volume::Union{Nothing,Float64} - Plasma volume [m³]
• betat::Union{Nothing,Float64} - Toroidal beta
• betan::Union{Nothing,Float64} - Normalized beta
• betaj::Union{Nothing,Float64} - Total beta
• betap1::Union{Nothing,Float64} - Poloidal beta (definition 1)
• betap2::Union{Nothing,Float64} - Poloidal beta (definition 2)
• betap3::Union{Nothing,Float64} - Poloidal beta (definition 3)
• li1::Union{Nothing,Float64} - Internal inductance (definition 1)
• li2::Union{Nothing,Float64} - Internal inductance (definition 2)
• li3::Union{Nothing,Float64} - Internal inductance (definition 3)

FieldLineDerivParams

A struct to hold constant parameters for the ODE integration, making them easily accessible within the derivative function direct_fieldline_der!.

InverseRunInput(...)

A container struct for inputs to the inverse_run function.

Fields

• config::EquilibriumConfig - The equilibrium configuration object
• sq_in::CubicSeriesInterpolant - 1D profile spline for F, P, q
• rz_in_xs::Vector{Float64} - ψ coordinate grid for rz_in
• rz_in_ys::Vector{Float64} - θ coordinate grid for rz_in
• rz_in_R::CubicInterpolantND - R coordinate interpolant [m]
• rz_in_Z::CubicInterpolantND - Z coordinate interpolant [m]
• ro::Float64 - R-coordinate of magnetic axis [m]
• zo::Float64 - Z-coordinate of magnetic axis [m]
• psio::Float64 - Total flux difference |ψaxis - ψboundary| [Wb/rad]

LargeAspectRatioConfig(...)

A mutable struct holding parameters for the Large Aspect Ratio (LAR) plasma equilibrium model.

Fields:

• lar_r0: The major radius of the plasma [m].
• lar_a: The minor radius of the plasma [m].
• beta0: The beta value on axis (normalized pressure).
• q0: The safety factor on axis.
• p_pres: The exponent for the pressure profile, defined as p00 * (1 - (r / a)^2)^p_pres.
• p_sig: The exponent that determines the shape of the current-related function profile.
• sigma_type: The type of sigma profile, can be "default" or "wesson". If "wesson", the sigma profile is defined as sigma0 * (1 - (r / a)^2)^p_sig.
• mtau: The number of grid points in the poloidal direction.
• ma: The number of grid points in the radial direction.
• zeroth: If set to true, it neglects the Shafranov shift

Outer constructor for LargeAspectRatioConfig that enables a toml file interface for specifying the configuration settings

PlasmaEquilibrium(...)

The final, self-contained result of the equilibrium reconstruction. This object provides a complete representation of the processed plasma equilibrium in flux coordinates.

Fields

• config::EquilibriumConfig: The equilibrium configuration object used for the reconstruction.

• params::EquilibriumParameters: Computed equilibrium parameters and diagnostics.

• profiles::ProfileSplines: Named 1D profile splines (F, P, dV/dψ, q) on normalized psi grid. Access values at grid points via profiles.F_spline.y[i], etc. Access derivatives via profiles.F_deriv.y[i] or profiles.F_deriv(psi).

• Grid coordinates (shared by all rzphi/eqfun interpolants):

  • rzphi_xs::Vector{Float64}: ψ coordinates (length mpsi+1)
  • rzphi_ys::Vector{Float64}: θ coordinates (length mtheta+1)

• Geometric quantities (rzphi, 4 interpolants): 2D cubic interpolants for flux-coordinate mapping with periodic BC in theta.

  • x value: normalized ψ
  • y value: SFL poloidal angle ∈ [0, 1]
  • rzphi_rsquared::CubicInterpolantND: r_coord² = (R - ro)² + (Z - zo)²
  • rzphi_offset::CubicInterpolantND: η/(2π) - θₙₑw (angle offset)
  • rzphi_nu::CubicInterpolantND: ν in ϕ = 2πζ + ν(ψ, θ)
  • rzphi_jac::CubicInterpolantND: Jacobian

• Physics quantities (eqfun, 3 interpolants): 2D cubic interpolants storing local physics and geometric quantities.

  • x value: normalized ψ
  • y value: SFL poloidal angle θₙₑw
  • eqfun_B::CubicInterpolantND: Total magnetic field strength [T]
  • eqfun_metric1::CubicInterpolantND: (e₁⋅e₂ + q⋅e₃⋅e₁)/(J⋅B²)
  • eqfun_metric2::CubicInterpolantND: (e₂⋅e₃ + q⋅e₃⋅e₃)/(J⋅B²)

• ro::Float64: R-coordinate of the magnetic axis [m]

• zo::Float64: Z-coordinate of the magnetic axis [m]

• psio::Float64: Total flux difference |Ψaxis - Ψboundary| [Weber/radian]

ProfileSplines

Named 1D cubic spline interpolants for equilibrium profiles. Each profile is stored as a separate spline for code clarity.

Fields

• xs::Vector{Float64}: Shared x-axis (normalized psi)
• F_spline: 2π*F (toroidal flux function, where F = R * B_toroidal)
• P_spline: μ₀*P (plasma pressure × μ₀)
• dVdpsi_spline: dV/dψ (volume derivative)
• q_spline: q (safety factor)

Derivative Interpolants (for continuous derivative evaluation)

• F_deriv, P_deriv, dVdpsi_deriv, q_deriv: First derivative interpolants

Notes

• Node values at grid points: Access via spline.y[i]
• Derivative at any point: Call deriv(x) (derivative views are callable)
• Grid: Access via xs field or spline.cache.x

ProfileSplines(xs, F_vals, P_vals, dVdpsi_vals, q_vals; extrap=ExtendExtrap())

Create ProfileSplines from arrays of profile values. Uses CubicFit boundary conditions with extension extrapolation.

SolovevConfig(...)

A mutable struct holding parameters for the Solev'ev (SOL) plasma equilibrium model.

Fields:

• mr: number of radial grid zones
• mz: number of axial grid zones
• ma: number of flux grid zones
• e: elongation
• a: minor radius
• r0: major radius
• q0: safety factor at the o-point
• p0fac: scale on-axis pressure (P-> P+P0*p0fac. beta changes. Phi,q constant)
• b0fac: scale toroidal field at constant beta (sPhi,sf,s^2*P. bt changes. Shape,beta constant)
• f0fac: scale toroidal field at constant pressure (s*f. beta,q changes. Phi,p,bp constant)

Outer constructor for SolovevConfig that enables a toml file interface for specifying the configuration settings

read1dgfileformat(linesblock, numvalues)

Internal helper function to parse Fortran-style fixed-width numerical blocks from a vector of strings.

Arguments:

• lines_block: A Vector{String} containing the lines to parse.
• num_values: The total number of Float64 values to read from the block.

Returns:

• A Vector{Float64} containing the parsed values.

direct_fieldline_der!(dy, y, params, eta)

The derivative function for the field-line integration ODE. This is passed to the DifferentialEquations.jl solver. This is a Julia adaptation of the Fortran direct_fl_der subroutine, with an added safeguard against division by zero.

Arguments:

• dy: The derivative vector (output, modified in-place).
• y: The state vector [∫(dl/Bp), rfac, ∫(dl/(R²Bp)), ∫(jac*dl/Bp)].
• params: A FieldLineDerivParams struct with all necessary parameters.
• eta: The independent variable (geometric angle η).

direct_fieldline_int(psifac, raw_profile, ro, zo, rs2)

Performs the field-line integration for a single flux surface. This is a Julia adaptation of the Fortran direct_fl_int subroutine. Note that the array y_out is now indexed from 1:5 rather than 0:4 as in Fortran.

Arguments:

• psifac: normalized psi value for the surface (ψ_norm).
• raw_profile: DirectRunInput object containing splines and parameters.
• ro, zo: Coordinates of the magnetic axis [m].
• rs2: R-coordinate of the outboard separatrix [m].

Returns:

- `y_out`: A matrix containing the integrated quantities vs. the geometric angle `η`.
    - `y_out[:, 1]`: η (geometric poloidal angle)
    - `y_out[:, 2]`: ∫(dl/Bp)
    - `y_out[:, 3]`: rfac (radial distance from magnetic axis)
    - `y_out[:, 4]`: ∫(dl/(R²Bp))
    - `y_out[:, 5]`: ∫(jac*dl/Bp)

• bfield: A DirectBField object with values at the integration start point.

direct_get_bfield!(bf_out, r, z, psi_in, sq_in, sq_in_deriv, psio; derivs=0)

Calculates the magnetic field and its derivatives at a given (R,Z) point. The results are stored in-place in the bf_out object. This is equivalent to the direct_get_bfield subroutine in the Fortran code, adapted for the Julia spline implementation.

Arguments:

• bf_out: A mutable DirectBField struct to store the results
• r: R-coordinate to evaluate at
• z: Z-coordinate to evaluate at
• psi_in: 2D cubic interpolant for poloidal flux ψ(R,Z)
• sq_in: 1D cubic spline for profiles F(ψ_norm) and P(ψ_norm)
• sq_in_deriv: Pre-computed derivative view of sq_in
• psio: total toroidal flux
• derivs: An integer specifying number of derivatives to compute (0, 1, or 2)

direct_position!(raw_profile)

Finds the key geometric locations of the equilibrium: the magnetic axis (O-point) and the inboard/outboard separatrix crossings on the midplane. It also updates the spline representing the poloidal flux ψ(R,Z) based on the new magnetic axis location. This function performs the same overall function as the Fortran direct_position subroutine with better iteration control and error handling. We have also added a helper function for separatrix finding.

Arguments:

• raw_profile: A DirectRunInput object containing splines and parameters.

Returns:

• ro: R-coordinate of the magnetic axis [m].
• zo: Z-coordinate of the magnetic axis [m].
• rs1: R-coordinate of the inboard separatrix crossing [m].
• rs2: R-coordinate of the outboard separatrix crossing [m].
• psi_in_new : returns psi_in renormalized by * psio/psi(ro,zo)

direct_refine(rfac, eta, psi0, params)

Refines the radial distance rfac at a given angle eta to ensure the point lies exactly on the target flux surface psi0.

Arguments:

• rfac: The current guess for the radial distance from the magnetic axis.
• eta: The geometric poloidal angle.
• psi0: The target ψ value for the flux surface.
• params: A FieldLineDerivParams struct.

Returns:

• The refined rfac value.

equilibrium_global_parameters!(pe::PlasmaEquilibrium)

Computes and populates global equilibrium parameters in the PlasmaEquilibrium struct, such as rmean, amean, kappa, bt0, crnt, betat, betan, li1, etc. Performs the same function as equiloutglobal in the Fortran code.

equilibrium_gse!(equil::PlasmaEquilibrium)

Diagnoses the Grad-Shafranov solution by computing the residual of the Grad-Shafranov equation across the grid and writing diagnostic data to HDF5 files. Performs the same function as equiloutgse in the Fortran code.

equilibrium_qfind!(equil::PlasmaEquilibrium)

Finds the extrema of the safety factor profile q(ψ) in the plasma equilibrium and computes derived q-values such as q0, qmin, qmax, qa, and q95. Performs the same function as equiloutqfind in the Fortran code.

equilibrium_separatrix_find!(pe::PlasmaEquilibrium)

Finds the separatrix locations in the plasma equilibrium (rsep, zsep, rext, zext). Performs the same function as equiloutsep_find in the Fortran code.

equilibrium_solver(raw_profile)

The main driver for the direct equilibrium reconstruction. It orchestrates the entire process from finding the magnetic axis to integrating along field lines and constructing the final coordinate and physics quantity splines. This performs the same overall function as the Fortran direct_run subroutine, with better checks for numerical robustness.

Arguments:

• raw_profile: A DirectRunInput object containing the initial splines (psi_in, sq_in) and run parameters (equil_input).

Returns:

• A PlasmaEquilibrium object containing the final, processed equilibrium data, including the profile spline (sq), the coordinate mapping spline (rzphi), and the physics quantity spline (eqfun).

inverse_extrap(xx::Matrix{Float64}, ff::Matrix{Float64}, x::Float64) -> Vector{Float64}

Performs component-wise Lagrange extrapolation for a vector-valued function.

Arguments:

• xx: A (m × n) matrix where each row contains the x-values for each component.
• ff: A (m × n) matrix where each row contains function values at the corresponding xx.
• x: A scalar Float64 value at which to extrapolate.

Returns:

• A vector of length n representing the extrapolated function values at x.

lar_init_conditions(rmin, sigma_type, params)

Initializes the starting radius and state vector for solving the LAR ODE system. Also evaluates the initial derivative using the analytic model.

Arguments:

• rmin: Normalized starting radius (as a fraction of lar_a).
• lar_input: A LargeAspectRatioConfig object containing equilibrium parameters.

Returns:

• r: Physical radius corresponding to rmin * lar_a.
• y: Initial state vector of length 5.

read_chease_ascii(config)

Parses a ascii CHEASE file, creates initial 1D and 2D splines, finds magnetic axis, and bundles them into a InverseRunInput object.

Arguments:

• config: The EquilibriumConfig object containing the filename and parameters.

Returns:

• A InverseRunInput object ready for the inverse solver.

read_chease_binary(equil_config)

Parses a binary CHEASE file, creates initial 1D and 2D splines with proper normalization (R0, B0 scaling), and bundles them into a InverseRunInput object.

_read_efit(equil_in)

Parses an EFIT g-file, creates initial 1D and 2D splines, and bundles them into a DirectRunInput object.

Arguments:

• equil_in: The EquilInput object containing the filename and parameters.

Returns:

• A DirectRunInput object ready for the direct solver.

setup_equilibrium(eq_config::EquilibriumConfig)

The main public API for the Equilibrium module. It orchestrates the entire process of reading an equilibrium file, running the appropriate solver, and returning the final processed PlasmaEquilibrium object.

Arguments:

• eq_config: An EquilibriumConfig object containing all necessary setup parameters.

Returns:

• A PlasmaEquilibrium object containing the final result.

This function handles the Solovev analytical equilibrium model, transforming the input parameters into the necessary splines and scalar values for equilibrium construction. This is a Julia version of the Fortran code in sol.f, with no major differences except for arrays going from 0:n to 1:n+1.

Arguments:

• mr: Number of radial grid zones
• mz: Number of axial grid zones
• ma: Number of flux grid zones
• e: Elongation
• a: Minor radius
• r0: Major radius
• q0: Safety factor at the o-point
• p0fac: Scales on axis pressure (s*P. beta changes. Phi,q constant)
• b0fac: Scales on toroidal field (sPhi,sf,s^2*P. bt changes. Shape,beta constant)
• f0fac: Scales on toroidal field (s*f. bt,q changes. Phi,p,bp constant)

Returns:

• DirectRunInput object