Vacuum Module

GaussLegendreRule{N,T}

Allocation-free Gauss–Legendre nodes/weights on the canonical interval [-1, 1]. Stored as SVectors so tight loops can index them efficiently.

KernelWorkspace

Thread-local workspace for compute_3D_kernel_matrices! to enable parallel execution. Each thread gets its own workspace to avoid data races on temporary arrays.

KernelWorkspace(PATCH_DIM, RAD_DIM, ANG_DIM)

Create a new workspace with pre-allocated arrays for kernel matrix computation.

PlasmaGeometry

Struct holding plasma geometry data on the mtheta grid for vacuum calculations. Arrays are of length mtheta, where mtheta is the number of poloidal grid points and θ ∈ [0, 1).

Fields

• x::Vector{Float64}: Plasma surface R-coordinate on VACUUM theta grid
• z::Vector{Float64}: Plasma surface Z-coordinate on VACUUM theta grid
• ν::Vector{Float64}: Magnetic toroidal angle offset from geometric toroidal angle

PlasmaGeometry(inputs::VacuumInput) -> PlasmaGeometry

Initialize the plasma surface geometry based on the provided vacuum inputs. We interpolate the input plasma boundary arrays from the inputs struct onto the mtheta grid.

Arguments

• inputs::VacuumInput: Struct containing plasma boundary data

Returns

• PlasmaGeometry: Struct containing plasma surface coordinates, derivatives, and basis functions

PlasmaGeometry3D

3D toroidal surface geometry for vacuum boundary integral calculations.

Built by toroidally extruding a 2D poloidal contour (PlasmaGeometry) and computing Cartesian coordinates, tangent vectors, normals, and differential area elements.

Fields

• mtheta::Int: Number of poloidal grid points
• nzeta::Int: Number of toroidal grid points
• r::Matrix{Float64}: Surface points in Cartesian (X,Y,Z), shape (num_points, 3)
• dr_dθ::Matrix{Float64}: Poloidal tangent vector ∂r/∂θ × dθ, shape (num_gridpoints, 3)
• dr_dζ::Matrix{Float64}: Toroidal tangent vector ∂r/∂ζ × dζ, shape (num_points, 3)
• normal::Matrix{Float64}: Oriented normal vectors, shape (num_points, 3)
• normal_orient::Int: Forces normals to face out from vacuum region (+1 or -1)

PlasmaGeometry3D(inputs::VacuumInput)

Construct a 3D toroidal plasma surface from vacuum input data.

This constructor builds a PlasmaGeometry3D directly from the VacuumInput struct, handling both axisymmetric (2D boundary, nzeta_in == 1) and fully 3D input boundaries (nzeta_in > 1).

Axisymmetric input (inputs.nzeta_in == 1)

1. Build a 2D poloidal contour on the vacuum mtheta grid using PlasmaGeometry(inputs) to obtain R(theta), Z(theta), and nu(theta).
2. Toroidally extrude this contour onto a uniform nzeta grid using the SFL angle zeta = phi - nu(theta) and map to Cartesian coordinates: X = R(theta) * cos(zeta - nu(theta)), Y = R(theta) * sin(zeta - nu(theta)), Z = Z(theta).

Fully 3D input (inputs.nzeta_in > 1)

1. Interpolate the input (x, y, z) arrays from the original mthetain × nzetain grid onto the vacuum mtheta × nzeta grid using periodic bicubic interpolation in both angles. The inputs are assumed to already be equally spaced on the SFL angle grid.

Steps

1. Fit periodic bicubic splines to each Cartesian component on the (theta, zeta) grid.
2. Compute tangent vectors dr/dtheta and dr/dzeta from spline derivatives, scaled by the grid spacings.
3. Form oriented normals via the cross product n = (dr/dtheta) × (dr/dzeta) and enforce a consistent orientation (inward for the plasma surface).
4. Compute average poloidal/toroidal grid spacings and report the aspect ratio for diagnostics.

Arguments

• inputs::VacuumInput: Vacuum calculation inputs defining the boundary geometry and the desired mtheta, nzeta resolution.

Returns

• PlasmaGeometry3D: Complete 3D surface description on the mtheta × nzeta grid, including points, tangents, normals, and orientation.

PnQuadEntry

Cached sinh/cosh values for the 32-point Gaussian quadrature at a given toroidal mode number n. Each field is a 32-element vector indexed by Gauss node ig.

SingularQuadratureData

Precomputed data for singular correction quadrature following BIEST approach. Initialized once on first use.

Fields

