From f0aaf7202cd5add1d4c32a03f7c0138309ad6ed3 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Wed, 25 Feb 2026 12:00:12 -0500
Subject: [PATCH 01/10] VACUUM - WIP - adding in vtk dump for 3D surfaces

---
 .gitignore              |  1 +
 Project.toml            |  4 ++-
 src/Vacuum/Utilities.jl | 57 +++++++++++++++++++++++++++++++++++++++++
 src/Vacuum/Vacuum.jl    |  8 ++++++
 4 files changed, 69 insertions(+), 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index 873d5d11..446774fe 100644
--- a/.gitignore
+++ b/.gitignore
@@ -19,3 +19,4 @@ docs/build/
 anaconda_projects/
 modovmc
 pestotv
+*.vtu
diff --git a/Project.toml b/Project.toml
index 6438c943..4927f97e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,6 +25,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 AdaptiveArrayPools = "0.2.0"
@@ -32,8 +33,8 @@ DelimitedFiles = "1.9.1"
 DiffEqCallbacks = "4.9.0"
 Documenter = "1.14.1"
 FFTW = "1.9.0"
-FastInterpolations = ">=0.2.13"
 FastGaussQuadrature = "1.1.0"
+FastInterpolations = ">=0.2.13"
 HDF5 = "0.17.2"
 JLD2 = "0.6.3"
 LinearAlgebra = "1.11.0"
@@ -47,3 +48,4 @@ StaticArrays = "1.9.15"
 Statistics = "1.11.0"
 TOML = "1.0.3"
 Test = "1.11.0"
+WriteVTK = "1.21.2"
diff --git a/src/Vacuum/Utilities.jl b/src/Vacuum/Utilities.jl
index de4d220b..91ca09c0 100644
--- a/src/Vacuum/Utilities.jl
+++ b/src/Vacuum/Utilities.jl
@@ -152,3 +152,60 @@ Inline function for fast cross product of two 3D vectors at a given index.
         c[idx, 3] = a1*b2 - a2*b1
     end
 end
+
+"""
+    write_surface_to_vtk(surface_coords::Matrix{Float64}, mtheta::Int, nzeta::Int, filename::String)
+
+Write a 3D toroidal surface to a VTK unstructured grid file for visualization in ParaView.
+
+The surface is represented as a structured grid of quadrilaterals with mtheta poloidal points
+and nzeta toroidal points. The grid is periodic in both directions.
+
+# Arguments
+
+  - `surface_coords::Matrix{Float64}`: Array of shape (mtheta * nzeta, 3) containing (X, Y, Z) coordinates
+  - `mtheta::Int`: Number of poloidal grid points
+  - `nzeta::Int`: Number of toroidal grid points
+  - `filename::String`: Base filename for output (without extension)
+
+# Returns
+
+  - `String`: Full path to the created VTK file
+"""
+function write_surface_to_vtk(surface_coords::Matrix{Float64}, mtheta::Int, nzeta::Int, filename::String)
+    num_points = size(surface_coords, 1)
+    @assert num_points == mtheta * nzeta "Surface coordinates size ($num_points) must equal mtheta * nzeta ($(mtheta * nzeta))"
+
+    # WriteVTK expects points in format (dim, num_points), so transpose
+    points = transpose(surface_coords)  # (3, num_points)
+
+    # Create connectivity for quadrilateral cells
+    # Each cell connects 4 points: (i, j), (i+1, j), (i+1, j+1), (i, j+1)
+    # where i is poloidal index (0:mtheta-1) and j is toroidal index (0:nzeta-1)
+    # Handle periodicity: wrap around at boundaries
+    cells = Vector{MeshCell}()
+
+    for j in 0:(nzeta-1)
+        j_next = mod(j + 1, nzeta)
+        for i in 0:(mtheta-1)
+            i_next = mod(i + 1, mtheta)
+
+            # Indices for the 4 corners of the quadrilateral (1-based indexing)
+            idx1 = i + j * mtheta + 1  # (i, j)
+            idx2 = i_next + j * mtheta + 1  # (i+1, j)
+            idx3 = i_next + j_next * mtheta + 1  # (i+1, j+1)
+            idx4 = i + j_next * mtheta + 1  # (i, j+1)
+
+            # VTK_QUAD cell type
+            push!(cells, MeshCell(VTKCellTypes.VTK_QUAD, [idx1, idx2, idx3, idx4]))
+        end
+    end
+
+    # Create VTK grid
+    vtk = vtk_grid(filename, points, cells)
+
+    # Write the file
+    vtk_save(vtk)
+
+    return "$filename.vtu"
+end
diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl
index 8d0816ee..408ac565 100644
--- a/src/Vacuum/Vacuum.jl
+++ b/src/Vacuum/Vacuum.jl
@@ -6,6 +6,7 @@ using FastGaussQuadrature: gausslegendre
 using StaticArrays: SVector
 using SparseArrays
 using AdaptiveArrayPools
+using WriteVTK
 
 # Import parent modules
 import ..Equilibrium
@@ -152,6 +153,13 @@ It computes both interior (grri) and exterior (grre) Green's functions for GPEC
         if inputs.nzeta > 1 # 3D
             plasma_pts .= plasma_surf.r
             wall_pts .= wall.r
+            # Export 3D surfaces to VTK files
+            plasma_vtk_file = write_surface_to_vtk(plasma_surf.r, inputs.mtheta, inputs.nzeta, "plasma_surface_3D")
+            println("  Plasma surface written to: $plasma_vtk_file")
+            if !wall.nowall
+                wall_vtk_file = write_surface_to_vtk(wall.r, inputs.mtheta, inputs.nzeta, "wall_surface_3D")
+                println("  Wall surface written to: $wall_vtk_file")
+            end
         else # 2D
             @views begin
                 plasma_pts[:, 1] .= plasma_surf.x

From ef5dcae3026ecbd41ebb8ba615c1290d34260bcc Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Thu, 26 Feb 2026 12:52:24 -0500
Subject: [PATCH 02/10] VACUUM - WIP - makign 3D grids use the SFL toroidal
 coordinate

