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.

Thread-local scratch for compute_3D_kernel_matrices!. One instance per thread so the parallel observer loop does not race.

KernelWorkspace(PATCH_DIM, RAD_DIM, ANG_DIM)

Preallocated patch, polar, and grid buffers for one observer row.

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 polar singular-correction quadrature. P2G maps polar samples to the Cartesian patch (grid = P2G * polar, polar = P2G' * patch). Gpou/Ppou are the Cartesian and polar partitions of unity; Gpou = -χ.

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, RAD_DIM, INTERP_ORDER)

Build the polar quadrature, partitions of unity, and Lagrange interpolant for the singular patch. ANG_DIM = 2 * RAD_DIM. INTERP_ORDER must be ≤ 2 * PATCH_RAD + 1.

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

StellSymBasis

Unitary change of basis U for residue class k making U†D̂ₖU real, and block-diagonal when the class is self-conjugate. Both Laplace kernels depend only on |r_obs − r_src| and n_src·(r_obs − r_src), so a stellarator-symmetric boundary gives

D̂ₖ[σp, σq] = ω^{k(a_p − a_q)} · conj(D̂ₖ[p, q]),   ω = exp(-2πi/nfp),

with σ the within-period reflection (θ, ζ) → (−θ, −ζ) and a_p = 0 on the ζ = 0 plane, 1 elsewhere. Pairing p with σp into a symmetric and an antisymmetric column makes the result real; a self-conjugate class (is_self_conjugate) has real D̂ₖ and splits by the sign of ω^{k a_p}.

Fields

• σ_map: σ as a grid-index map of 1:num_points_per_fp
• σ_phase: ω^{k a_p}, relating a pair's two operator rows
• pair_reps: one grid index per reflection pair (partner is σ_map[p])
• columns, pair_columns: basis columns, and the range each pair contributes
• block_sizes: columns per surface in each block (length is 1 or 2)

StellSymBasis(σ_map, mtheta, k, nfp)

Build the symmetry-adapted basis for toroidal residue class k.

Column cp·e_p + cq·e_q of the symmetry-adapted basis, at slot of block. A grid point the reflection fixes has q == p and cq == 0.

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 and in mod(n, nfp) for 3D
• 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. Stored without the μ₀/4π² normalization: the physical surface inductance is μ₀(2π)²·I_v⁻¹ (see PerturbedEquilibrium.calc_surface_inductance). In 3D it differences two solves, so it needs a finer toroidal grid than wv to converge.
• 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 have mtheta * nzeta rows, ordered with the poloidal index fastest.

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, plasma_surf::PlasmaGeometry3D, wall_settings::WallShapeSettings) -> WallGeometry3D

Constructor to initialize the 3D wall geometry based on the provided vacuum inputs, plasma surface and wall shape settings. This is the 3D counterpart of WallGeometry and selects the shape the same way: an axisymmetric boundary builds the 2D poloidal contour and revolves it, so every shape WallGeometry offers is available; a non-axisymmetric boundary offsets the plasma surface.

Expects a full-torus boundary on the same mtheta/nzeta grid as plasma_surf, and nfp == 1; expand_field_periods produces both from a per-period non-axisymmetric boundary.

Arguments

• inputs::VacuumInput: Struct containing vacuum calculation parameters
• plasma_surf::PlasmaGeometry3D: Plasma surface the wall is built around (used by the conformal wall)
• wall_settings::WallShapeSettings: Struct specifying wall shape and parameters

Returns

• WallGeometry3D: Struct containing wall surface coordinates and derivatives

Notes

• Axisymmetric boundaries (nzeta_in == 1) support nowall, conformal, elliptical, dee, moddee and fromfile; non-axisymmetric boundaries support nowall and conformal
• The non-axisymmetric conformal wall displaces each plasma point along its own normal, so wall grid index (i, j) is the closest wall point to plasma index (i, j) — the correspondence the near-field patch of compute_3D_kernel_matrices! assumes. The revolved branch does not share it: the wall sits on the geometric angle ϕ while PlasmaGeometry3D places the plasma on ζ = ϕ − ν(θ), and equal_arc_wall re-parameterizes poloidally on top of that
• Rejects a conformal offset that folds the surface, and warns when the plasma–wall gap is smaller than one cell of the coarser grid (the double-layer near-field quadrature cannot resolve it). Both fold tests are local, so a global self-intersection between distant lobes — where the offset surface meets itself with aligned normals and a positive area element — passes them

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 non-axisymmetric boundary (nzeta_in > 1) supports "nowall" and "conformal"; the others are poloidal contours that get revolved and so need nzeta_in == 1.

• a::Float64: Distance of wall from plasma in units of the minor radius 0.5(max R - min R) (conformal), or shape parameter (others). On a non-axisymmetric boundary the extrema are taken over the whole torus, so a scales with half the global R-extent rather than a cross-section minor radius; it reduces to the axisymmetric definition when the boundary is axisymmetric.

• 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). Re-parameterizing the contour breaks the plasma/wall grid index correspondence that the singular quadrature relies on; it is ignored outright for nzeta_in > 1 walls.

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

_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, inputs, wall_settings; compute_Iv=false)

3D (inputs.nzeta > 1) vacuum response via block-circulant field-period reduction.