- `qx::Vector{Float64}`: Radial quadrature points in [0,1]
- `qw::Vector{Float64}`: Radial quadrature weights
- `Gpou::Matrix{Float64}`: Partition of unity on Cartesian grid (PATCH_DIM × PATCH_DIM)
- `Ppou::Matrix{Float64}`: Partition of unity on polar grid (RAD_DIM × ANG_DIM)
- `P2G::SparseMatrixCSC{Float64,Int}`: Sparse interpolation matrix (Ngrid × Npolar) mapping polar quadrature points to Cartesian grid
    - Forward (patch→polar): `polar = P2G' * patch`
    - Backward (polar→grid): `grid = P2G * polar`.
- `PATCH_DIM::Int`: Patch dimension (odd integer)
- `PATCH_RAD::Int`: Patch radius (number of points adjacent to source point treated as singular)
- `ANG_DIM::Int`: Number of angular quadrature points
- `RAD_DIM::Int`: Number of radial quadrature points
- `INTERP_ORDER::Int`: Lagrange interpolation order

SingularQuadratureData(PATCH_RAD::Int, RAD_DIM::Int, INTERP_ORDER::Int)

Constructor which initializes quadrature points, weights, partition-of-unity functions, and interpolation matrices for singular correction based on input parameters. Follows BIEST's approach.

Arguments

• PATCH_RAD::Int: Number of points adjacent to source point to treat as singular
• RAD_DIM::Int: Radial quadrature order
• INTERP_ORDER::Int: Lagrange interpolation order

Returns

• SingularQuadratureData: Precomputed quadrature data

VacuumInput

Struct holding plasma boundary and mode data as provided from ForceFreeStates namelist and computed quantities. For an axisymmetric boundary, nzetain = 1 and only the x and z arrays need to be provided - the code can then be run with nzeta = 1 for 2D vacuum calculation or nzeta > 1 for 3D vacuum calculation. For a non-axisymmetric boundary, nzetain > 1 and the x, y, and z arrays need to be provided - the code can then be run with nzeta = 1 for 2D vacuum calculation or nzeta > 1 for 3D vacuum calculation. The arrays should be for a single field period only, with excluded endpoints.

Fields

• x::Vector{Float64}: Plasma boundary X-coordinate (length mthetain * nzetain)
• y::Vector{Float64}: Plasma boundary Y-coordinate (length mthetain * nzetain)
• z::Vector{Float64}: Plasma boundary Z-coordinate (length mthetain * nzetain)
• ν::Vector{Float64}: Free parameter in specifying toroidal angle, ζ = ϕ + ν(θ), on input theta grid (axisymmetric only, length mtheta_in)
• mtheta_in::Int: Number of input poloidal grid points
• nzeta_in::Int: Number of input toroidal grid points (1 for axisymmetric, > 1 for non-axisymmetric)
• m_modes::Vector{Int}: Vector of poloidal mode numbers. E.g. collect(mlow:mhigh) for a contiguous range.
• n_modes::Vector{Int}: Vector of toroidal mode numbers. E.g. collect(nlow:nhigh) for a contiguous range or collect(nlow:n_stride:nhigh) for a strided list for stellarator mode-family calculations.
• mtheta::Int: Number of vacuum calculation poloidal grid points
• nzeta::Int: Number of vacuum calculation toroidal grid points (1 for 2D vacuum calculation, > 1 for 3D vacuum calculation)
• nfp::Int: Number of field periods

VacuumInput(
    equil::Equilibrium.PlasmaEquilibrium,
    ψ::Float64,
    mtheta::Int,
    nzeta::Int,
    m_modes::AbstractVector{<:Integer},
    n_modes::AbstractVector{<:Integer};
) -> VacuumInput

Constructor to create a VacuumInput struct for computing Green's functions at arbitrary flux surface. Extracts plasma geometry from equilibrium at the given flux surface and packages it into VacuumInput format.

Arguments

• equil: Equilibrium solution
• ψ: Normalized flux coordinate
• mtheta: Number of vacuum calculation poloidal points
• nzeta: Number of vacuum calculation toroidal points (1 for 2D, >1 for 3D)
• m_modes: Poloidal mode numbers (e.g. mlow:mhigh)
• n_modes: Toroidal mode numbers (e.g. [n] for a single mode, or nlow:nhigh)

Returns

VacuumInput structure ready for computevacuumresponse()

VacuumResponse

Output of compute_vacuum_response: the vacuum energy matrix and the surface data the boundary-integral solve produces along the way.

Fields