---
 src/Utilities/FourierTransforms.jl | 26 ++++++++++++++------------
 src/Vacuum/DataTypes.jl            | 14 ++++++--------
 src/Vacuum/Vacuum.jl               | 12 +++++-------
 3 files changed, 25 insertions(+), 27 deletions(-)

diff --git a/src/Utilities/FourierTransforms.jl b/src/Utilities/FourierTransforms.jl
index 49318060..73c8c148 100644
--- a/src/Utilities/FourierTransforms.jl
+++ b/src/Utilities/FourierTransforms.jl
@@ -88,33 +88,35 @@ function compute_fourier_coefficients(
     npert::Int,
     nlow::Int;
     n_2D::Union{Nothing, Int}=nothing,
-    ν::Vector{Float64}=zeros(Float64, mtheta)
+    ν::Union{Nothing, Vector{Float64}}=nothing
 )
-    # Validate inputs
-    @assert length(ν) == mtheta "ν must have length mtheta"
-    @assert !(nzeta > 1 && n_2D !== nothing) "If nzeta > 1, don't set n_2D since we use the full nlow - nhigh range for the toroidal modes"
-    @assert !(n_2D === nothing && nzeta == 1) "If nzeta == 1 (2D), set n_2D to the toroidal mode number"
 
     # Uniform theta grid: [0, 2π)
     θ_grid = range(; start=0, length=mtheta, step=2π/mtheta)
 
