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.

KernelParams2D(n::Int)

Parameter struct for 2D vacuum kernel dispatch. Holds the toroidal mode number n.

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

Parameter struct for 3D vacuum kernel dispatch. Holds singular quadrature parameters.

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. Note that the gradient/area elements are scaled by dθ and dζ.

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.

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)
• mlow::Int: Lower poloidal mode number
• mpert::Int: Number of poloidal modes
• nlow::Int: Lower toroidal mode number
• npert::Int: Number of toroidal modes
• 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)
• force_wv_symmetry::Bool: Boolean flag to enforce symmetry in the vacuum response matrix

VacuumInput(
    equil::Equilibrium.PlasmaEquilibrium,
    ψ::Float64,
    mtheta::Int,
    mpert::Int,
    mlow::Int,
    n::Int,
    force_wv_symmetry::Bool = true
) -> 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
• mpert: Number of perturbing poloidal modes
• mlow: Lowest poloidal mode number
• n: Toroidal mode number
• force_wv_symmetry::Bool: Boolean flag to enforce symmetry in the vacuum response matrix (default: true)

Returns

VacuumInput structure ready for computevacuumresponse()

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(inputs::VacuumInput, wall_settings::WallShapeSettings;
    green_only=false)

Compute the vacuum response matrix and both Green's functions using provided vacuum inputs.

Single entry point for vacuum calculations.

• For 3D (inputs.nzeta > 1), computes the full coupled response across all (m,n) modes defined by inputs.(mlow, mpert, nlow, npert).

• For 2D geometry (inputs.nzeta == 1), supports either:

  • single-n (inputs.npert == 1): computes (m,n) response for n = inputs.nlow
  • multi-n (inputs.npert > 1): loops over n = inputs.nlow:(inputs.nlow+inputs.npert-1) and returns blocks of the full response matrices with one block per toroidal mode number.

This is the pure Julia implementation that replaces the Fortran mscvac function. It computes both interior (grri) and exterior (grre) Green's functions for GPEC response calculations.

Arguments

• inputs: VacuumInput struct with mode numbers, grid resolution, and boundary info.
• wall_settings::WallShapeSettings: Wall geometry configuration.
• green_only: If true, skip building the response matrix wv and return zeros for wv and xzpts.

Returns

• wv: Complex vacuum response matrix.

  • 2D single-n: mpert × mpert
  • 2D multi-n: (mpert*npert) × (mpert*npert) (block diagonal)
  • 3D: num_modes × num_modes (full coupled)

• grri: Interior Green's function matrix.

• grre: Exterior Green's function matrix.

• xzpts: Coordinate array (mtheta×4 for 2D, mtheta*nzeta×4 for 3D) [Rplasma, Zplasma, Rwall, Zwall].

_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

kernel!(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 each time this function is called, but grad_greenfunction is not 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).

Arguments

• grad_greenfunction: Double-layer kernel matrix (Nobs × Nsrc) filled in place

• greenfunction: Single-layer kernel matrix (Nobs × Nsrc) filled in place

• observer: Observer geometry (PlasmaGeometry3D)

• source: Source geometry (PlasmaGeometry3D)

• PATCH_RAD: Number of points adjacent to source point to treat as singular

  • Total patch size in # of gridpoints = (2 * PATCHRAD + 1) x (2 * PATCHRAD + 1)

• RAD_DIM: Polar radial quadrature order. Angular order = 2 * RAD_DIM

• INTERP_ORDER: Lagrange interpolation order

  • Must be ≤ (2 * PATCH_RAD + 1)

Threading

This function automatically uses all available threads (Threads.nthreads()). Start Julia with julia -t auto or set JULIA_NUM_THREADS to enable multi-threading.

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

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,
    inputs::VacuumInput,
    wall_settings::WallShapeSettings;
    green_only=false)

In-place variant that computes the vacuum response and directly populates the arrays stored in vac_data.

The vac_data argument is expected to provide the following writable fields with compatible sizes:

• wv::AbstractMatrix{ComplexF64} – vacuum response matrix
• grri::AbstractMatrix{Float64} – interior Green's functions
• grre::AbstractMatrix{Float64} – exterior Green's functions
• plasma_pts::AbstractMatrix{Float64} – plasma surface coordinates
• wall_pts::AbstractMatrix{Float64} – wall surface coordinates

This is designed to work with ForceFreeStates.VacuumData but does not depend on its concrete type (duck-typed on field names only).

compute_vacuum_response(
    inputs::VacuumInput,
    wall_settings::WallShapeSettings;
    green_only=false)

Allocate and return the vacuum response matrix and Green's functions for the given vacuum inputs.

This is a thin allocating wrapper around the in‑place [compute_vacuum_response!] implementation. For performance‑critical paths that already own preallocated storage (e.g. ForceFreeStates.VacuumData), prefer the in‑place method to avoid extra heap allocations.

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

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 + β * C.

Arguments

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

kernel!(grad_greenfunction, greenfunction, observer, source, params::KernelParams3D)

Dispatch wrapper for 3D kernel that forwards to compute_3D_kernel_matrices! with params.

laplace_double_layer(x_obs, x_src, n_src) -> Float64

Evaluate the Laplace double-layer (DxU) kernel between a point and a surface element. Returns 0.0 if the observation point coincides with the source point to avoid singularity. Allocation-free scalar arithmetic is used for maximum performance.

The double-layer kernel K is the normal derivative of the fundamental solution:

K(x_obs, x_src, n_src) = ∇_{x_src} φ · n_src = (x_obs - x_src) · n_src / |x_obs - x_src|³

Arguments

• x_obs: Observation point (3D Cartesian coordinates, any AbstractVector)
• x_src: Source point on surface (3D Cartesian coordinates, any AbstractVector)
• n_src: Outward UNIT normal at source point (must be normalized!, any AbstractVector)

Returns

• Float64: Kernel value K(xobs, xsrc, n_src)

laplace_single_layer(x_obs, x_src) -> Float64

Evaluate the Laplace single-layer (FxU) kernel between two 3D points. Returns 0.0 if the observation point coincides with the source point to avoid singularity.

The single-layer kernel φ is the fundamental solution to Laplace's equation:

φ(x_obs, x_src) = 1 / |x_obs - x_src|

Arguments

• x_obs: Observation point (3D Cartesian coordinates, any AbstractVector)
• x_src: Source point (3D Cartesian coordinates, any AbstractVector)

Returns

• Float64: Kernel value φ(xobs, xsrc)

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;
    green_only=false)

Allocate and return the vacuum response matrix and Green's functions for the given vacuum inputs.

This is a thin allocating wrapper around the in‑place [compute_vacuum_response!] implementation. For performance‑critical paths that already own preallocated storage (e.g. ForceFreeStates.VacuumData), prefer the in‑place method to avoid extra heap allocations.

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)