• wv::Matrix{ComplexF64}: Vacuum energy matrix Wᵛ (num_modes × num_modes), block-diagonal in n for 2D
• I_v::Matrix{ComplexF64}: Vacuum surface-current matrix Iᵛ (num_modes × num_modes), left zeroed unless compute_vacuum_response is called with compute_Iv=true (2D only). Stored without the μ₀/4π² normalization: the physical surface inductance is μ₀(2π)²·I_v⁻¹ (see PerturbedEquilibrium.calc_surface_inductance).
• plasma_pts, wall_pts::Matrix{Float64}: Cartesian surface coordinates (num_points × 3)

VacuumResponse(inputs::VacuumInput) -> VacuumResponse

Allocate zeroed output arrays sized for inputs (full torus, so nfp field periods).

WallGeometry

Struct holding wall geometry data for vacuum calculations. Arrays are of length mtheta, where mtheta is the number of poloidal grid points and θ ∈ [0, 1).

Fields

• nowall::Bool: Boolean flag indicating if there is no wall
• x::Vector{Float64}: Wall R-coordinates
• z::Vector{Float64}: Wall Z-coordinates

WallGeometry(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_settings::WallShapeSettings) -> WallGeometry

Constructor to initialize the wall geometry based on the provided vacuum inputs and wall shape settings.

This performs functionality similar to portions of the arrays function in the original Fortran VACUUM code. It returns a WallGeometry struct containing the necessary wall surface data for vacuum calculations.

Arguments

• inputs::VacuumInput: Struct containing vacuum calculation parameters
• plasma_surf::PlasmaGeometry: Struct with plasma surface geometry (used for reference)
• wall_settings::WallShapeSettings: Struct specifying wall shape and parameters

Returns

• WallGeometry: Struct containing wall surface coordinates and derivatives

Notes

• Supports multiple wall shapes: nowall, conformal, elliptical, dee, moddee, fromfile
• Optionally redistributes wall points to equal arc length spacing if equal_arc_wall=true

WallGeometry3D

Struct holding wall geometry data for vacuum calculations. Arrays are of length mtheta, where mtheta is the number of poloidal grid points and θ ∈ [0, 1).

Fields

• nowall::Bool: Boolean flag indicating if there is no wall
• is_closed_toroidal::Bool: Boolean flag indicating if the wall is a closed toroidal surface
• mtheta::Int: Number of poloidal grid points
• nzeta::Int: Number of toroidal grid points
• r::Matrix{Float64}: (x, y, z) wall coordinates at each grid point
• dr_dθ::Matrix{Float64}: Derivative dR/dθ at wall
• dr_dζ::Matrix{Float64}: Derivative dR/dζ at wall
• normal::Matrix{Float64}: Outward normal vectors at wall

WallGeometry3D(inputs::VacuumInput, wall_settings::WallShapeSettings)

Constructor to initialize the 3D wall geometry based on the provided vacuum inputs and wall shape settings. Currently only works for axisymmetric walls generated by toroidal extrusion of 2D poloidal contours.

Arguments

• inputs::VacuumInput: Struct containing vacuum calculation parameters
• wall_settings::WallShapeSettings: Struct specifying wall shape and parameters

Returns

• WallGeometry3D: Struct containing wall surface coordinates and derivatives

WallShapeSettings

Struct containing input settings for vacuum wall geometry.

Fields

• shape::String: String selecting wall shape. Options are:

  + `"nowall"`: No wall
  + `"conformal"`: Wall conformal to plasma surface at distance `a`
  + `"elliptical"`: Elliptical wall
  + `"dee"`: Dee-shaped wall
  + `"mod_dee"`: Modified Dee-shaped wall
  + `"filepath"`: Custom wall shape from the file you specify

• a::Float64: Distance of wall from plasma in units of major radius (conformal), or shape parameter (others)

• aw::Float64: Half-thickness parameter for Dee-shaped walls

• bw::Float64: Elongation parameter for wall shapes

• cw::Float64: Offset of the center of the wall from the major radius

• dw::Float64: Triangularity parameter for wall shapes

• tw::Float64: Sharpness of the corners of the wall (try 0.05 as initial value)

  Core shape selection

• equal_arc_wall::Bool: Flag to enforce equal arc length distribution of nodes on the wall (recommended unless wall is very close to plasma)

Pn_minus_half_1997(s, n)

Compute the Legendre function of the first kind of order -1/2, P^n_{-1/2}(s), recursively using Chance 1997 equations (47)-(50).

The implementation follows the original Fortran code. Note: equation (50) in the paper has a typo where the exponent should be -1/4 instead of +1/2.

Arguments

• s::Real: Legendre function parameter (s > 1)
• n::Int: Maximum order n (n ≥ 0)

Returns

• P::Vector{Float64}: Array of values P^0{-1/2}(s) through P^{n+1}{-1/2}(s)