-    if !isnothing(n_2D)
+    if nzeta == 1
+        @assert n_2D !== nothing "n_2D must be set for 2D"
+        @assert ν !== nothing "ν must be set for 2D"
+        @assert length(ν) == mtheta "ν must have length mtheta"
+
         # In 2D, we only use one toroidal mode at a time
         # Compute sin(mθ - nν) and cos(mθ - nν)
         sin_mn_basis = sin.((mlow .+ (0:(mpert-1))') .* θ_grid .- n_2D .* ν)
         cos_mn_basis = cos.((mlow .+ (0:(mpert-1))') .* θ_grid .- n_2D .* ν)
-    else # nzeta > 1
+    else # 3D
+        @assert (n_2D === nothing && ν === nothing) "n_2D and ν should be nothing for 3D"
+
         # In 3D, we need to compute the basis for all modes and grid points
-        # Compute sin(mθ - nν - nϕ) and cos(mθ - nν - nϕ)
-        ϕ_grid = range(; start=0, length=nzeta, step=2π/nzeta)
+        # Compute sin(mθ - nζ) and cos(mθ - nζ)
+        ζ_grid = range(; start=0, length=nzeta, step=2π/nzeta)
         sin_mn_basis = zeros(mtheta * nzeta, mpert * npert)
         cos_mn_basis = zeros(mtheta * nzeta, mpert * npert)
         for idx_n in 1:npert, idx_m in 1:mpert
             n = nlow + idx_n - 1
             m = mlow + idx_m - 1
-            for (j, ϕ) in enumerate(ϕ_grid), (i, θ) in enumerate(θ_grid)
-                cos_mn_basis[i+(j-1)*mtheta, idx_m+(idx_n-1)*mpert] = cos(m * θ - n * (ν[i] + ϕ))
-                sin_mn_basis[i+(j-1)*mtheta, idx_m+(idx_n-1)*mpert] = sin(m * θ - n * (ν[i] + ϕ))
+            for (j, ζ) in enumerate(ζ_grid), (i, θ) in enumerate(θ_grid)
+                cos_mn_basis[i+(j-1)*mtheta, idx_m+(idx_n-1)*mpert] = cos(m * θ - n * ζ)
+                sin_mn_basis[i+(j-1)*mtheta, idx_m+(idx_n-1)*mpert] = sin(m * θ - n * ζ)
             end
         end
     end
diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl
index 48d2ae67..f76ce351 100644
--- a/src/Vacuum/DataTypes.jl
+++ b/src/Vacuum/DataTypes.jl
@@ -229,7 +229,6 @@ that the gradient/area elements are scaled by dθ and dζ.
   - `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)
-  - `ν::Vector{Float64}`: Magnetic toroidal angle offset from geometric toroidal angle (same as 2D PlasmaGeometry)
   - `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)
@@ -239,7 +238,6 @@ struct PlasmaGeometry3D
     mtheta::Int
     nzeta::Int
     r::Matrix{Float64}
-    ν::Vector{Float64}
     dr_dθ::Matrix{Float64}
     dr_dζ::Matrix{Float64}
     normal::Matrix{Float64}
@@ -276,7 +274,7 @@ function PlasmaGeometry3D(inputs::VacuumInput)
     dθ = 2π / mtheta
     dζ = 2π / nzeta
     θ_grid = range(; start=0, length=mtheta, step=dθ)
-    ϕ_grid = range(; start=0, length=nzeta, step=dζ)
+    ζ_grid = range(; start=0, length=nzeta, step=dζ)
 
     # Allocate output arrays
     r = zeros(num_points, 3)
@@ -288,14 +286,15 @@ function PlasmaGeometry3D(inputs::VacuumInput)
     surf_2D = PlasmaGeometry(inputs)
 
     # Build 3D surface point-by-point from 2D contour
-    for i in 1:mtheta, (j, ϕ) in enumerate(ϕ_grid)
-        r[i+mtheta*(j-1), 1] = surf_2D.x[i] * cos(ϕ)
-        r[i+mtheta*(j-1), 2] = surf_2D.x[i] * sin(ϕ)
+    for i in 1:mtheta, (j, ζ) in enumerate(ζ_grid)
+        # Our 3D grids are the SFL angle ζ = ϕ - ν
+        r[i+mtheta*(j-1), 1] = surf_2D.x[i] * cos(ζ - surf_2D.ν[i])
+        r[i+mtheta*(j-1), 2] = surf_2D.x[i] * sin(ζ - surf_2D.ν[i])
         r[i+mtheta*(j-1), 3] = surf_2D.z[i]
     end
 
     # Compute tangent vectors and normal vectors via periodic bicubic splines
-    itps = [cubic_interp((θ_grid, ϕ_grid), reshape(r[:, k], mtheta, nzeta); bc=(PeriodicBC(; endpoint=:exclusive), PeriodicBC(; endpoint=:exclusive))) for k in 1:3]
+    itps = [cubic_interp((θ_grid, ζ_grid), reshape(r[:, k], mtheta, nzeta); bc=(PeriodicBC(; endpoint=:exclusive), PeriodicBC(; endpoint=:exclusive))) for k in 1:3]
     for i in 1:mtheta, j in 1:nzeta
         idx = i + (j - 1) * mtheta
         for k in 1:3
@@ -322,7 +321,6 @@ function PlasmaGeometry3D(inputs::VacuumInput)
         mtheta,
         nzeta,
         r,
-        surf_2D.ν,
         dr_dθ,
         dr_dζ,
         normal,
diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl
index 408ac565..6edde279 100644
--- a/src/Vacuum/Vacuum.jl
+++ b/src/Vacuum/Vacuum.jl
@@ -78,7 +78,8 @@ It computes both interior (grri) and exterior (grre) Green's functions for GPEC
     wall = inputs.nzeta > 1 ? WallGeometry3D(inputs, wall_settings) : WallGeometry(inputs, plasma_surf, wall_settings)
 
     # Compute Fourier basis coefficients
-    cos_mn_basis, sin_mn_basis = compute_fourier_coefficients(inputs.mtheta, inputs.mpert, inputs.mlow, inputs.nzeta, inputs.npert, inputs.nlow; n_2D=n_override, ν=plasma_surf.ν)
+    ν = hasproperty(plasma_surf, :ν) ? plasma_surf.ν : nothing
+    cos_mn_basis, sin_mn_basis = compute_fourier_coefficients(inputs.mtheta, inputs.mpert, inputs.mlow, inputs.nzeta, inputs.npert, inputs.nlow; n_2D=n_override, ν=ν)
     num_points_surf, num_modes = size(cos_mn_basis)
 
     # Create kernel parameters structs used to dispatch to the correct kernel
@@ -154,12 +155,9 @@ It computes both interior (grri) and exterior (grre) Green's functions for GPEC
             plasma_pts .= plasma_surf.r
             wall_pts .= wall.r
             # Export 3D surfaces to VTK files
-            plasma_vtk_file = write_surface_to_vtk(plasma_surf.r, inputs.mtheta, inputs.nzeta, "plasma_surface_3D")
-            println("  Plasma surface written to: $plasma_vtk_file")
-            if !wall.nowall
-                wall_vtk_file = write_surface_to_vtk(wall.r, inputs.mtheta, inputs.nzeta, "wall_surface_3D")
-                println("  Wall surface written to: $wall_vtk_file")
-            end
+            println("Writing 3D surfaces to VTK files...")
+            write_surface_to_vtk(plasma_surf.r, inputs.mtheta, inputs.nzeta, "plasma_surface_3D")
+            !wall.nowall && write_surface_to_vtk(wall.r, inputs.mtheta, inputs.nzeta, "wall_surface_3D")
         else # 2D
             @views begin
                 plasma_pts[:, 1] .= plasma_surf.x

From 4a65443da44f529f089c0c273215ee2fb6184d14 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Thu, 26 Feb 2026 19:12:20 -0500
Subject: [PATCH 03/10] VACUUM - IMPROVEMENT - adding in an interface for 3D
 equilibria

---
 src/Vacuum/DataTypes.jl | 69 +++++++++++++++++++++++++++++------------
 1 file changed, 49 insertions(+), 20 deletions(-)

diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl
index f76ce351..15d8f2b2 100644
--- a/src/Vacuum/DataTypes.jl
+++ b/src/Vacuum/DataTypes.jl
@@ -2,18 +2,25 @@
     VacuumInput
 
 Struct holding plasma boundary and mode data as provided from ForceFreeStates namelist and computed quantities.
+For an axisymmetric boundary, nzeta_in = 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,
+nzeta_in > 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
 
-  - `r::Vector{Float64}`: Plasma boundary R-coordinate as a function of poloidal angle
-  - `z::Vector{Float64}`: Plasma boundary Z-coordinate as a function of poloidal angle
-  - `ν::Vector{Float64}`: Free parameter in specifying toroidal angle, ϕ = 2πζ + ν(ψ, θ), on theta grid
+  - `x::Vector{Float64}`: Plasma boundary X-coordinate (length mtheta_in * nzeta_in)
+  - `y::Vector{Float64}`: Plasma boundary Y-coordinate (length mtheta_in * nzeta_in)
+  - `z::Vector{Float64}`: Plasma boundary Z-coordinate (length mtheta_in * nzeta_in)
+  - `ν::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
+  - `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
 
 # Notes
@@ -21,13 +28,16 @@ Struct holding plasma boundary and mode data as provided from ForceFreeStates na
   - This is a mutable struct because we need to be able to modify the ν vector in n<0 cases.
 """
 @kwdef mutable struct VacuumInput
-    r::Vector{Float64} = Float64[]
+    x::Vector{Float64} = Float64[]
+    y::Vector{Float64} = Float64[]
     z::Vector{Float64} = Float64[]
     ν::Vector{Float64} = Float64[]
-    mlow::Int = 0
-    mpert::Int = 0
-    nlow::Int = 0
-    npert::Int = 0
+    mtheta_in::Int = 0
+    nzeta_in::Int = 1
+    mlow::Int = 1
+    mpert::Int = 1
+    nlow::Int = 1
+    npert::Int = 1
     mtheta::Int = 1
     nzeta::Int = 1
     force_wv_symmetry::Bool = true
@@ -76,9 +86,10 @@ function VacuumInput(
     r, z, ν = extract_plasma_surface_at_psi(equil, ψ)
 
     return VacuumInput(;
-        r=reverse(r),
+        x=reverse(r),
         z=reverse(z),
         ν=reverse(ν),
+        mtheta_in=length(r),
         mlow=mlow,
         mpert=mpert,
         nlow=nlow,
@@ -203,12 +214,15 @@ We interpolate the input plasma boundary arrays from the inputs struct onto the
 """
 function PlasmaGeometry(inputs::VacuumInput)
 
+    @assert !isempty(inputs.ν) "ν must be specified for 2D calculations"
+    @assert length(inputs.x) == length(inputs.z) == length(inputs.ν) "x, z, and ν must have the same length"
+
     # Interpolate arrays from input onto mtheta grid
-    θ_in = range(0.0, 2π; length=length(inputs.r)) # VacuumInput uses [0, 2π] grid
+    θ_in = range(0.0, 2π; length=inputs.mtheta_in) # VacuumInput uses [0, 2π] grid
     θ_out = range(; start=0, length=inputs.mtheta, step=2π/inputs.mtheta) # VACUUM uses [0, 2π) grid
 
     # Use one-shot API with PeriodicBC
-    x = cubic_interp(θ_in, inputs.r, θ_out; bc=PeriodicBC()) # no endpoint handling needed!
+    x = cubic_interp(θ_in, inputs.x, θ_out; bc=PeriodicBC()) # no endpoint handling needed!
     z = cubic_interp(θ_in, inputs.z, θ_out; bc=PeriodicBC())
     ν = cubic_interp(θ_in, inputs.ν, θ_out; bc=PeriodicBC())
 
@@ -249,6 +263,8 @@ end
 
 Construct a 3D axisymmetric toroidal surface from a 2D poloidal contour.
 
+# TODO: update this
+
 # Algorithm
 
  0. Interpolate 2D arrays onto mtheta grid
@@ -283,14 +299,27 @@ function PlasmaGeometry3D(inputs::VacuumInput)
     normal = zeros(num_points, 3)
 
     # Interpolate arrays from input onto mtheta grid (same as 2D)
-    surf_2D = PlasmaGeometry(inputs)
-
-    # Build 3D surface point-by-point from 2D contour
-    for i in 1:mtheta, (j, ζ) in enumerate(ζ_grid)
-        # Our 3D grids are the SFL angle ζ = ϕ - ν
-        r[i+mtheta*(j-1), 1] = surf_2D.x[i] * cos(ζ - surf_2D.ν[i])
-        r[i+mtheta*(j-1), 2] = surf_2D.x[i] * sin(ζ - surf_2D.ν[i])
-        r[i+mtheta*(j-1), 3] = surf_2D.z[i]
+    if inputs.nzeta_in == 1
+        # Build 3D surface point-by-point from 2D contour
+        surf_2D = PlasmaGeometry(inputs)
+        for i in 1:mtheta, (j, ζ) in enumerate(ζ_grid)
+            # Our 3D grids are the SFL angle ζ = ϕ - ν
+            r[i+mtheta*(j-1), 1] = surf_2D.x[i] * cos(ζ - surf_2D.ν[i])
+            r[i+mtheta*(j-1), 2] = surf_2D.x[i] * sin(ζ - surf_2D.ν[i])
+            r[i+mtheta*(j-1), 3] = surf_2D.z[i]
+        end
+    else
+        # TODO: make this better
+        # Interpolate inputs onto vacuum grid (there's gotta be a better way to do this)
+        θ_in = range(0.0, 2π; length=inputs.mtheta_in)
+        ζ_in = range(0.0, 2π; length=inputs.nzeta_in)
+        θ_flat = repeat(collect(θ_grid); inner=nzeta)
+        ζ_flat = repeat(collect(ζ_grid); outer=mtheta)
+        grid_points = (θ_flat, ζ_flat)
+        for (k, data) in enumerate((inputs.x, inputs.y, inputs.z))
+            itp = cubic_interp((θ_in, ζ_in), reshape(data, inputs.mtheta_in, inputs.nzeta_in); bc=(PeriodicBC(), PeriodicBC()))
+            r[:, k] = itp(grid_points)
+        end
     end
 
     # Compute tangent vectors and normal vectors via periodic bicubic splines

From befa0ad35da904b829e24480257f461e60ff4eb6 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Thu, 5 Mar 2026 14:58:56 -0500
Subject: [PATCH 04/10] VACUUM - BUGFIX - incorrect shape for nonaxisymmetric
 surfs

---
 src/Vacuum/DataTypes.jl | 4 ++--
 src/Vacuum/Kernel3D.jl  | 8 ++++----
 src/Vacuum/Vacuum.jl    | 2 +-
 3 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl
index 15d8f2b2..366132c9 100644
--- a/src/Vacuum/DataTypes.jl
+++ b/src/Vacuum/DataTypes.jl
@@ -313,8 +313,8 @@ function PlasmaGeometry3D(inputs::VacuumInput)
         # Interpolate inputs onto vacuum grid (there's gotta be a better way to do this)
         θ_in = range(0.0, 2π; length=inputs.mtheta_in)
         ζ_in = range(0.0, 2π; length=inputs.nzeta_in)
-        θ_flat = repeat(collect(θ_grid); inner=nzeta)
-        ζ_flat = repeat(collect(ζ_grid); outer=mtheta)
+        θ_flat = repeat(collect(θ_grid), nzeta)
+        ζ_flat = repeat(collect(ζ_grid); inner=mtheta)
         grid_points = (θ_flat, ζ_flat)
         for (k, data) in enumerate((inputs.x, inputs.y, inputs.z))
             itp = cubic_interp((θ_in, ζ_in), reshape(data, inputs.mtheta_in, inputs.nzeta_in); bc=(PeriodicBC(), PeriodicBC()))
diff --git a/src/Vacuum/Kernel3D.jl b/src/Vacuum/Kernel3D.jl
index bf49966d..2511be86 100644
--- a/src/Vacuum/Kernel3D.jl
+++ b/src/Vacuum/Kernel3D.jl
@@ -422,11 +422,11 @@ function compute_3D_kernel_matrices!(
 
     # Initialize quadrature data
     # This allows the code to run at lower resolution without erroring out, but will warn the user.
-    if PATCH_RAD > (observer.nzeta - 1) ÷ 2
-        warn(
-            "PATCH_RAD is greater than the number of points in the toroidal direction, which is not supported. Setting PATCH_RAD to $((observer.nzeta - 1) ÷ 2). Double check that you are converged."
+    if PATCH_RAD > (min(observer.nzeta, observer.mtheta) - 1) ÷ 2
+        @warn(
+            "PATCH_RAD is greater than the number of points in the toroidal or poloidal direction, which is not supported. Setting PATCH_RAD to $((min(observer.nzeta, observer.mtheta) - 1) ÷ 2). Double check that you are converged."
         )
-        PATCH_RAD = (observer.nzeta - 1) ÷ 2
+        PATCH_RAD = (min(observer.nzeta, observer.mtheta) - 1) ÷ 2
     end
     quad_data = get_singular_quadrature(PATCH_RAD, RAD_DIM, INTERP_ORDER)
     (; PATCH_DIM, PATCH_RAD, ANG_DIM, RAD_DIM, Ppou, Gpou, P2G) = quad_data
diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl
index 6edde279..8a177b01 100644
--- a/src/Vacuum/Vacuum.jl
+++ b/src/Vacuum/Vacuum.jl
@@ -156,7 +156,7 @@ It computes both interior (grri) and exterior (grre) Green's functions for GPEC
             wall_pts .= wall.r
             # Export 3D surfaces to VTK files
             println("Writing 3D surfaces to VTK files...")
-            write_surface_to_vtk(plasma_surf.r, inputs.mtheta, inputs.nzeta, "plasma_surface_3D")
+            write_surface_to_vtk(plasma_surf.r, inputs.mtheta, inputs.nzeta, "plasma_surface_mt_$(inputs.mtheta)_nz_$(inputs.nzeta)")
             !wall.nowall && write_surface_to_vtk(wall.r, inputs.mtheta, inputs.nzeta, "wall_surface_3D")
         else # 2D
             @views begin

From 62af08ba5a14cc98ed8085a8b0e98ac7e0a457f3 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Thu, 12 Mar 2026 11:49:59 -0400
Subject: [PATCH 05/10] VACUUM - DOCS - updating PlasmaGeometry3D docstring

---
 src/Vacuum/DataTypes.jl | 56 ++++++++++++++++++++++++++++-------------
 1 file changed, 39 insertions(+), 17 deletions(-)

diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl
index e8244945..46fa4ff2 100644
--- a/src/Vacuum/DataTypes.jl
+++ b/src/Vacuum/DataTypes.jl
@@ -255,28 +255,52 @@ struct PlasmaGeometry3D
 end
 
 """
-    PlasmaGeometry3D(plasma_2d::PlasmaGeometry, nzeta::Int)
+    PlasmaGeometry3D(inputs::VacuumInput)
 
-Construct a 3D axisymmetric toroidal surface from a 2D poloidal contour.
+Construct a 3D toroidal plasma surface from vacuum input data.
 
-# TODO: update this
+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`).
 
-# Algorithm
+## Axisymmetric input (inputs.nzeta_in == 1)
 
- 0. Interpolate 2D arrays onto mtheta grid
- 1. Map 2D (R, Z, ν) to 3D Cartesian: X = R cos(ϕ+ν), Y = R sin(ϕ+ν), Z = Z
- 2. Fit periodic bicubic splines to (X, Y, Z) on (θ, ϕ) grid
- 3. Compute tangent vectors via spline gradients
- 4. Compute normals via cross product: n = ∂r/∂θ × ∂r/∂ζ
+ 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).
 
-# Arguments
+## Fully 3D input (inputs.nzeta_in > 1)
 
-  - `plasma_2d`: 2D poloidal plasma geometry
-  - `nzeta`: Number of toroidal grid points
+ 1. Interpolate the input (x, y, z) arrays from the original
+    mtheta_in × nzeta_in grid onto the vacuum mtheta × nzeta grid using
+    periodic bicubic interpolation in both angles.
 
-# Returns
+## 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
+  - `PlasmaGeometry3D`: Complete 3D surface description on the
+    `mtheta × nzeta` grid, including points, tangents, normals, and
+    orientation.
 """
 function PlasmaGeometry3D(inputs::VacuumInput)
 
@@ -294,7 +318,6 @@ function PlasmaGeometry3D(inputs::VacuumInput)
     dr_dζ = zeros(num_points, 3)
     normal = zeros(num_points, 3)
 
-    # Interpolate arrays from input onto mtheta grid (same as 2D)
     if inputs.nzeta_in == 1
         # Build 3D surface point-by-point from 2D contour
         surf_2D = PlasmaGeometry(inputs)
@@ -305,8 +328,7 @@ function PlasmaGeometry3D(inputs::VacuumInput)
             r[i+mtheta*(j-1), 3] = surf_2D.z[i]
         end
     else
-        # TODO: make this better
-        # Interpolate inputs onto vacuum grid (there's gotta be a better way to do this)
+        # Interpolate 3D inputs onto vacuum grid
         θ_in = range(0.0, 2π; length=inputs.mtheta_in)
         ζ_in = range(0.0, 2π; length=inputs.nzeta_in)
         θ_flat = repeat(collect(θ_grid), nzeta)

From 23ef615f3c2400ea0a65cb3268d24c14a8ddf13d Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Thu, 12 Mar 2026 12:15:36 -0400
Subject: [PATCH 06/10] TEST - BUGFIX - fixing vacuum unit tests

---
 src/Utilities/FourierTransforms.jl |  1 -
 test/runtests_vacuum.jl            | 43 ++++++++++++++++++++----------
 2 files changed, 29 insertions(+), 15 deletions(-)

diff --git a/src/Utilities/FourierTransforms.jl b/src/Utilities/FourierTransforms.jl
index f2332e8f..36a53e0b 100644
--- a/src/Utilities/FourierTransforms.jl
+++ b/src/Utilities/FourierTransforms.jl
@@ -125,7 +125,6 @@ function compute_fourier_coefficients(
                     sin_mn_basis[idx, col] = s
                 end
             end
-            end
         end
     end
 
diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl
index a8136a1f..ca694e4e 100644
--- a/test/runtests_vacuum.jl
+++ b/test/runtests_vacuum.jl
@@ -6,13 +6,16 @@
         @testset "VacuumInput" begin
             @testset "default constructor" begin
                 vac = VacuumInput()
-                @test vac.r == Float64[]
+                @test vac.x == Float64[]
+                @test vac.y == Float64[]
                 @test vac.z == Float64[]
                 @test vac.ν == Float64[]
-                @test vac.mlow == 0
-                @test vac.mpert == 0
-                @test vac.nlow == 0
-                @test vac.npert == 0
+                @test vac.mtheta_in == 0
+                @test vac.nzeta_in == 1
+                @test vac.mlow == 1
+                @test vac.mpert == 1
+                @test vac.nlow == 1
+                @test vac.npert == 1
                 @test vac.mtheta == 1
                 @test vac.nzeta == 1
                 @test vac.force_wv_symmetry == true
@@ -50,7 +53,9 @@
         @testset "PlasmaGeometry" begin
             @testset "from VacuumInput" begin
                 inputs = VacuumInput(
-                    r=[1.0, 1.1, 1.2, 1.1, 1.0],
+                    mtheta_in=5,
+                    nzeta_in=1,
+                    x=[1.0, 1.1, 1.2, 1.1, 1.0],
                     z=[0.0, 0.1, 0.0, -0.1, 0.0],
                     ν=zeros(5),
                     mtheta=5,
@@ -72,8 +77,10 @@
             @testset "edge: mtheta larger than input length" begin
                 # Periodic spline requires at least 4 points
                 inputs = VacuumInput(
-                    r=[1.0, 1.1, 1.2, 1.1, 1.0],
+                    x=[1.0, 1.1, 1.2, 1.1, 1.0],
                     z=[0.0, 0.1, 0.0, -0.1, 0.0],
+                    mtheta_in=5,
+                    nzeta_in=1,
                     ν=zeros(5),
                     mtheta=8,
                     mpert=1,
@@ -92,7 +99,9 @@
         # -------------------------------------------------------------------------
         @testset "WallGeometry" begin
             _circle_inputs(mtheta) = VacuumInput(
-                r=1.7 .+ 0.3 .* cos.(range(0, 2π, length=mtheta)),
+                mtheta_in=mtheta,
+                nzeta_in=1,
+                x=1.7 .+ 0.3 .* cos.(range(0, 2π, length=mtheta)),
                 z=0.3 .* sin.(range(0, 2π, length=mtheta)),
                 ν=zeros(mtheta),
                 mtheta=mtheta,
@@ -149,7 +158,9 @@
                 inputs = _circle_inputs(16)
                 plasma_surf_near_zero = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry(
                     VacuumInput(
-                        r=0.05 .+ 0.03 .* cos.(range(0, 2π, length=16)),
+                        mtheta_in=16,
+                        nzeta_in=1,
+                        x=0.05 .+ 0.03 .* cos.(range(0, 2π, length=16)),
                         z=0.03 .* sin.(range(0, 2π, length=16)),
                         ν=zeros(16),
                         mtheta=16,
@@ -170,7 +181,9 @@
                 inputs = _circle_inputs(16)
                 plasma_surf_near_zero = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry(
                     VacuumInput(
-                        r=0.05 .+ 0.03 .* cos.(range(0, 2π, length=16)),
+                        mtheta_in=16,
+                        nzeta_in=1,
+                        x=0.05 .+ 0.03 .* cos.(range(0, 2π, length=16)),
                         z=0.03 .* sin.(range(0, 2π, length=16)),
                         ν=zeros(16),
                         mtheta=16,
@@ -360,7 +373,9 @@
         # -------------------------------------------------------------------------
         @testset "compute_vacuum_response" begin
             _make_inputs(; mtheta=128, mtheta_eq=17, mpert=2, nlow=1, npert=1) = VacuumInput(
-                r=collect(1.7 .+ 0.3 .* cos.(range(0, 2π, length=mtheta_eq))),
+                mtheta_in=mtheta_eq,
+                nzeta_in=1,
+                x=collect(1.7 .+ 0.3 .* cos.(range(0, 2π, length=mtheta_eq))),
                 z=collect(0.3 .* sin.(range(0, 2π, length=mtheta_eq))),
                 ν=zeros(mtheta_eq),
                 mlow=1,
@@ -437,7 +452,9 @@
     # Kernel requires mtheta, nzeta >= PATCH_DIM (23 for default KernelParams3D(11, 20, 5))
     @testset "Vacuum.jl (3D)" begin
         _make_3d_inputs(; mtheta=32, mtheta_eq=17, mpert=2, nlow=0, npert=2, nzeta=32) = VacuumInput(
-            r=collect(1.7 .+ 0.3 .* cos.(range(0, 2π, length=mtheta_eq))),
+            mtheta_in=mtheta_eq,
+            nzeta_in=1,
+            x=collect(1.7 .+ 0.3 .* cos.(range(0, 2π, length=mtheta_eq))),
             z=collect(0.3 .* sin.(range(0, 2π, length=mtheta_eq))),
             ν=zeros(mtheta_eq),
             mlow=1,
@@ -471,8 +488,6 @@
             @test size(surf.dr_dθ) == (num_points, 3)
             @test size(surf.dr_dζ) == (num_points, 3)
             @test size(surf.normal) == (num_points, 3)
-            # ν is from 2D contour: one value per poloidal point (length mtheta)
-            @test length(surf.ν) == inputs.mtheta
             @test surf.normal_orient in (1, -1)
             @test all(isfinite, surf.r)
             @test all(isfinite, surf.normal)

From 43ccc05cb2d6b73ed4649e8c8de9167f85cf27a9 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Thu, 12 Mar 2026 12:24:08 -0400
Subject: [PATCH 07/10] TEST - IMPROVEMENT - addition of nonaxisymmetric unit
 tests

---
 test/runtests_vacuum.jl | 74 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 74 insertions(+)

diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl
index ca694e4e..debb8299 100644
--- a/test/runtests_vacuum.jl
+++ b/test/runtests_vacuum.jl
@@ -465,6 +465,43 @@
             mtheta=mtheta
         )
 
+        # Helper: simple nonaxisymmetric (3D) plasma boundary built from SFL-style (θ, ζ) coordinates.
+        _make_3d_nonaxis_inputs(; mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12, mpert=2, nlow=0, npert=2) = begin
+            θ_in = range(0, 2π, length=mtheta_in)
+            ζ_in = range(0, 2π, length=nzeta_in)
+
+            X = zeros(mtheta_in, nzeta_in)
+            Y = zeros(mtheta_in, nzeta_in)
+            Z = zeros(mtheta_in, nzeta_in)
+
+            R0 = 1.7
+            a = 0.3
+            ε = 0.05
+
+            for (i, θ) in enumerate(θ_in), (j, ζ) in enumerate(ζ_in)
+                # Base circular cross‑section with a small toroidal corrugation
+                R = R0 + a * cos(θ) + ε * cos(2ζ) * cos(θ)
+                Z_ij = 0.3 * sin(θ) + ε * sin(2ζ) * sin(θ)
+                X[i, j] = R * cos(ζ)
+                Y[i, j] = R * sin(ζ)
+                Z[i, j] = Z_ij
+            end
+
+            VacuumInput(
+                x=vec(X),
+                y=vec(Y),
+                z=vec(Z),
+                mtheta_in=mtheta_in,
+                nzeta_in=nzeta_in,
+                mlow=1,
+                mpert=mpert,
+                nlow=nlow,
+                npert=npert,
+                mtheta=mtheta,
+                nzeta=nzeta
+            )
+        end
+
         @testset "VacuumInput nzeta > 1" begin
             vac = VacuumInput(mtheta=32, nzeta=24, mpert=2, npert=2)
             @test vac.nzeta == 24
@@ -496,6 +533,22 @@
             @test !isapprox(surf.r[1, 1], surf.r[1+32, 1]) || !isapprox(surf.r[1, 2], surf.r[1+32, 2])
         end
 
+        @testset "PlasmaGeometry3D nonaxisymmetric input" begin
+            inputs = _make_3d_nonaxis_inputs(mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12)
+            surf = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs)
+            num_points = inputs.mtheta * inputs.nzeta
+
+            @test surf.mtheta == 24
+            @test surf.nzeta == 24
+            @test size(surf.r) == (num_points, 3)
+            @test size(surf.normal) == (num_points, 3)
+            @test surf.normal_orient in (1, -1)
+            @test all(isfinite, surf.r)
+            @test all(isfinite, surf.normal)
+            # Nonaxisymmetric boundary should have genuine 3D structure (non‑trivial Y variation)
+            @test maximum(abs, surf.r[:, 2]) > 0
+        end
+
         @testset "WallGeometry3D nowall" begin
             inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17)
             wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, WallShapeSettings(shape="nowall"))
@@ -541,6 +594,27 @@
             @test isapprox(wv, wv', rtol=1e-12)
         end
 
+        @testset "compute_vacuum_response 3D nonaxisymmetric boundary" begin
+            inputs = _make_3d_nonaxis_inputs(mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12, mpert=2, nlow=0, npert=2)
+            wall_settings = WallShapeSettings(shape="nowall")
+            wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings)
+
+            numpoints = inputs.mtheta * inputs.nzeta
+            num_modes = inputs.mpert * inputs.npert
+
+            @test size(wv) == (num_modes, num_modes)
+            @test eltype(wv) == ComplexF64
+            @test all(isfinite, wv)
+            @test size(grri) == (2 * numpoints, 2 * num_modes)
+            @test size(grre) == (2 * numpoints, 2 * num_modes)
+            @test all(isfinite, grri)
+            @test all(isfinite, grre)
+            @test size(plasma_pts) == (numpoints, 3)
+            @test all(isfinite, plasma_pts)
+            @test size(wall_pts) == (numpoints, 3)
+            @test isapprox(wv, wv', rtol=1e-12)
+        end
+
         @testset "compute_vacuum_response 3D conformal wall" begin
             inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17)
             wall_settings = WallShapeSettings(shape="conformal", a=0.3)

From 0e1ca2ce5cdd3de12b0946469e9f0d2de4f8745e Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Thu, 12 Mar 2026 12:45:02 -0400
Subject: [PATCH 08/10] GPEC - MINOR - change from develop that didn't get
 propagated to 3D

---
 examples/Solovev_ideal_example_3D/gpec.toml | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml
index 68de3c34..d9e8b663 100644
--- a/examples/Solovev_ideal_example_3D/gpec.toml
+++ b/examples/Solovev_ideal_example_3D/gpec.toml
@@ -48,9 +48,7 @@ ktw = 50.0                    # Parameter for ktanh_flag (default: 50.0)
 ion_flag = true               # Include ion dW_k when kin_flag is true
 electron_flag = false         # Include electron dW_k when kin_flag is true
 
-tol_nr = 1e-6                 # Relative tolerance of dynamic integration steps away from rationals
-tol_r = 1e-7                  # Relative tolerance of dynamic integration steps near rationals
-crossover = 1e-2              # Fractional distance from rational q at which tol switches
+eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations
 singfac_min = 1e-4            # Fractional distance from rational q at which ideal jump enforced
 ucrit = 1e3                   # Maximum fraction of solutions allowed before re-normalized
 force_wv_symmetry = true      # Forces vacuum energy matrix symmetry

From dc4e75accde67c3c02000b5cf66855d43d8fda81 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Thu, 12 Mar 2026 12:48:29 -0400
Subject: [PATCH 09/10] VACUUM - MINOR - removing VTK writing for a future PR

---
 .gitignore              |  1 -
 Project.toml            |  2 --
 src/Vacuum/Utilities.jl | 57 -----------------------------------------
 src/Vacuum/Vacuum.jl    |  7 ++---
 4 files changed, 2 insertions(+), 65 deletions(-)

diff --git a/.gitignore b/.gitignore
index 446774fe..873d5d11 100644
--- a/.gitignore
+++ b/.gitignore
@@ -19,4 +19,3 @@ docs/build/
 anaconda_projects/
 modovmc
 pestotv
-*.vtu
diff --git a/Project.toml b/Project.toml
index 7508741a..7d74e308 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,7 +27,6 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 julia = "1.10"
@@ -53,4 +52,3 @@ StaticArrays = "1.9.15"
 Statistics = "1"
 TOML = "1"
 Test = "1"
-WriteVTK = "1.21.2"
diff --git a/src/Vacuum/Utilities.jl b/src/Vacuum/Utilities.jl
index 91ca09c0..de4d220b 100644
--- a/src/Vacuum/Utilities.jl
+++ b/src/Vacuum/Utilities.jl
@@ -152,60 +152,3 @@ Inline function for fast cross product of two 3D vectors at a given index.
         c[idx, 3] = a1*b2 - a2*b1
     end
 end
-
-"""
-    write_surface_to_vtk(surface_coords::Matrix{Float64}, mtheta::Int, nzeta::Int, filename::String)
-
-Write a 3D toroidal surface to a VTK unstructured grid file for visualization in ParaView.
-
-The surface is represented as a structured grid of quadrilaterals with mtheta poloidal points
-and nzeta toroidal points. The grid is periodic in both directions.
-
-# Arguments
-
-  - `surface_coords::Matrix{Float64}`: Array of shape (mtheta * nzeta, 3) containing (X, Y, Z) coordinates
-  - `mtheta::Int`: Number of poloidal grid points
-  - `nzeta::Int`: Number of toroidal grid points
-  - `filename::String`: Base filename for output (without extension)
-
-# Returns
-
-  - `String`: Full path to the created VTK file
-"""
-function write_surface_to_vtk(surface_coords::Matrix{Float64}, mtheta::Int, nzeta::Int, filename::String)
-    num_points = size(surface_coords, 1)
-    @assert num_points == mtheta * nzeta "Surface coordinates size ($num_points) must equal mtheta * nzeta ($(mtheta * nzeta))"
-
-    # WriteVTK expects points in format (dim, num_points), so transpose
-    points = transpose(surface_coords)  # (3, num_points)
-
-    # Create connectivity for quadrilateral cells
-    # Each cell connects 4 points: (i, j), (i+1, j), (i+1, j+1), (i, j+1)
-    # where i is poloidal index (0:mtheta-1) and j is toroidal index (0:nzeta-1)
-    # Handle periodicity: wrap around at boundaries
-    cells = Vector{MeshCell}()
-
-    for j in 0:(nzeta-1)
-        j_next = mod(j + 1, nzeta)
-        for i in 0:(mtheta-1)
-            i_next = mod(i + 1, mtheta)
-
-            # Indices for the 4 corners of the quadrilateral (1-based indexing)
-            idx1 = i + j * mtheta + 1  # (i, j)
-            idx2 = i_next + j * mtheta + 1  # (i+1, j)
-            idx3 = i_next + j_next * mtheta + 1  # (i+1, j+1)
-            idx4 = i + j_next * mtheta + 1  # (i, j+1)
-
-            # VTK_QUAD cell type
-            push!(cells, MeshCell(VTKCellTypes.VTK_QUAD, [idx1, idx2, idx3, idx4]))
-        end
-    end
-
-    # Create VTK grid
-    vtk = vtk_grid(filename, points, cells)
-
-    # Write the file
-    vtk_save(vtk)
-
-    return "$filename.vtu"
-end
diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl
index de1dd2de..2f6335d6 100644
--- a/src/Vacuum/Vacuum.jl
+++ b/src/Vacuum/Vacuum.jl
@@ -6,7 +6,6 @@ using FastGaussQuadrature: gausslegendre
 using StaticArrays: SVector
 using SparseArrays
 using AdaptiveArrayPools
-using WriteVTK
 
 # Import parent modules
 import ..Equilibrium
@@ -58,7 +57,9 @@ It computes both interior (grri) and exterior (grre) Green's functions for GPEC
       + 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) [R_plasma, Z_plasma, R_wall, Z_wall].
 """
 @with_pool pool function _compute_vacuum_response_single!(
@@ -154,10 +155,6 @@ It computes both interior (grri) and exterior (grre) Green's functions for GPEC
         if inputs.nzeta > 1 # 3D
             plasma_pts .= plasma_surf.r
             wall_pts .= wall.r
-            # Export 3D surfaces to VTK files
-            println("Writing 3D surfaces to VTK files...")
-            write_surface_to_vtk(plasma_surf.r, inputs.mtheta, inputs.nzeta, "plasma_surface_mt_$(inputs.mtheta)_nz_$(inputs.nzeta)")
-            !wall.nowall && write_surface_to_vtk(wall.r, inputs.mtheta, inputs.nzeta, "wall_surface_3D")
         else # 2D
             @views begin
                 plasma_pts[:, 1] .= plasma_surf.x

From 1fd10c8111abe231a998e003a025365132fb66bf Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Thu, 12 Mar 2026 13:04:48 -0400
Subject: [PATCH 10/10] VACUUM - MINOR - claude suggestions

---
 src/Vacuum/DataTypes.jl | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl
index 46fa4ff2..63d029b4 100644
--- a/src/Vacuum/DataTypes.jl
+++ b/src/Vacuum/DataTypes.jl
@@ -113,11 +113,19 @@ Struct containing input settings for vacuum wall geometry.
         + `"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)
 """
@@ -277,7 +285,8 @@ struct, handling both axisymmetric (2D boundary, `nzeta_in == 1`) and fully
 
  1. Interpolate the input (x, y, z) arrays from the original
     mtheta_in × nzeta_in grid onto the vacuum mtheta × nzeta grid using
-    periodic bicubic interpolation in both angles.
+    periodic bicubic interpolation in both angles. The inputs are assumed
+    to already be equally spaced on the SFL angle grid.
 
 ## Steps
 
@@ -329,6 +338,8 @@ function PlasmaGeometry3D(inputs::VacuumInput)
         end
     else
         # Interpolate 3D inputs onto vacuum grid
+        @assert !isempty(inputs.y) "y must be provided for 3D (nzeta_in > 1) inputs"
+        @assert length(inputs.x) == length(inputs.y) == length(inputs.z) == inputs.mtheta_in * inputs.nzeta_in "x, y, z must have length mtheta_in * nzeta_in"
         θ_in = range(0.0, 2π; length=inputs.mtheta_in)
         ζ_in = range(0.0, 2π; length=inputs.nzeta_in)
         θ_flat = repeat(collect(θ_grid), nzeta)
@@ -604,6 +615,8 @@ function WallGeometry3D(inputs::VacuumInput, wall_settings::WallShapeSettings)
         )
     end
 
+    inputs.nzeta_in > 1 && error("3D wall geometry not yet implemented for non-axisymmetric inputs")
+
     # Plasma surface coordinates (2D)
     surf_2D = PlasmaGeometry(inputs)
     wall_2D = WallGeometry(inputs, surf_2D, wall_settings)