The nfp-periodic boundary makes the layer operators block-circulant in the field-period index, so the problem block-diagonalizes by toroidal residue class k = mod(n, nfp):

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

with D_d, S_d coupling observers in field period 0 to sources in period d. Each class is one solve wv[k] = (4π²/M)·Ẽᴴ·(D̂ₖ \ Ŝₖ)|_plasma·Ẽ, the 3D analogue of the 2D loop over decoupled n. The phase sum is folded into the kernel write, so only the reduced operator is stored.

Real D_d, S_d give D̂₋ₖ = conj(D̂ₖ), so class mod(nfp - k, nfp) reuses this factorization with a conjugated mode basis. A stellarator-symmetric boundary gets a StellSymBasis, making the class operator real and splitting a self-conjugate class into two half-size blocks.

_solve_projected_rhs!(grre, grri, Ẽ, op, scratch)

Form S Ẽᴴ and solve D \ · into grre (exterior) and optionally grri (interior). op = (; green, lu_ext, lu_int) with lu_int === nothing when I_v is not needed. scratch = (; basis_real, rhs_real, rhs_real_int) holds the [Re Im] workspace for a real operator. A real green uses that split because BLAS has no mixed real/complex product.

_unsplit!(dest, split)

Recombine a real [Re Im] column pair into the complex dest.

_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_blocks, green_blocks, observer, source, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym_basis=nothing)

3D single- and double-layer kernels with Malhotra's polar singular correction (PPCF 2019 024004).

Far field: trapezoidal rule. Near field: polar quadrature blended by a partition of unity.

A source in field period d is accumulated onto the period-0 column with weight phases[d+1], so the block-circulant reduction is written in and the full-torus source blocks are never stored. phases = [1.0] is the single-period (real) case.

grad_blocks is ∇_{x_src} φ · n_src; green_blocks is φ, filled only when source is the plasma (Chance 1997 eqs. 26–27). sym_basis evaluates one point per reflection pair and emits the real reduced operator; omit it for the untransformed matrix.

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

n = ∂r/∂θ × ∂r/∂ζ at the polar nodes. Re-apply normal_orient: these normals are rebuilt from interpolated tangents, so they do not inherit the stored surface orientation.

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

Compute the vacuum response into vac_data. 2D (inputs.nzeta == 1) and 3D otherwise. Pass compute_Iv=true to also fill the surface-current matrix I_v.

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

Allocating version of compute_vacuum_response!.

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 a periodically wrapped PATCH_DIM × PATCH_DIM patch centered at (idx_pol_center, idx_tor_center).

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_conjugate_groups(classes, nfp) -> Vector{Vector{Int}}

Each class is a mode family k = mod(n, nfp). Group each with its conjugate family mod(nfp - k, nfp) when that class is present. Self-conjugate families (is_self_conjugate) never pair.

Returns a vector of groups; each group is a vector of indices into classes, with the representative first and its conjugate partner second when paired ([i] or [i, j]).

get_pn_quad_cache(n::Int) -> PnQuadEntry

Return cached sinh/cosh values for toroidal mode n, computing on first access. Takes a lock, so call it once per kernel call and pass the entry down.

get_singular_quadrature(PATCH_RAD, RAD_DIM, INTERP_ORDER)

Return the cached SingularQuadratureData, rebuilding it if the parameters changed.

green(x_obs, z_obs, x_source, z_source, dx_dtheta, dz_dtheta, n; gamma_prefactor, pn_cache, 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.
• pn_cache: get_pn_quad_cache(n); takes a lock, so tight-loop callers look it up once and pass it in.
• 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 a Cartesian patch onto polar quadrature nodes: polar = P2G' * patch.

Whether residue class k is its own conjugate (ω^k = ±1), so D̂ₖ and the field-period phases are real.

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

Fused Laplace kernels: single = 1/r, double = (Δx·n)/r³. Shares √(r²); returns (0, 0) at coincidence (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.

Grid index of (θ, ζ) → (−θ, −ζ) for the linear index p = i_θ + mtheta·(i_ζ − 1).

stell_sym_map(plasma, wall, nfp) -> Union{Vector{Int},Nothing}

Within-period map of (θ, ζ) → (−θ, −ζ) if both surfaces are stellarator symmetric to 1e-10 of the largest coordinate, nothing otherwise. A nearly symmetric boundary (e.g. one extracted from an equilibrium, off by ~1e-3) is rejected and takes the general path; this test is the only guard on the symmetric path.

store_kernel_row!(blocks, sym_basis, kernel_row, scratch, i_pair, row_index, col_index)

Write one accumulated kernel row into blocks. With a StellSymBasis, reconstruct the partner row from the reflection identity and emit U†DU; sym_basis === nothing stores the row as-is. i_pair indexes the observer list (every period-0 point, or one representative per pair). row_index/col_index are 1 = plasma, 2 = wall. scratch needs 2·npts when a basis is present.

transform_mode_basis!(E_out, E, sym_basis; conjugate=false)

Write E_out = E·U into E_out[b] per block. U is unitary, so substituting E_out for E leaves the RHS, the wv projection and I_v unchanged. conjugate=true gives conj(E)·U for this class's conjugate partner. sym_basis === nothing is the identity.

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

Allocating version of compute_vacuum_response!.