Notes

• Uses recursive relation from Chance 1997 eq. (47)
• Base cases computed from eqs. (48)-(50) using elliptic integrals

Pn_minus_half_2007(s, n)

Compute the Legendre function of the first kind of order -1/2, P^n_{-1/2}(s), using methods from Chance J. Comp. Phys 221 (2007) 330-348.

This implementation uses:

1. Bulirsch's algorithm for elliptic integrals (more accurate than polynomial approximations)
2. Gaussian integration for large mode numbers (nrhohat >= 0.1) where rhohat = 1/√(2y*w)
3. Upward recurrence for small mode numbers

Arguments

• s::Real: Legendre function parameter (s > 1)
• n::Int: Maximum order n (n ≥ 0)

Returns

• P::Vector{Float64}: Array of values P^0{-1/2}(s) through P^{n+1}{-1/2}(s)

Notes

• This version is more accurate than Pnminushalf_1997 for large n
• Expected to diverge from 1997 version at large nloc
• Reference: JCP 221 (2007) 330-348 # Constants

_calculate_potential_chi(R_obs, Z_obs, inputs, plasma_surf, grri, Bn_real, Bn_imag)

Calculate the magnetic scalar potential chi at a single observation point (Robs, Zobs).

This is the Julia version of the Fortran chi subroutine. The potential is computed by integrating the Green's function response with the source perturbation at the plasma surface.

Arguments

• R_obs::Float64: R-coordinate of observation point
• Z_obs::Float64: Z-coordinate of observation point
• inputs::VacuumInput: Struct containing vacuum calculation parameters
• plasma_surf::PlasmaGeometry: Struct with plasma surface geometry
• grri::Matrix{Float64}: Inverted Green's function response matrix
• Bn_real::Vector{Float64}: Real part of normal field Fourier harmonics
• Bn_imag::Vector{Float64}: Imaginary part of normal field Fourier harmonics

Returns

• chi_real::Float64: Real part of the magnetic scalar potential at (Robs, Zobs)
• chi_imag::Float64: Imaginary part of the magnetic scalar potential at (Robs, Zobs)

Notes

• The potential is computed via Fourier series over poloidal modes
• Includes coupling term from Green's function derivative
• Factor of -0.5 * dtheta applied from Fortran convention

_compute_vacuum_response_2d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false)

2D (axisymmetric) vacuum response calculation.

Each toroidal mode n decouples in 2D geometry, so the routine loops over inputs.n_modes, building the double-/single-layer operators, solving the exterior system for wv, and optionally the interior system to build I_v when compute_Iv=true. Green's functions are internal scratch only.

_compute_vacuum_response_3d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false)

3D (inputs.nzeta > 1) vacuum response via block-circulant field-period reduction. For nfp == 1 the block-circulant assembly and residue-class loop are skipped in favour of a more efficient direct real solve on the built kernel matrices (equivalent to 2D).

The nfp-periodic boundary makes the single-/double-layer operators S, D block-circulant in the field-period index. Writing the response as wv = (4π²/N)·Eᴴ·D⁻¹S·E (E the complex Fourier basis, N = mtheta·nzeta_full), the structure block-diagonalizes the problem by toroidal residue class k = mod(n, nfp): modes with different k do not couple, and within a class

D̂ₖ = Σ_d D_d ω^{k d},   Ŝₖ = Σ_d S_d ω^{k d},   ω = exp(-2πi/nfp),

with D_d, S_d the blocks coupling observers in field period 0 to sources in period d. Each class needs one solve wv[class k] = (4π²/M)·E_localᴴ·(D̂ₖ \ Ŝₖ)|_plasma·E_local. Only the first block-row of the operators is built, so the kernel cost drops by nfp and the dense O(N³) factorization is replaced by per-class O(M³) solves (M = N/nfp).

Only wv is produced currently; I_v is left zeroed when compute_Iv=true (surface-current / inductance not yet supported in 3D). Extension point: per residue class, apply the per-period basis to D̂ₖ⁻¹Ŝₖ for the exterior columns and the interior variant -D + 2I for the interior columns, then scatter back into the [2N × 2·num_modes] arrays.

_create_pickup_grid(R_grid, Z_grid)

Create a flattened 1D list of (R, Z) coordinates from 2D grid vectors.

This is the Julia version of the Fortran loops subroutine.

Arguments

• R_grid::AbstractVector: Vector of R-coordinates defining the grid
• Z_grid::AbstractVector: Vector of Z-coordinates defining the grid

Returns

• R_points::Vector{Float64}: Flattened array of R-coordinates (length nx*nz)
• Z_points::Vector{Float64}: Flattened array of Z-coordinates (length nx*nz)

Notes

• Grid points are ordered as: [(R[1],Z[1]), (R[1],Z[2]), ..., (R[1],Z[nz]), (R[2],Z[1]), ...]

_pickup_field(inputs, plasma_surf, grri, Bn_real, Bn_imag, R_grid, Z_grid)

Calculate the magnetic field on a specified grid using finite differencing of the magnetic scalar potential chi.

This is the Julia version of the Fortran pickup routine. It computes the vacuum magnetic field perturbation at a set of grid points given the plasma surface perturbation.

Arguments

• inputs::VacuumInput: Struct containing vacuum calculation parameters
• plasma_surf::PlasmaGeometry: Struct with plasma surface geometry and basis functions
• grri::Matrix{Float64}: Inverted Green's function response matrix from vaccal!
• Bn_real::Vector{Float64}: Real part of normal field Fourier harmonics at plasma surface
• Bn_imag::Vector{Float64}: Imaginary part of normal field Fourier harmonics at plasma surface
• R_grid::AbstractVector: R-coordinates for output field evaluation
• Z_grid::AbstractVector: Z-coordinates for output field evaluation

Returns

• B_R::Matrix{ComplexF64}: R-component of magnetic field on grid (nx × nz)
• B_Z::Matrix{ComplexF64}: Z-component of magnetic field on grid (nx × nz)
• B_phi::Matrix{ComplexF64}: Toroidal component of magnetic field on grid (nx × nz)
• grid_info::Matrix{Int}: Grid point classification (1=inside plasma, 0=outside)

Notes

• Uses 5-point finite difference stencil for computing field from potential
• Field components computed as: BR = -∂χ/∂R, BZ = -∂χ/∂Z, B_φ = inχ/R

_warn_and_symmetrize!(mat, name)

Replace mat by its Hermitian part in place, warning first if the anti-Hermitian residual exceeds _HERMITICITY_WARN_TOL.

The matrices passed here are Hermitian in exact arithmetic — δW_v = ξ†Wᵛξ is a real energy, and Iᵛ is the inverse of the Hermitian surface inductance up to a real scalar — so any anti-Hermitian part is a discretization artifact that vanishes as the vacuum grid is refined.

compute_2D_kernel_matrices!(grad_greenfunction, greenfunction, observer, source, n)

Compute kernels of integral equation for Laplace's equation in a torus. WARNING: This kernel only supports closed toroidal walls currently. The residue calculation needs to be updated for open walls.

Arguments

• grad_greenfunction: Gradient Green's function matrix (output)
• greenfunction: Green's function matrix (output)
• observer: Observer geometry struct (PlasmaGeometry or WallGeometry)
• source: Source geometry struct (PlasmaGeometry or WallGeometry)
• n: Toroidal mode number

Returns

Modifies grad_greenfunction and greenfunction in place. Note that greenfunction is zeroed only when the source is plasma; grad_greenfunction is not zeroed since it fills a different block of the (2 * mtheta, 2 * mtheta) depending on the source/observer.

Notes

• Uses Simpson's rule for integration away from singular points
• Uses Gaussian quadrature near singular points for improved accuracy
• Implements analytical singularity removal [Chance Phys. Plasmas 1997 2161]

compute_3D_kernel_matrices!(grad_greenfunction, greenfunction, observer, source, PATCH_RAD, RAD_DIM, INTERP_ORDER)

Compute boundary integral kernel matrices for 3D geometries with the singular correction algorithm from [Malhotra Plasma Phys. and Cont. Fusion 2019 024004]. Uses multi-threading for parallel computation over observer points.

• Far regions: Rectangle rule with uniform weights (1/N)
• Singular regions: Polar quadrature with partition-of-unity blending

gradgreenfunction is the double-layer kernel matrix, where each entry is ∇{xsrc} φ(xobs, xsrc) · nsrc, and greenfunction is the single-layer kernel matrix, where each entry is φ(xobs, xsrc).

Takes advantage of field periodicity to evaluate the kernel only over a single field period.

Arguments

• grad_greenfunction: Double-layer kernel matrix (Nobs × Nsrc) filled in place
• greenfunction: Single-layer kernel matrix (Nobs × Nsrc); filled only when source is plasma
• observer: Observer geometry (PlasmaGeometry3D)
• source: Source geometry (PlasmaGeometry3D)
• PATCH_RAD: Number of points adjacent to source point to treat as singular
• RAD_DIM: Polar radial quadrature order. Angular order = 2 * RAD_DIM
• INTERP_ORDER: Lagrange interpolation order, must be ≤ (2 * PATCH_RAD + 1)

compute_polar_normal!(n_polar, dr_dθ_polar, dr_dζ_polar, normal_orient)

Compute normal vector (= ∂r/∂θ × ∂r/∂ζ) at polar quadrature points from interpolated tangent vectors. We already scaled the normals by normal_orient in the geometry construction, so we need to reapply that here since we are recomputing the normals from the derivatives.

Arguments

• n_polar: Preallocation unit normal vector at each polar point (RADDIM × ANGDIM × 3)
• dr_dθ_polar: Interpolated ∂r/∂θ at polar points (RADDIM × ANGDIM × 3)
• dr_dζ_polar: Interpolated ∂r/∂ζ at polar points (RADDIM × ANGDIM × 3)
• normal_orient: Multiplier applied to normals to make them orient out of vacuum region (+1 or -1)

compute_vacuum_field(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall::WallGeometry,
       Bn::Vector{<:Number}, R_grid::AbstractVector, Z_grid::AbstractVector)

Calculate the perturbed magnetic field in the vacuum region resulting from a normal magnetic field perturbation (Bn) at the plasma surface. Replaces mscfld from Fortran.

This function orchestrates the vacuum field calculation by:

1. Calling vaccal! to compute the vacuum response kernel (grri)
2. Defining a grid of points (R_grid, Z_grid) where the field is to be calculated
3. Calling _pickup_field to compute the magnetic field components on that grid using the kernel and the source perturbation Bn

Arguments

• inputs::VacuumInput: Struct containing vacuum calculation parameters (n, mpert, mtheta, etc.)
• plasma_surf::PlasmaGeometry: Struct with plasma surface geometry and basis functions
• wall::WallGeometry: Struct with wall geometry
• Bn::Vector{<:Number}: Complex vector of Fourier harmonics of the normal magnetic field perturbation at the plasma surface, B_n = B_n_real + i*B_n_imag. Length must be mpert.
• R_grid::AbstractVector: Vector of R coordinates for the output field grid
• Z_grid::AbstractVector: Vector of Z coordinates for the output field grid

Returns

• B_R::Matrix{ComplexF64}: The R-component of the magnetic field on the grid
• B_Z::Matrix{ComplexF64}: The Z-component of the magnetic field on the grid
• B_phi::Matrix{ComplexF64}: The toroidal component of the magnetic field on the grid
• grid_info::Matrix{Int}: Information about the grid points (1=inside plasma, 0=outside)

compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false)

In-place variant that populates the arrays of an existing VacuumResponse. Dispatches on dimensionality only: 2D (inputs.nzeta == 1) routes to [_compute_vacuum_response_2d!], 3D to [_compute_vacuum_response_3d!].

compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) -> VacuumResponse

Compute the vacuum response for the given inputs. Allocating wrapper around compute_vacuum_response!; pass a preallocated VacuumResponse to that method instead when reusing storage across calls. Pass compute_Iv=true to additionally populate the surface-current matrix I_v (2D only).

cross3!(c, a, b, idx)

Inline function for fast cross product of two 3D vectors at a given index.

Arguments

• c::AbstractMatrix: Output matrix
• a::AbstractMatrix: First input vector
• b::AbstractMatrix: Second input vector
• idx: Index of the output vector

distribute_to_equal_arc_grid(xin, zin)

Given a set of points (xin, zin) that define a closed curve in 2D, redistribute these points to be equally spaced in terms of arc length along the curve. This is useful for ensuring that points on a wall or plasma surface are uniformly distributed in space rather than in the parameter θ. This performs the same function as eqarcw in the Fortran code. We now use FastInterpolations instead of a manual Lagrange interpolation.

Arguments

• xin::Vector{Float64}: x-coordinates of the original points defining the curve (endpoint not included)
• zin::Vector{Float64}: z-coordinates of the original points defining the curve (endpoint not included)

Returns

• xout::Vector{Float64}: x-coordinates of the redistributed points, equally spaced in arc length
• zout::Vector{Float64}: z-coordinates of the redistributed points, equally spaced in arc length

This function is different from elliptic integral E(k). Be careful.

Returns : E(1-m1)

This function is different from elliptic integral K(k). Be careful.

Returns : K(1-m1)

elliptic_integrals_bulirsch(m1; error=1e-8, maxit=10)

Compute complete elliptic integrals K(m1) and E(m1) using Bulirsch's algorithm. This is the Julia equivalent of the Fortran ek3 subroutine.

Arguments

• m1::Float64: Complementary parameter (1 - k²), where k is the elliptic modulus
• error::Float64: Convergence tolerance (default 1e-8)
• maxit::Int: Maximum iterations (default 10)

Returns

• K::Float64: Complete elliptic integral of the first kind K(m1)
• E::Float64: Complete elliptic integral of the second kind E(m1)
• convergence::Float64: Convergence metric
• iterations::Int: Number of iterations performed

Notes

• Based on Bulirsch's method as described in Numerical Recipes
• Precision is approximately error²
• Reference: JCP 221 (2007) 330-348

expand_field_periods(inputs::VacuumInput) -> VacuumInput

Reconstruct the full-torus boundary from a single field period using the device field periodicity. When inputs.nfp > 1 (and the boundary is 3D), the input x, y, z arrays describe one field period spanning ζ ∈ [0, 2π/nfp) on an mtheta_in × nzeta_in grid; the full torus is generated by nfp rigid rotations about the Z-axis by 2π·k/nfp. The toroidal counts nzeta_in and nzeta are scaled to their full-torus totals and nfp is reset to 1, so all downstream code operates on the full surface unchanged.

This is the layer-1 interface: it shrinks the Python→Julia payload to one field period while producing a full-torus geometry identical (to floating-point) to sending the whole torus.

Returns inputs unchanged for axisymmetric (nzeta_in == 1) or single-period (nfp <= 1) cases.

extract_patch!(patch, data, idx_pol_center, idx_tor_center, npol, ntor, PATCH_DIM)

Extract a PATCHDIM × PATCHDIM patch of data centered at (idxpolcenter, idxtorcenter) with periodic wrapping.

Arguments

• patch: Preallocated output array for data around the singular point (PATCHDIM × PATCHDIM × dof)
• data: Source data array (can be coordinates, normals, or area elements)
• idx_pol_center: Center poloidal index
• idx_tor_center: Center toroidal index
• npol: Number of poloidal points
• ntor: Number of toroidal points
• PATCH_DIM: Patch size (must be odd)

extract_plasma_surface_at_psi(equil::Equilibrium.PlasmaEquilibrium, ψ::Float64) -> (r, z, ν)

Extract plasma surface geometry from equilibrium at specified flux coordinate.

Evaluates equilibrium bicubic spline around the flux surface to get R, Z coordinates and computes the toroidal angle offset ν for vacuum calculations.

Arguments

• equil: Equilibrium solution with rzphi bicubic spline
• ψ: Normalized flux coordinate (0 at axis, 1 at edge)

Returns

Tuple of:

• r::Vector{Float64}: R-coordinates around flux surface [mtheta]
• z::Vector{Float64}: Z-coordinates around flux surface [mtheta]
• ν::Vector{Float64}: Toroidal angle offset ν [mtheta]

Implementation

The equilibrium bicubic spline rzphi stores:

• rzphi.f[1] = r² (or rfac²)
• rzphi.f[2] = offset (straight-fieldline poloidal angle offset from geometric poloidal angle)
• rzphi.f[3] = ν (toroidal angle offset from geometric toroidal angle)
• rzphi.f[4] = jac (Jacobian)

From these we compute:

• r_minor = √(rzphi.f[1])
• θcyl = 2π*(θsfl + rzphi.f[2])
• R = R₀ + rminor * cos(θcyl)
• Z = Z₀ + rminor * sin(θcyl)

Reference

[Chance Phys. Plasmas 1997 2161] Matches GPEC's ahgwrite and gpvacuumflxsurf approach (gpvacuum.f line 291-296)

get_pn_quad_cache(n::Int) -> PnQuadEntry

Return cached sinh/cosh values for toroidal mode n, computing on first access. Works for any n ≥ 1 with no upper limit.

The fast path (same n as last call) is lock-free: one atomic read + comparison.

get_singular_quadrature(PATCH_RAD::Int, RAD_DIM::Int, INTERP_ORDER::Int)

Get cached singular quadrature data, initializing if necessary. Returns cached data if parameters match the cached initialization; reinitializes if parameters differ. This allows the user to change quadrature parameters between calls, but prevents redundant reinitialization when parameters are unchanged.

green(x_obs, z_obs, x_source, z_source, dx_dtheta, dz_dtheta, n; gamma_prefactor, uselegacygreenfunction=false)

Compute the Green's function and related quantities for axisymmetric geometry according to equations (36)-(42) of Chance 1997. Replaces green from Fortran code.

Arguments

• x_obs: Observation point R-coordinate (Float64)
• z_obs: Observation point Z-coordinate (Float64)
• x_source: Source point R-coordinate (Float64)
• z_source: Source point Z-coordinate (Float64)
• dx_dtheta: Derivative ∂R'/∂θ at source point (Float64)
• dz_dtheta: Derivative ∂Z'/∂θ at source point (Float64)
• n: Toroidal mode number (Int)
• gamma_prefactor: Precomputed value of 2√π · Γ(1/2 - n) [Chance Phys. Plasmas 1997 eq. 40]. Constant for a given n; callers in tight loops should compute this once and pass it in. Defaults to 2 * sqrt(π) * gamma(0.5 - n) if omitted.
• uselegacygreenfunction::Bool: Flag to use the 1997 version of the Legendre function (default false, uses 2007 version)

Returns

• G_n: 2π𝒢ⁿ(θ,θ′) — Green's function value
• coupling_n: 𝒥 ∇'𝒢ⁿ∇'ℒ — Coupling term for mode n
• coupling_0: 1/(2π) 𝒥 ∇'𝒢⁰∇'ℒ — Coupling term for mode 0

Notes

• Uses Legendre functions P^n_{-1/2}(s) computed via elliptic integrals
• Implements analytical derivatives from Chance 1997 equations
• The coupling terms include the Jacobian factor from the coordinate transformation
• By default uses the 2007 Legendre function implementation (Bulirsch + Gaussian integration)

interpolate_to_polar!(polar_data, patch, quad_data)

Interpolate Cartesian patch data to polar quadrature points using sparse matrix multiply. Overwrites polar_data using mul! function arguments, mul!(C, A, B) -> C where C = A * B.

Arguments

• polar_data: Preallocated output array for polar data (RADDIM × ANGDIM × dof)
• patch: Patch data (PATCHDIM × PATCHDIM × dof)
• P2G: Sparse interpolation matrix

laplace_kernel(ox, oy, oz, sx, sy, sz, nx, ny, nz) -> (single, double)

Fused scalar-argument Laplace kernels for the 3D vacuum BIE.

Returns a tuple (single, double) where:

• single = 1/r is the single-layer kernel
• double = (Δx⋅n)/r^3 is the double-layer kernel

This is used when compute_3D_kernel_matrices! needs both kernels for the same pair, so the distance computation (sqrt(r²)) is shared. Returns (0.0, 0.0) when r² < 1e-30.

periodic_wrap(x, n) -> Int

Inline function for fast periodic wrapping for indices near the valid range [1, n]. Equivalent to mod1(x, n) but avoids division for small offsets. Only valid when x is within one period of the valid range (i.e., 1-n < x < 2n).

precompute_lagrange_stencils(gaussian_points)

Precompute 5-point Lagrange interpolation stencils for Gaussian quadrature points.

Returns a tuple (left, right) where each entry is a Vector of SVector{5,Float64} containing the stencil weights for points on the left/right panel.

compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) -> VacuumResponse

Compute the vacuum response for the given inputs. Allocating wrapper around compute_vacuum_response!; pass a preallocated VacuumResponse to that method instead when reusing storage across calls. Pass compute_Iv=true to additionally populate the surface-current matrix I_v (2D only).

compute_vacuum_field(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall::WallGeometry,
       Bn::Vector{<:Number}, R_grid::AbstractVector, Z_grid::AbstractVector)

Calculate the perturbed magnetic field in the vacuum region resulting from a normal magnetic field perturbation (Bn) at the plasma surface. Replaces mscfld from Fortran.

This function orchestrates the vacuum field calculation by:

1. Calling vaccal! to compute the vacuum response kernel (grri)
2. Defining a grid of points (R_grid, Z_grid) where the field is to be calculated
3. Calling _pickup_field to compute the magnetic field components on that grid using the kernel and the source perturbation Bn

Arguments

• inputs::VacuumInput: Struct containing vacuum calculation parameters (n, mpert, mtheta, etc.)
• plasma_surf::PlasmaGeometry: Struct with plasma surface geometry and basis functions
• wall::WallGeometry: Struct with wall geometry
• Bn::Vector{<:Number}: Complex vector of Fourier harmonics of the normal magnetic field perturbation at the plasma surface, B_n = B_n_real + i*B_n_imag. Length must be mpert.
• R_grid::AbstractVector: Vector of R coordinates for the output field grid
• Z_grid::AbstractVector: Vector of Z coordinates for the output field grid

Returns

• B_R::Matrix{ComplexF64}: The R-component of the magnetic field on the grid
• B_Z::Matrix{ComplexF64}: The Z-component of the magnetic field on the grid
• B_phi::Matrix{ComplexF64}: The toroidal component of the magnetic field on the grid
• grid_info::Matrix{Int}: Information about the grid points (1=inside plasma, 0=outside)