diff --git a/CLAUDE.md b/CLAUDE.md
index 272e784d..042c77be 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -78,8 +78,8 @@ JPEC consists of four main modules, each organized in `src/`:
    - `Vacuum_math.jl` - Mathematical utilities
    - `fortran/` - Legacy Fortran code compiled to shared library
 
-4. **DCON** (`src/DCON/`) - Stability analysis
-   - `DconStructs.jl` - Core data structures
+4. **ForceFreeStates** (`src/ForceFreeStates/`) - Stability analysis
+   - `ForceFreeStatesStructs.jl` - Core data structures
    - `Main.jl` - Main entry points
    - `Ode.jl` - ODE solver for stability equations
    - `Sing.jl` - Singular point handling
@@ -126,13 +126,13 @@ CODE - TAG - Detailed message
 ```
 
 Where:
-- CODE: Module name (EQUIL, DCON, VAC, VACUUM, etc.)
+- CODE: Module name (EQUIL, ForceFreeStates, VAC, VACUUM, etc.)
 - TAG: Type descriptor (WIP, MINOR, IMPROVEMENT, BUG FIX, NEW FEATURE, etc.)
 
 Examples:
 - `VAC - WIP - Converting vaccal wall code to Julia`
 - `EQUIL - BUG FIX - Fixed separatrix finding for high kappa`
-- `DCON - NEW FEATURE - Added Mercier criterion calculation`
+- `ForceFreeStates - NEW FEATURE - Added Mercier criterion calculation`
 
 This format is used for compiling release notes, so tags should be human-readable and descriptive.
 
diff --git a/README.md b/README.md
index b15a9b58..f33334d1 100644
--- a/README.md
+++ b/README.md
@@ -30,7 +30,7 @@ To assist with the process of compiling release notes, please confirm to the com
 ```
 CODE - TAG - Detailed message
 ```
-where CODE is EQUIL, DCON, VAC, etc. and TAGs are short descriptors of the type of commit. Example tags might be WIP (work in progress), MINOR, IMPROVEMENT, BUG FIX, NEW FEATURE, etc. Look through old commits for examples of what tags are commonly used.
+where CODE is EQUIL, ForceFreeStates, VAC, etc. and TAGs are short descriptors of the type of commit. Example tags might be WIP (work in progress), MINOR, IMPROVEMENT, BUG FIX, NEW FEATURE, etc. Look through old commits for examples of what tags are commonly used.
 
 Again, this is currently used for the by-hand compilation of release notes. The tags thus need to be human readable but are not strictly inforced to be within some limited set. The objective is to allow a lead developer to skim through commits and pick out only the key new features / bug fixes / improvements to note in the release while not having to read all the work-in-progress or minor changes.
 
diff --git a/benchmarks/DIIID_ideal_example/dcon.toml b/benchmarks/DIIID_ideal_example/ForceFreeStates.toml
similarity index 99%
rename from benchmarks/DIIID_ideal_example/dcon.toml
rename to benchmarks/DIIID_ideal_example/ForceFreeStates.toml
index b332bab6..71bad754 100644
--- a/benchmarks/DIIID_ideal_example/dcon.toml
+++ b/benchmarks/DIIID_ideal_example/ForceFreeStates.toml
@@ -1,4 +1,4 @@
-[DCON_CONTROL]
+[ForceFreeStates_CONTROL]
 bal_flag = false              # Ideal MHD ballooning criterion for short wavelengths
 mat_flag = true               # Construct coefficient matrices for diagnostic purposes
 ode_flag = true              # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
diff --git a/benchmarks/DIIID_ideal_example/euler_fortran_julia_comparison.ipynb b/benchmarks/DIIID_ideal_example/euler_fortran_julia_comparison.ipynb
index 427e347b..6f43da85 100644
--- a/benchmarks/DIIID_ideal_example/euler_fortran_julia_comparison.ipynb
+++ b/benchmarks/DIIID_ideal_example/euler_fortran_julia_comparison.ipynb
@@ -328,7 +328,7 @@
     "# Run DCON in Julia\n",
     "Pkg.activate(\"../..\")\n",
     "using JPEC\n",
-    "JPEC.DCON.Main(\"./\") # \"./\" tells us to obtain inputs and direct outputs to our current folder"
+    "JPEC.ForceFreeStates.Main(\"./\") # \"./\" tells us to obtain inputs and direct outputs to our current folder"
    ]
   },
   {
diff --git a/benchmarks/DIIID_ideal_example/vacuum_accuracy_benchmark.jl b/benchmarks/DIIID_ideal_example/vacuum_accuracy_benchmark.jl
index de4a1caf..5068b531 100644
--- a/benchmarks/DIIID_ideal_example/vacuum_accuracy_benchmark.jl
+++ b/benchmarks/DIIID_ideal_example/vacuum_accuracy_benchmark.jl
@@ -44,7 +44,7 @@ if benchmark_n
         println("Testing n = $n_test")
 
         # Compute Fortran solution
-        JPEC.Vacuum.set_dcon_params(mtheta_eq, vac_inputs.mlow, vac_inputs.mhigh,
+        JPEC.Vacuum.set_surface_params(mtheta_eq, vac_inputs.mlow, vac_inputs.mhigh,
                                     n_test, vac_inputs.qa, vac_inputs.r,
                                     vac_inputs.z, vac_inputs.delta)
 
@@ -53,7 +53,7 @@ if benchmark_n
                           complex_flag, kernelsign, wall_flag,
                           farwall_flag, grri, xzpts)
 
-        JPEC.Vacuum.unset_dcon_params()
+        JPEC.Vacuum.unset_surface_params()
 
         # Compute Julia solution
         vac_inputs_test = deepcopy(vac_inputs)
@@ -123,7 +123,7 @@ if benchmark_m
 
     # Fortran reference
     mpert_ref = mhigh_reference - vac_inputs.mlow
-    JPEC.Vacuum.set_dcon_params(mtheta_eq, vac_inputs.mlow, mhigh_reference,
+    JPEC.Vacuum.set_surface_params(mtheta_eq, vac_inputs.mlow, mhigh_reference,
                                 1, vac_inputs.qa, vac_inputs.r,
                                 vac_inputs.z, vac_inputs.delta)
     wv_block_ref_fortran_m = zeros(ComplexF64, mpert_ref, mpert_ref)
@@ -131,7 +131,7 @@ if benchmark_m
     JPEC.Vacuum.mscvac(wv_block_ref_fortran_m, mpert_ref, mtheta_eq, mthvac,
                        complex_flag, kernelsign, wall_flag,
                        farwall_flag, grri_ref, xzpts, ahg_file, "$(@__DIR__)/../../examples/DIIID-like_ideal_example")
-    JPEC.Vacuum.unset_dcon_params()
+    JPEC.Vacuum.unset_surface_params()
 
     println("Reference solutions computed.\n")
 
@@ -142,7 +142,7 @@ if benchmark_m
         mpert_test = mhigh_test - vac_inputs.mlow
 
         # Compute Fortran solution
-        JPEC.Vacuum.set_dcon_params(mtheta_eq, vac_inputs.mlow, mhigh_test,
+        JPEC.Vacuum.set_surface_params(mtheta_eq, vac_inputs.mlow, mhigh_test,
                                     1, vac_inputs.qa, vac_inputs.r,
                                     vac_inputs.z, vac_inputs.delta)
 
@@ -152,7 +152,7 @@ if benchmark_m
                           complex_flag, kernelsign, wall_flag,
                           farwall_flag, grri_test, xzpts, ahg_file, "$(@__DIR__)/../../examples/DIIID-like_ideal_example")
 
-        JPEC.Vacuum.unset_dcon_params()
+        JPEC.Vacuum.unset_surface_params()
 
         # Compute Julia solution
         vac_inputs_test_m = deepcopy(vac_inputs)
diff --git a/benchmarks/DIIID_ideal_example/vacuum_speed_benchmark.jl b/benchmarks/DIIID_ideal_example/vacuum_speed_benchmark.jl
index 04726363..9c46a527 100644
--- a/benchmarks/DIIID_ideal_example/vacuum_speed_benchmark.jl
+++ b/benchmarks/DIIID_ideal_example/vacuum_speed_benchmark.jl
@@ -36,7 +36,7 @@ if benchmark_n
         println("Benchmarking n = $n_test")
         
         # Benchmark Fortran version
-        JPEC.Vacuum.set_dcon_params(mtheta_eq, vac_inputs.mlow, vac_inputs.mhigh, 
+        JPEC.Vacuum.set_surface_params(mtheta_eq, vac_inputs.mlow, vac_inputs.mhigh, 
                                     n_test, vac_inputs.qa, vac_inputs.r, 
                                     vac_inputs.z, vac_inputs.delta)
         
@@ -49,7 +49,7 @@ if benchmark_n
         push!(fortran_allocs, median(b_fortran).allocs)
         push!(fortran_memory, median(b_fortran).memory / 1e6)
         
-        JPEC.Vacuum.unset_dcon_params()
+        JPEC.Vacuum.unset_surface_params()
         
         # Benchmark Julia version
         vac_inputs_test = deepcopy(vac_inputs)
@@ -106,7 +106,7 @@ if benchmark_m
         mpert_test = mhigh_test - vac_inputs.mlow
         
         # Benchmark Fortran version
-        JPEC.Vacuum.set_dcon_params(mtheta_eq, vac_inputs.mlow, mhigh_test, 
+        JPEC.Vacuum.set_surface_params(mtheta_eq, vac_inputs.mlow, mhigh_test, 
                                     1, vac_inputs.qa, vac_inputs.r, 
                                     vac_inputs.z, vac_inputs.delta)
         
@@ -119,7 +119,7 @@ if benchmark_m
         push!(fortran_mhigh_allocs, median(b_fortran).allocs)
         push!(fortran_mhigh_memory, median(b_fortran).memory / 1e6)
         
-        JPEC.Vacuum.unset_dcon_params()
+        JPEC.Vacuum.unset_surface_params()
         
         # Benchmark Julia version
         vac_inputs_test = deepcopy(vac_inputs)
diff --git a/benchmarks/Solovev_ideal_example/vacuum_accuracy_benchmark.jl b/benchmarks/Solovev_ideal_example/vacuum_accuracy_benchmark.jl
index 4db6c7bd..3555cfff 100644
--- a/benchmarks/Solovev_ideal_example/vacuum_accuracy_benchmark.jl
+++ b/benchmarks/Solovev_ideal_example/vacuum_accuracy_benchmark.jl
@@ -48,7 +48,7 @@ if benchmark_n
         println("Testing n = $n_test")
 
         # Compute Fortran solution
-        JPEC.Vacuum.set_dcon_params(mtheta_eq, vac_inputs.mlow, vac_inputs.mhigh,
+        JPEC.Vacuum.set_surface_params(mtheta_eq, vac_inputs.mlow, vac_inputs.mhigh,
                                     n_test, vac_inputs.qa, vac_inputs.r,
                                     vac_inputs.z, vac_inputs.delta)
 
@@ -57,7 +57,7 @@ if benchmark_n
                           complex_flag, kernelsign, wall_flag,
                           farwall_flag, grri, xzpts)
 
-        JPEC.Vacuum.unset_dcon_params()
+        JPEC.Vacuum.unset_surface_params()
 
         # Compute Julia solution
         vac_inputs_test = deepcopy(vac_inputs)
@@ -122,7 +122,7 @@ if benchmark_a
         println("Testing a = $a_test")
 
         # Compute Fortran solution
-        JPEC.Vacuum.set_dcon_params(mtheta_eq, vac_inputs.mlow, vac_inputs.mhigh,
+        JPEC.Vacuum.set_surface_params(mtheta_eq, vac_inputs.mlow, vac_inputs.mhigh,
                                     vac_inputs.n, vac_inputs.qa, vac_inputs.r,
                                     vac_inputs.z, vac_inputs.delta)
 
@@ -136,7 +136,7 @@ if benchmark_a
                           complex_flag, kernelsign, wall_flag,
                           farwall_flag, grri, xzpts)
 
-        JPEC.Vacuum.unset_dcon_params()
+        JPEC.Vacuum.unset_surface_params()
 
         # Compute Julia solution
         new_wall = JPEC.Vacuum.WallShapeSettings(
diff --git a/docs/src/equilibrium.md b/docs/src/equilibrium.md
index ab4b486c..7feb00bc 100644
--- a/docs/src/equilibrium.md
+++ b/docs/src/equilibrium.md
@@ -19,7 +19,7 @@ Key responsibilities of the module:
 	(global parameters, q-profile finding, separatrix finding, GSE checks).
 
 The module exposes a small public API that covers setup, configuration,
-and common analyses used by other JPEC components (e.g. `DCON`, vacuum
+and common analyses used by other JPEC components (e.g. `ForceFreeStates`, vacuum
 interfaces).
 
 ## API Reference
@@ -65,7 +65,7 @@ Basic example: read a TOML config and build an equilibrium
 using JPEC
 
 # Build from a TOML file (searches relative paths if needed)
-pe = JPEC.Equilibrium.setup_equilibrium("docs/examples/dcon.toml")
+pe = JPEC.Equilibrium.setup_equilibrium("docs/examples/ForceFreeStates.toml")
 
 println("Magnetic axis: ", pe.params.r0, ", ", pe.params.z0)
 println("q(0) = ", pe.params.q0)
diff --git a/docs/src/examples/vacuum.md b/docs/src/examples/vacuum.md
index 0cad1192..3351253e 100644
--- a/docs/src/examples/vacuum.md
+++ b/docs/src/examples/vacuum.md
@@ -4,14 +4,14 @@ This page demonstrates the usage of the Vacuum module for magnetostatic calculat
 
 ## Basic Vacuum Field Calculation
 
-### Setting up DCON Parameters
+### Setting up Surface Parameters
 
-The vacuum module requires initialization with DCON (Displacement CONtinuum) parameters:
+The vacuum module requires initialization with surface parameters:
 
 ```julia
 using JPEC
 
-# Define DCON parameters
+# Define surface parameters
 mthin = Int32(4)     # Number of theta points
 lmin = Int32(1)      # Minimum poloidal mode number
 lmax = Int32(4)      # Maximum poloidal mode number
@@ -24,8 +24,8 @@ xin = rand(Float64, n_modes)      # Radial coordinates
 zin = rand(Float64, n_modes)      # Vertical coordinates
 deltain = rand(Float64, n_modes)  # Displacement data
 
-# Initialize DCON interface
-JPEC.VacuumMod.set_dcon_params(mthin, lmin, lmax, nnin, qa1in, xin, zin, deltain)
+# Initialize surface parameters interface
+JPEC.VacuumMod.set_surface_params(mthin, lmin, lmax, nnin, qa1in, xin, zin, deltain)
 ```
 
 ### Vacuum Matrix Calculation
@@ -137,7 +137,7 @@ println("Most unstable eigenvalue: ", eigenvals[max_growth_idx])
 # xin, zin, deltain = process_boundary_data(boundary_data)
 
 # 4. Perform vacuum calculation
-# JPEC.VacuumMod.set_dcon_params(mthin, lmin, lmax, nnin, qa1in, xin, zin, deltain)
+# JPEC.VacuumMod.set_surface_params(mthin, lmin, lmax, nnin, qa1in, xin, zin, deltain)
 # ... vacuum calculation ...
 
 # 5. Analyze stability
@@ -148,7 +148,7 @@ println("Most unstable eigenvalue: ", eigenvals[max_growth_idx])
 
 ### Common Issues
 
-1. **Initialization Error**: Ensure `set_dcon_params` is called before `mscvac`
+1. **Initialization Error**: Ensure `set_surface_params` is called before `mscvac`
 2. **Memory Issues**: Large `mtheta`/`mthvac` values require significant memory
 3. **Convergence**: Check that geometry arrays are properly normalized
 4. **Complex Arithmetic**: Ensure `complex_flag=true` for stability analysis
diff --git a/docs/src/index.md b/docs/src/index.md
index 4cbd24e7..34e42584 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -46,7 +46,7 @@ To assist with release note compilation, please follow the commit message format
 ```
 CODE - TAG - Detailed message
 ```
-where CODE is EQUIL, DCON, VAC, etc. and TAGs are descriptors like WIP, MINOR, IMPROVEMENT, BUG FIX, NEW FEATURE, etc.
+where CODE is Equil, ForceFreeStates, Vacuum, etc. and TAGs are descriptors like WIP, MINOR, IMPROVEMENT, BUG FIX, NEW FEATURE, etc.
 
 Additionally, please see [this](https://docs.google.com/document/d/1XAOTz1IV8ErZAAk-iSuEuddNOLB5XcoVZsAbPKRUUuA) google doc for more details on using the GitHub.
 
diff --git a/docs/src/vacuum.md b/docs/src/vacuum.md
index 038046fc..44ce0a5f 100644
--- a/docs/src/vacuum.md
+++ b/docs/src/vacuum.md
@@ -7,7 +7,7 @@ Refactored/interfaced from/with VACUUM by M.S. Chance.
 
 The module includes:
 
-- Interface to Fortran vacuum field calculations (`mscvac`, `set_dcon_params`)
+- Interface to Fortran vacuum field calculations (`mscvac`, `set_surface_params`)
 - Pure Julia implementation of vacuum response calculations (`compute_vacuum_response`, `compute_vacuum_field`)
 - Support for various wall geometries and configurations
 
@@ -35,9 +35,9 @@ Modules = [JPEC.Vacuum]
 
 ## Functions
 
-### set_dcon_params
+### set_surface_params
 ```@docs
-JPEC.Vacuum.set_dcon_params
+JPEC.Vacuum.set_surface_params
 ```
 
 ### mscvac
@@ -62,15 +62,15 @@ JPEC.Vacuum.compute_vacuum_field
 ```julia
 using JPEC
 
-# Set DCON parameters
+# Set surface parameters
 mtheta, lmin, lmax, nnin = Int32(4), Int32(1), Int32(4), Int32(2)
 qa1in = 1.23
 xin = rand(Float64, lmax - lmin + 1)
 zin = rand(Float64, lmax - lmin + 1)
 deltain = rand(Float64, lmax - lmin + 1)
 
-# Initialize DCON interface
-JPEC.Vacuum.set_dcon_params(mtheta, lmin, lmax, nnin, qa1in, xin, zin, deltain)
+# Initialize surface parameters interface
+JPEC.Vacuum.set_surface_params(mtheta, lmin, lmax, nnin, qa1in, xin, zin, deltain)
 
 # Set up vacuum calculation parameters
 mpert = 5
@@ -127,7 +127,7 @@ wv, grri, xzpts = JPEC.Vacuum.compute_vacuum_response(inputs, wall_settings)
 
 ## Notes
 
-- Requires proper initialization of DCON parameters before using the Fortran interface
+- Requires proper initialization of surface parameters before using the Fortran interface
 - The pure Julia implementation (`compute_vacuum_response`) provides equivalent functionality
 - For n=0 modes with closed walls, automatic regularization is applied
 
diff --git a/examples/DIIID-like_ideal_example/dcon.toml b/examples/DIIID-like_ideal_example/ForceFreeStates.toml
similarity index 99%
rename from examples/DIIID-like_ideal_example/dcon.toml
rename to examples/DIIID-like_ideal_example/ForceFreeStates.toml
index cb5e2142..3cdb3e74 100644
--- a/examples/DIIID-like_ideal_example/dcon.toml
+++ b/examples/DIIID-like_ideal_example/ForceFreeStates.toml
@@ -1,4 +1,4 @@
-[DCON_CONTROL]
+[ForceFreeStates_CONTROL]
 bal_flag = false              # Ideal MHD ballooning criterion for short wavelengths
 mat_flag = true               # Construct coefficient matrices for diagnostic purposes
 ode_flag = true              # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
diff --git a/examples/DIIID-like_ideal_example/run_and_analyze.ipynb b/examples/DIIID-like_ideal_example/run_and_analyze.ipynb
index 7422391e..713f4022 100644
--- a/examples/DIIID-like_ideal_example/run_and_analyze.ipynb
+++ b/examples/DIIID-like_ideal_example/run_and_analyze.ipynb
@@ -43,7 +43,7 @@
     "# Run DCON in Julia\n",
     "Pkg.activate(\"../..\")\n",
     "using JPEC\n",
-    "JPEC.DCON.Main(\"./\") # \"./\" tells us to obtain inputs and direct outputs to our current folder"
+    "JPEC.ForceFreeStates.Main(\"./\") # \"./\" tells us to obtain inputs and direct outputs to our current folder"
    ]
   },
   {
diff --git a/examples/Solovev_ideal_example/dcon.toml b/examples/Solovev_ideal_example/ForceFreeStates.toml
similarity index 99%
rename from examples/Solovev_ideal_example/dcon.toml
rename to examples/Solovev_ideal_example/ForceFreeStates.toml
index dde9cd83..ac866e92 100644
--- a/examples/Solovev_ideal_example/dcon.toml
+++ b/examples/Solovev_ideal_example/ForceFreeStates.toml
@@ -1,4 +1,4 @@
-[DCON_CONTROL]
+[ForceFreeStates_CONTROL]
 bal_flag = false              # Ideal MHD ballooning criterion for short wavelengths
 mat_flag = true               # Construct coefficient matrices for diagnostic purposes
 ode_flag = true               # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
diff --git a/examples/Solovev_ideal_example/run_and_analyze.ipynb b/examples/Solovev_ideal_example/run_and_analyze.ipynb
index e18c47df..160ccf9b 100644
--- a/examples/Solovev_ideal_example/run_and_analyze.ipynb
+++ b/examples/Solovev_ideal_example/run_and_analyze.ipynb
@@ -43,7 +43,7 @@
     "# Run DCON in Julia\n",
     "Pkg.activate(\"../..\")\n",
     "using JPEC\n",
-    "JPEC.DCON.Main(\"./\") # \"./\" tells us to obtain inputs and direct outputs to our current folder"
+    "JPEC.ForceFreeStates.Main(\"./\") # \"./\" tells us to obtain inputs and direct outputs to our current folder"
    ]
   },
   {
diff --git a/examples/Solovev_ideal_example/vacuum_wall_debug.jl b/examples/Solovev_ideal_example/vacuum_wall_debug.jl
index 5b9fc41a..fe0093c3 100644
--- a/examples/Solovev_ideal_example/vacuum_wall_debug.jl
+++ b/examples/Solovev_ideal_example/vacuum_wall_debug.jl
@@ -32,7 +32,7 @@ if !use_wall_arg
 end
 
 
-JPEC.Vacuum.set_dcon_params(mtheta_eq, mlow, mhigh, 
+JPEC.Vacuum.set_surface_params(mtheta_eq, mlow, mhigh, 
                                     n, vac_inputs.qa, vac_inputs.r, 
                                     vac_inputs.z, vac_inputs.delta)
         
@@ -40,7 +40,7 @@ wv_fortran = JPEC.Vacuum.mscvac(wv_block, mpert, mtheta_eq, mthvac,
                          complex_flag, kernelsign, wall_flag, 
                          farwall_flag, grri, xzpts, ahg_file, "$(@__DIR__)")
 
-JPEC.Vacuum.unset_dcon_params()
+JPEC.Vacuum.unset_surface_params()
 
 # Now deepcopy the julia vac_inputs and set the new ms
 vac_inputs = deepcopy(vac_inputs)
diff --git a/examples/Solovev_ideal_example_multi_n/dcon.toml b/examples/Solovev_ideal_example_multi_n/ForceFreeStates.toml
similarity index 99%
rename from examples/Solovev_ideal_example_multi_n/dcon.toml
rename to examples/Solovev_ideal_example_multi_n/ForceFreeStates.toml
index a2e05d3f..4cca117d 100644
--- a/examples/Solovev_ideal_example_multi_n/dcon.toml
+++ b/examples/Solovev_ideal_example_multi_n/ForceFreeStates.toml
@@ -1,4 +1,4 @@
-[DCON_CONTROL]
+[ForceFreeStates_CONTROL]
 bal_flag = false              # Ideal MHD ballooning criterion for short wavelengths
 mat_flag = true               # Construct coefficient matrices for diagnostic purposes
 ode_flag = true               # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
diff --git a/examples/Solovev_ideal_example_multi_n/run_and_analyze.ipynb b/examples/Solovev_ideal_example_multi_n/run_and_analyze.ipynb
index 695b6d68..58dd0da1 100644
--- a/examples/Solovev_ideal_example_multi_n/run_and_analyze.ipynb
+++ b/examples/Solovev_ideal_example_multi_n/run_and_analyze.ipynb
@@ -42,9 +42,9 @@
     "# Run DCON in Julia\n",
     "Pkg.activate(\"../..\")\n",
     "using JPEC\n",
-    "JPEC.DCON.Main(\"./\") # \"./\" tells us to obtain inputs and direct outputs to our current folder\n",
-    "JPEC.DCON.Main(\"single_n_1\")\n",
-    "JPEC.DCON.Main(\"single_n_2\")"
+    "JPEC.ForceFreeStates.Main(\"./\") # \"./\" tells us to obtain inputs and direct outputs to our current folder\n",
+    "JPEC.ForceFreeStates.Main(\"single_n_1\")\n",
+    "JPEC.ForceFreeStates.Main(\"single_n_2\")"
    ]
   },
   {
diff --git a/examples/Solovev_ideal_example_multi_n/single_n_1/dcon.toml b/examples/Solovev_ideal_example_multi_n/single_n_1/ForceFreeStates.toml
similarity index 99%
rename from examples/Solovev_ideal_example_multi_n/single_n_1/dcon.toml
rename to examples/Solovev_ideal_example_multi_n/single_n_1/ForceFreeStates.toml
index 3ea698b1..42c1ee56 100644
--- a/examples/Solovev_ideal_example_multi_n/single_n_1/dcon.toml
+++ b/examples/Solovev_ideal_example_multi_n/single_n_1/ForceFreeStates.toml
@@ -1,4 +1,4 @@
-[DCON_CONTROL]
+[ForceFreeStates_CONTROL]
 bal_flag = false              # Ideal MHD ballooning criterion for short wavelengths
 mat_flag = true               # Construct coefficient matrices for diagnostic purposes
 ode_flag = true               # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
diff --git a/examples/Solovev_ideal_example_multi_n/single_n_2/dcon.toml b/examples/Solovev_ideal_example_multi_n/single_n_2/ForceFreeStates.toml
similarity index 99%
rename from examples/Solovev_ideal_example_multi_n/single_n_2/dcon.toml
rename to examples/Solovev_ideal_example_multi_n/single_n_2/ForceFreeStates.toml
index 59d2840d..956cd48f 100644
--- a/examples/Solovev_ideal_example_multi_n/single_n_2/dcon.toml
+++ b/examples/Solovev_ideal_example_multi_n/single_n_2/ForceFreeStates.toml
@@ -1,4 +1,4 @@
-[DCON_CONTROL]
+[ForceFreeStates_CONTROL]
 bal_flag = false              # Ideal MHD ballooning criterion for short wavelengths
 mat_flag = true               # Construct coefficient matrices for diagnostic purposes
 ode_flag = true               # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
diff --git a/notebooks/vacuum_fortran_example.ipynb b/notebooks/vacuum_fortran_example.ipynb
index 6f0da82c..ff997391 100644
--- a/notebooks/vacuum_fortran_example.ipynb
+++ b/notebooks/vacuum_fortran_example.ipynb
@@ -32,9 +32,9 @@
     "zin = rand(Float64, lmax - lmin + 1)\n",
     "deltain = rand(Float64, lmax - lmin + 1)\n",
     "\n",
-    "println(\"Testing set_dcon_params…\")\n",
-    "JPEC.VacuumMod.set_dcon_params(mthin, lmin, lmax, nnin, qa1in, xin, zin, deltain)\n",
-    "println(\"set_dcon_params OK!\")"
+    "println(\"Testing set_surface_params…\")\n",
+    "JPEC.VacuumMod.set_surface_params(mthin, lmin, lmax, nnin, qa1in, xin, zin, deltain)\n",
+    "println(\"set_surface_params OK!\")"
    ]
   },
   {
diff --git a/src/DCON/FixedBoundaryStability.jl b/src/ForceFreeStates/FixedBoundaryStability.jl
similarity index 100%
rename from src/DCON/FixedBoundaryStability.jl
rename to src/ForceFreeStates/FixedBoundaryStability.jl
diff --git a/src/DCON/DCON.jl b/src/ForceFreeStates/ForceFreeStates.jl
similarity index 90%
rename from src/DCON/DCON.jl
rename to src/ForceFreeStates/ForceFreeStates.jl
index 59b21be5..ed846378 100644
--- a/src/DCON/DCON.jl
+++ b/src/ForceFreeStates/ForceFreeStates.jl
@@ -1,4 +1,4 @@
-module DCON
+module ForceFreeStates
 
 # All imports and includes for the DCON module
 using LinearAlgebra
@@ -16,8 +16,7 @@ using Printf
 import StaticArrays: @MVector, @MMatrix
 
 # Include all necessary files
-include("DconStructs.jl")
-include("Main.jl")
+include("ForceFreeStatesStructs.jl")
 include("Mercier.jl")
 include("Ode.jl")
 include("Sing.jl")
diff --git a/src/DCON/DconStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl
similarity index 99%
rename from src/DCON/DconStructs.jl
rename to src/ForceFreeStates/ForceFreeStatesStructs.jl
index e3411c39..53d477ab 100644
--- a/src/DCON/DconStructs.jl
+++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl
@@ -55,7 +55,7 @@ A mutable struct containing settings for debugging and benchmarking output.
 end
 
 """
-    DconInternal
+    ForceFreeStatesInternal
 
 A mutable struct holding internal state variables for stability calculations.
 
@@ -84,7 +84,7 @@ A mutable struct holding internal state variables for stability calculations.
   - `q1lim::Float64` - Safety factor derivative at psilim
   - `locstab::Spl.CubicSpline{Float64}` - Spline for local stability analysis
 """
-@kwdef mutable struct DconInternal
+@kwdef mutable struct ForceFreeStatesInternal
     dir_path::String = ""
     mlow::Int = 0
     mhigh::Int = 0
@@ -111,7 +111,7 @@ A mutable struct holding internal state variables for stability calculations.
 end
 
 """
-    DconControl
+    ForceFreeStatesControl
 
 A mutable struct containing control parameters for stability analysis, set by the user in dcon.toml.
 
@@ -171,7 +171,7 @@ A mutable struct containing control parameters for stability analysis, set by th
   - `HDF5_filename::String` - Name of HDF5 output file
   - `force_wv_symmetry::Bool` - Boolean flag to enforce symmetry in the vacuum response matrix
 """
-@kwdef mutable struct DconControl
+@kwdef mutable struct ForceFreeStatesControl
     verbose::Bool = true
     bal_flag::Bool = false
     mat_flag::Bool = false
diff --git a/src/DCON/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl
similarity index 99%
rename from src/DCON/Fourfit.jl
rename to src/ForceFreeStates/Fourfit.jl
index fe6edae6..bdcc765c 100644
--- a/src/DCON/Fourfit.jl
+++ b/src/ForceFreeStates/Fourfit.jl
@@ -142,7 +142,7 @@ function make_metric(equil::Equilibrium.PlasmaEquilibrium; mband::Int, fft_flag:
 end
 
 """
-    make_matrix(metric::MetricData, equil::Equilibrium.PlasmaEquilibrium, intr::DconInternal) -> FourFitVars
+    make_matrix(metric::MetricData, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal) -> FourFitVars
 
 Constructs main DCON matrices for a given toroidal mode number and returns
 them as a new `FourFitVars` object. See the appendix of the 2016 Glasser
@@ -174,7 +174,7 @@ and does not affect the actual matrix sizes, they are all dense.
 Add kinetic metric tensor components for kin_flag = true
 Set powers if necessary
 """
-function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::DconInternal, metric::MetricData)
+function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, metric::MetricData)
 
     # --- Extract inputs ---
     sq = equil.sq
diff --git a/src/DCON/Free.jl b/src/ForceFreeStates/Free.jl
similarity index 91%
rename from src/DCON/Free.jl
rename to src/ForceFreeStates/Free.jl
index 6a2004a4..7f243836 100644
--- a/src/DCON/Free.jl
+++ b/src/ForceFreeStates/Free.jl
@@ -1,5 +1,5 @@
 """
-    free_run!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::DconInternal)
+    free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
 
 Compute the free boundary energies using VACUUM. Performs the same function as `free_run`
 in the Fortran code, except now all data is passed in memory instead of via files. This
@@ -10,7 +10,7 @@ and data dumping.
 ### TODOs
 Check if normalize is ever false, currently always true, and if not, remove related logic
 """
-function free_run!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::DconInternal, wall_settings::Vacuum.WallShapeSettings)
+function free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal, wall_settings::Vacuum.WallShapeSettings)
 
     # TODO: this is always true in fortran - just get rid of it?
     normalize = true
@@ -72,7 +72,7 @@ function free_run!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaE
         # Store block in full wv matrix
         @views vac.wv[(ipert_n-1)*intr.mpert+1 : ipert_n*intr.mpert, (ipert_n-1)*intr.mpert+1 : ipert_n*intr.mpert] .= wv_block
 
-        Vacuum.unset_dcon_params()
+        Vacuum.unset_surface_params()
     end
 
     # Compute complex energy eigenvalues and vectors
@@ -142,7 +142,7 @@ function free_run!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaE
 end
 
 """
-    set_vacuum_inputs(psifac::Float64, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, intr::DconInternal)
+    set_vacuum_inputs(psifac::Float64, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal)
 
 Prepare and write the necessary parameters and boundary shape to VACUUM for computing the vacuum response matrix.
 Performs the same function as `free_write_msc` in the Fortran code, except we will always use in-memory communication.
@@ -152,10 +152,10 @@ Performs the same function as `free_write_msc` in the Fortran code, except we wi
   - `psifac`: Flux surface value at the plasma boundary (Float64)
   - `n`: Toroidal mode number (Int)
   - `equil`: Plasma equilibrium data (Equilibrium.PlasmaEquilibrium)
-  - `intr`: Internal DCON parameters (DconInternal)
+  - `intr`: Internal DCON parameters (ForceFreeStatesInternal)
 
 """
-function set_vacuum_inputs(psifac::Float64, n::Int, equil::Equilibrium.PlasmaEquilibrium, intr::DconInternal, ctrl::DconControl)
+function set_vacuum_inputs(psifac::Float64, n::Int, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl)
 
     # Allocations
     theta_norm = Vector(equil.rzphi.ys)
@@ -185,7 +185,7 @@ function set_vacuum_inputs(psifac::Float64, n::Int, equil::Equilibrium.PlasmaEqu
     end
 
     # Pass all required values to VACUUM
-    Vacuum.set_dcon_params(equil.config.control.mtheta, intr.mlow, intr.mhigh, n, qa,
+    Vacuum.set_surface_params(equil.config.control.mtheta, intr.mlow, intr.mhigh, n, qa,
         reverse(r), reverse(z), reverse(delta))
 
     # For input to the Julia vacuum code
@@ -205,13 +205,13 @@ function set_vacuum_inputs(psifac::Float64, n::Int, equil::Equilibrium.PlasmaEqu
 end
 
 """
-    free_compute_wv_spline(ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, intr::DconInternal)
+    free_compute_wv_spline(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal)
 
 Compute a spline of vacuum response matrices over the range of psi from 'ctrl.psi_edge' to
 `intr.qlim`. This is used for fast evaluation of wt during `ode_record_edge`. Performs the
 same function as `free_wvmats` in the Fortran code.
 """
-function free_compute_wv_spline(ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, intr::DconInternal)
+function free_compute_wv_spline(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal)
 
     # Number of psi grid points for the spline: 4 per q-window minimum
     # TODO: 4 spline points is arbitrary - is there a better way?
@@ -275,7 +275,7 @@ function free_compute_wv_spline(ctrl::DconControl, equil::Equilibrium.PlasmaEqui
             @views wv_array[i, (ipert_n-1)*intr.mpert+1 : ipert_n*intr.mpert, (ipert_n-1)*intr.mpert+1 : ipert_n*intr.mpert] .= wv_block
 
             # Free VACUUM memory
-            Vacuum.unset_dcon_params()
+            Vacuum.unset_surface_params()
         end
     end
 
@@ -283,7 +283,7 @@ function free_compute_wv_spline(ctrl::DconControl, equil::Equilibrium.PlasmaEqui
 end
 
 """
-    free_compute_total(equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::DconInternal, odet::OdeState) -> ComplexF64
+    free_compute_total(equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal, odet::OdeState) -> ComplexF64
 
 Compute total complex energy eigenvalue (total1). This is a trimmed down version of `free_run`
 that only computes the total energy eigenvalue for the mode unstable mode, used in `ode_record_edge_dW`
@@ -291,7 +291,7 @@ which calls this function at each step in the psiedge -> psilim region of integr
 the same function as `free_test` in the Fortran code, except we have moved the creation of the
 wv matrix spline to `free_compute_wv_spline` and pass it in `odet`.wvmat_spline.
 """
-function free_compute_total(equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::DconInternal, odet::OdeState)
+function free_compute_total(equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal, odet::OdeState)
 
     normalize = true
 
diff --git a/src/DCON/Mercier.jl b/src/ForceFreeStates/Mercier.jl
similarity index 100%
rename from src/DCON/Mercier.jl
rename to src/ForceFreeStates/Mercier.jl
diff --git a/src/DCON/Ode.jl b/src/ForceFreeStates/Ode.jl
similarity index 92%
rename from src/DCON/Ode.jl
rename to src/ForceFreeStates/Ode.jl
index c07dad80..0c58ab4f 100644
--- a/src/DCON/Ode.jl
+++ b/src/ForceFreeStates/Ode.jl
@@ -1,5 +1,5 @@
 """
-    `ode_run(ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, intr::DconInternal)`
+    `ode_run(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal)`
 
 Main driver for integrating plasma equilibrium and detecting singular surfaces.
 Has the same functionality as `ode_run` in the Fortran code, with the addition of
@@ -20,7 +20,7 @@ restype functionality if we decide to do this
 
 An OdeState struct containing the final state of the ODE solver after integration is complete.
 """
-function ode_run(ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::DconInternal)
+function ode_run(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
 
     # Initialization
     odet = OdeState(intr.numpert_total, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
@@ -83,7 +83,7 @@ function ode_run(ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::
 end
 
 """
-    ode_axis_init!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, intr::DconInternal)
+    ode_axis_init!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal)
 
 Initialize the OdeState struct for the case of sing_start = 0 (axis initialization). This includes
 determining `psifac`, `psimax`, `ising`, `singfac`, and initializing `u`.
@@ -92,7 +92,7 @@ determining `psifac`, `psimax`, `ising`, `singfac`, and initializing `u`.
 
 Support for `kin_flag`
 """
-function ode_axis_init!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, intr::DconInternal)
+function ode_axis_init!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal)
 
     # Shorthand to evaluate q/q1 inside newton iteration
     qval = psi -> Spl.spline_eval!(equil.sq, psi)[4]
@@ -175,7 +175,7 @@ function ode_sing_init()
 end
 
 """
-    ode_ideal_cross!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::DconInternal)
+    ode_ideal_cross!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
 
 Handle the crossing of a rational surface during ODE integration if `kin_flag` is false.
 Performs the same function as `ode_ideal_cross` in the Fortran code. Differences mainly in integration data
@@ -187,7 +187,7 @@ location and parameters of the next singular surface and writes outputs as desir
 
 Remove while true logic
 """
-function ode_ideal_cross!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::DconInternal)
+function ode_ideal_cross!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
 
     # Fixup solution at singular surface
     ode_unorm!(odet.u, odet, ctrl, intr, true)
@@ -272,7 +272,7 @@ function ode_kin_cross()
 end
 
 """
-    ode_step!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::DconInternal)
+    ode_step!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
 
 Integrate the Euler-Lagrange equations to the next rational surface or edge.
 Performs the same function as `ode_step` in the Fortran code, with the addition of
@@ -290,7 +290,7 @@ the solution at the new point.
 Check sensitivity of results to tolerances, currently using same logic as Fortran
 Check absolute tolerances, currently only relative tolerances are updated
 """
-function ode_step!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::DconInternal)
+function ode_step!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
 
     # Callback to be run at every step, handles fixups, tolerances, and data storage
     cb = DiscreteCallback((u, t, integrator) -> true, integrator_callback!)
@@ -343,7 +343,7 @@ function integrator_callback!(integrator)
 end
 
 """
-    compute_tols(ctrl::DconControl, intr::DconInternal, odet::OdeState)
+    compute_tols(ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, odet::OdeState)
 
 Compute relative and absolute tolerances for the ODE solver based on proximity
 to singular surfaces and magnitude of the solution vectors. In Fortran, this was
@@ -390,7 +390,7 @@ function compute_tols(ctrl, intr, odet)
 end
 
 """
-    ode_unorm!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::DconControl, intr::DconInternal, sing_flag::Bool)
+    ode_unorm!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, sing_flag::Bool)
 
 Computes norms of the solution vectors of the array `u` and normalizes them
 if this is not the first call after a fixup. Throws an error if any vector
@@ -414,7 +414,7 @@ operate on `u` directly without extra copies.
 
 Add resizing logic for unorm arrays when ifix exceeds allocated size
 """
-function ode_unorm!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::DconControl, intr::DconInternal, sing_flag::Bool)
+function ode_unorm!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, sing_flag::Bool)
     
     # Compute norms of first solution vectors, abort if any are zero
     odet.unorm .= norm.(eachcol(u[:, :, 1]))
@@ -445,7 +445,7 @@ function ode_unorm!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::DconControl, i
 end
 
 """
-    ode_fixup!(u::Array{ComplexF64,3}, odet::OdeState, intr::DconInternal, sing_flag::Bool)
+    ode_fixup!(u::Array{ComplexF64,3}, odet::OdeState, intr::ForceFreeStatesInternal, sing_flag::Bool)
 
 Applies Gaussian reduction to orthogonalize solution vectors in `u`. Performs
 the same function as `ode_fixup` in the Fortran code, except now relevant
@@ -456,7 +456,7 @@ description of `ode_unorm!` for more details on the benefits of in-place `u`
 updates.
 
 """
-function ode_fixup!(u::Array{ComplexF64,3}, odet::OdeState, intr::DconInternal, sing_flag::Bool)
+function ode_fixup!(u::Array{ComplexF64,3}, odet::OdeState, intr::ForceFreeStatesInternal, sing_flag::Bool)
 
     # Store data for the current fixup
     ifix = odet.ifix
@@ -492,7 +492,7 @@ function ode_fixup!(u::Array{ComplexF64,3}, odet::OdeState, intr::DconInternal,
 end
 
 """
-    findmax_dW_edge!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::DconInternal)
+    findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
 
 Records the total dW in the integration region between `ctrl.psiedge` and 
 `ctrl.psilim`. This performs the same function as `ode_record_edge` in the
@@ -509,7 +509,7 @@ calculation into `free_compute_wv_spline` and `free_compute_total` respectively
 for clarity. We create the wv matrix spline once prior to the loop.
 
 """
-function findmax_dW_edge!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::DconInternal)
+function findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
 
     # Since we search for the maximum dW, initialize to -Infinity
     fill!(odet.dW_edge, -Inf * (1 + im))
@@ -531,7 +531,7 @@ function findmax_dW_edge!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.
 end
 
 """
-    transform_u!(odet::OdeState, intr::DconInternal)
+    transform_u!(odet::OdeState, intr::ForceFreeStatesInternal)
 
 Constructs the transformation matrices to form the true solution vectors. Effectively
 "undoes" the Gaussian reduction applied during fixups throughout the integration, such that
@@ -540,7 +540,7 @@ Performs a similar function as `idcon_transform` + `idcon_build` in the Fortran
 except we separate the building of the transformation matrices and determining the coefficients
 for a chosen force-free solution, which can be done in postprocessing.
 """
-function transform_u!(odet::OdeState, intr::DconInternal)
+function transform_u!(odet::OdeState, intr::ForceFreeStatesInternal)
 
     # Gaussian reduction matrices for each fixup
     gauss = Array{ComplexF64,3}(undef, intr.numpert_total, intr.numpert_total, odet.ifix)
diff --git a/src/DCON/Sing.jl b/src/ForceFreeStates/Sing.jl
similarity index 93%
rename from src/DCON/Sing.jl
rename to src/ForceFreeStates/Sing.jl
index f55ecb42..c85a44ed 100644
--- a/src/DCON/Sing.jl
+++ b/src/ForceFreeStates/Sing.jl
@@ -1,23 +1,23 @@
 """
-    sing_scan!(intr::DconInternal, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars)
+    sing_scan!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars)
 
 Scan all singular surfaces and calculate asymptotic vmat and mmat matrices
 and Mericer criterion. Performs the same function as `sing_scan` in the Fortran code.
 """
-function sing_scan!(intr::DconInternal, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars)
+function sing_scan!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars)
     for ising in 1:intr.msing
         sing_vmat!(intr, ctrl, equil, ffit, ising)
     end
 end
 
 """
-    sing_find!(intr::DconInternal, equil::Equilibrium.PlasmaEquilibrium)
+    sing_find!(intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEquilibrium)
 
 Locate singular rational q-surfaces (q = m/nn) using a bisection method
 between extrema of the q-profile, and store their properties in `intr.sing`.
 Performs the same function as `sing_find` in the Fortran code.
 """
-function sing_find!(intr::DconInternal, equil::Equilibrium.PlasmaEquilibrium)
+function sing_find!(intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEquilibrium)
 
     # Loop over all toroidal mode numbers
     for n in intr.nlow:intr.nhigh
@@ -73,7 +73,7 @@ function sing_find!(intr::DconInternal, equil::Equilibrium.PlasmaEquilibrium)
 end
 
 """
-    sing_lim!(ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, intr::DconInternal)
+    sing_lim!(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal)
 
 Compute and set integration ψ, q, and q' limits by handling cases where user truncates
 before the last singular surface. Performs a similar function to `sing_lim`
@@ -91,7 +91,7 @@ If `qlim < qmax`, a Newton iteration is performed to find the corresponding
 Note that the Newton iteration will be triggered if either `set_psilim_via_dmlim` is true
 or `ctrl.qhigh < equil.params.qmax`. Otherwise, the equilibrium edge values are used.
 """
-function sing_lim!(intr::DconInternal, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium)
+function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium)
 
     # Initial guesses based on equilibrium
     intr.qlim = min(equil.params.qmax, ctrl.qhigh) # equilibrium solve only goes up to qmax, so we're capped there
@@ -142,7 +142,7 @@ function sing_lim!(intr::DconInternal, ctrl::DconControl, equil::Equilibrium.Pla
 end
 
 """
-    sing_vmat!(intr::DconInternal, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, ising::Int)
+    sing_vmat!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, ising::Int)
 
 Calculate asymptotic vmat and mmat matrices and Mercier criterion for
 singular surface `ising`. Performs the same function as `sing_vmat` in the Fortran code.
@@ -157,7 +157,7 @@ the 2016 Glasser DCON paper for the mathematical details.
 
 Check logic on typing of di
 """
-function sing_vmat!(intr::DconInternal, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, ising::Int)
+function sing_vmat!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, ising::Int)
 
     # Allocations
     singp = intr.sing[ising]
@@ -233,7 +233,7 @@ function sing_vmat!(intr::DconInternal, ctrl::DconControl, equil::Equilibrium.Pl
 end
 
 """
-    sing_mmat!(intr::DconInternal, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, ising::Int)
+    sing_mmat!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, ising::Int)
 
 Calculate asymptotic mmat matrix for singular surface `ising`. Performs the same
 function as `sing_mmat` in the Fortran code. Main differences are 1-indexing for
@@ -262,7 +262,7 @@ Better way to unpack the cubic splines
 Rename variables to be more intuitive? I don't like ff - maybe f and f_fact instead of f_lower
 Add a spline for F directly instead of the lower triangular factorization to avoid complexity?
 """
-function sing_mmat!(intr::DconInternal, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, ising::Int)
+function sing_mmat!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, ising::Int)
 
     # Initial allocations
     singp = intr.sing[ising]
@@ -510,7 +510,7 @@ See equation 47 in the Glass 2016 DCON paper. Identical to the Fortran
   - `singp::SingType`: The singular surface data structure containing all relevant matrices and parameters.
   - `k::Int`: The current order in the power series expansion.
 """
-function sing_solve!(singp::SingType, intr::DconInternal, k::Int)
+function sing_solve!(singp::SingType, intr::ForceFreeStatesInternal, k::Int)
     # TODO: rename this solver_higher_order_vmat?
     # Compute ∑Mₗvₖ₋ₗ
     for l in 1:k
@@ -575,7 +575,7 @@ function sing_matmul(a::Array{ComplexF64,3}, b::Array{ComplexF64,3})
 end
 
 """
-    sing_get_ua(ctrl::DconControl, intr::DconInternal, odet::OdeState)
+    sing_get_ua(ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, odet::OdeState)
 
 Compute the asymptotic series solution for a given singular surface.
 Fills and returns `ua` with the asymptotic solution vmat computed in
@@ -583,7 +583,7 @@ Fills and returns `ua` with the asymptotic solution vmat computed in
 2016 DCON paper. Performs the same function as `sing_get_ua` in the
 Fortran code.
 """
-function sing_get_ua(ctrl::DconControl, intr::DconInternal, odet::OdeState)
+function sing_get_ua(ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, odet::OdeState)
 
     singp = intr.sing[odet.ising]
     r1 = singp.r1
@@ -621,13 +621,13 @@ function sing_get_ua(ctrl::DconControl, intr::DconInternal, odet::OdeState)
 end
 
 """
-    sing_get_ca(ctrl::DconControl, intr::DconInternal, odet::OdeState)
+    sing_get_ca(ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, odet::OdeState)
 
 Compute the asymptotic expansion coefficients according to equation
 50 in Glasser 2016 DCON paper. Performs the same function as
 `sing_get_ca` in the Fortran code.
 """
-function sing_get_ca(ctrl::DconControl, intr::DconInternal, odet::OdeState)
+function sing_get_ca(ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, odet::OdeState)
 
     ua = sing_get_ua(ctrl, intr, odet)
 
@@ -656,7 +656,7 @@ end
     sing_der!(
         du::Array{ComplexF64,3},
         u::Array{ComplexF64,3},
-        params::Tuple{DconControl, Equilibrium.PlasmaEquilibrium, FourFitVars, DconInternal, OdeState},
+        params::Tuple{ForceFreeStatesControl, Equilibrium.PlasmaEquilibrium, FourFitVars, ForceFreeStatesInternal, OdeState},
         psieval::Float64
     )
 
@@ -683,7 +683,7 @@ more simplistic code with similar performance.
 
   - `du::Array{ComplexF64,3}`: Pre-allocated array to hold the derivative result, shape (mpert, mpert, 2), updated in-place
   - `u::Array{ComplexF64,3}`: Current state array, shape (mpert, mpert, 2)
-  - `params::Tuple{DconControl, Equilibrium.PlasmaEquilibrium, FourFitVars, DconInternal, OdeState}`: Tuple of relevant structs
+  - `params::Tuple{ForceFreeStatesControl, Equilibrium.PlasmaEquilibrium, FourFitVars, ForceFreeStatesInternal, OdeState}`: Tuple of relevant structs
   - `psieval::Float64`: Current psi value at which to evaluate the derivative
 
 ### TODOs
@@ -691,8 +691,8 @@ more simplistic code with similar performance.
 Implement kin_flag functionality
 """
 function sing_der!(du::Array{ComplexF64,3}, u::Array{ComplexF64,3},
-    params::Tuple{DconControl,Equilibrium.PlasmaEquilibrium,
-        FourFitVars,DconInternal,OdeState},
+    params::Tuple{ForceFreeStatesControl,Equilibrium.PlasmaEquilibrium,
+        FourFitVars,ForceFreeStatesInternal,OdeState},
     psieval::Float64)
 
     # Unpack structs and initialize
diff --git a/src/DCON/Utils.jl b/src/ForceFreeStates/Utils.jl
similarity index 100%
rename from src/DCON/Utils.jl
rename to src/ForceFreeStates/Utils.jl
diff --git a/src/ForcingTerms/ForcingTerms.jl b/src/ForcingTerms/ForcingTerms.jl
new file mode 100644
index 00000000..3982a025
--- /dev/null
+++ b/src/ForcingTerms/ForcingTerms.jl
@@ -0,0 +1,3 @@
+module ForcingTerms
+# Module for external field specification and forcing term calculations
+end
\ No newline at end of file
diff --git a/src/JPEC.jl b/src/JPEC.jl
index db61868d..0b4259b6 100644
--- a/src/JPEC.jl
+++ b/src/JPEC.jl
@@ -13,9 +13,20 @@ include("Vacuum/Vacuum.jl")
 import .Vacuum as Vacuum
 export Vacuum
 
-include("DCON/DCON.jl")
-import .DCON as DCON
-export DCON
+include("ForceFreeStates/ForceFreeStates.jl")
+import .ForceFreeStates as ForceFreeStates
+export ForceFreeStates
+
+include("ForcingTerms/ForcingTerms.jl")
+import .ForcingTerms as ForcingTerms
+export ForcingTerms
+
+include("PerturbedEquilibrium/PerturbedEquilibrium.jl")
+import .PerturbedEquilibrium as PerturbedEquilibrium
+export PerturbedEquilibrium
+
+include("Main.jl")
+export Main
 
 include(joinpath(@__DIR__, "..", "deps", "build_helpers.jl"))
 export build_fortran, build_spline_fortran, build_vacuum_fortran
diff --git a/src/DCON/Main.jl b/src/Main.jl
similarity index 90%
rename from src/DCON/Main.jl
rename to src/Main.jl
index 36de50a2..c246296d 100644
--- a/src/DCON/Main.jl
+++ b/src/Main.jl
@@ -1,3 +1,16 @@
+using TOML
+using Printf
+using LinearAlgebra
+using HDF5
+
+# Import necessary submodules
+import .ForceFreeStates: ForceFreeStatesInternal, ForceFreeStatesControl, DebugSettings, OdeState, VacuumData, FourFitVars
+import .ForceFreeStates: sing_lim!, sing_find!, sing_scan!, mercier_scan!, make_metric, make_matrix
+import .ForceFreeStates: ode_run, free_run!
+import .Equilibrium
+import .Vacuum
+import .Spl
+
 function Main(path::String="./")
 
     println("DCON START")
@@ -5,10 +18,10 @@ function Main(path::String="./")
     start_time = time()
 
     # Read input data and set up data structures
-    intr = DconInternal(; dir_path=path)
-    # TODO: leaving DCON_CONTROL as a part of the toml file, eventually can combine equil, gpec, etc. into one input file?
+    intr = ForceFreeStatesInternal(; dir_path=path)
+    # TODO: leaving ForceFreeStates_CONTROL as a part of the toml file, eventually can combine equil, gpec, etc. into one input file?
     inputs = TOML.parsefile(joinpath(intr.dir_path, "dcon.toml"))
-    ctrl = DconControl(; (Symbol(k) => v for (k, v) in inputs["DCON_CONTROL"])...)
+    ctrl = ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates_CONTROL"])...)
     equil = Equilibrium.setup_equilibrium(joinpath(intr.dir_path, "equil.toml"))
     if "WALL" in keys(inputs)
         wall_settings = Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["WALL"])...)
@@ -61,7 +74,7 @@ function Main(path::String="./")
 
     # Determine toroidal mode numbers
     if ctrl.nn_low == 0 && ctrl.nn_high == 0
-        error("Either nn_low or nn_high must be set in DCON_CONTROL (both are 0)")
+        error("Either nn_low or nn_high must be set in ForceFreeStates_CONTROL (both are 0)")
     elseif ctrl.nn_low == 0
         ctrl.nn_low = ctrl.nn_high
     elseif ctrl.nn_high == 0
@@ -179,7 +192,7 @@ function Main(path::String="./")
 end
 
 """
-    write_outputs_to_HDF5(ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, intr::DconInternal, odet::OdeState)
+    write_outputs_to_HDF5(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, odet::OdeState)
 
 Helper function to write the HDF5 output file with relevant run and equilibrium parameters.
 This combines the functionality of several pieces of the Fortran code in `ode_output.f`,
@@ -192,13 +205,13 @@ vacuum data if `vac_flag` is true.
 Combine spline unpacking if possible, too many extra lines
 
 """
-function write_outputs_to_HDF5(ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, intr::DconInternal, odet::OdeState, vac::Union{VacuumData, Nothing})
+function write_outputs_to_HDF5(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, odet::OdeState, vac::Union{VacuumData, Nothing})
 
     h5open(joinpath(intr.dir_path, ctrl.HDF5_filename), "w") do out_h5
 
         # Store input parameters
-        for (key, val) in zip(fieldnames(DconControl), getfield.(Ref(ctrl), fieldnames(DconControl)))
-            out_h5["input/DCON_CONTROL/$key"] = val
+        for (key, val) in zip(fieldnames(ForceFreeStatesControl), getfield.(Ref(ctrl), fieldnames(ForceFreeStatesControl)))
+            out_h5["input/ForceFreeStates_CONTROL/$key"] = val
         end
         for (key, val) in zip(fieldnames(Equilibrium.EquilibriumControl), getfield.(Ref(equil.config.control), fieldnames(Equilibrium.EquilibriumControl)))
             out_h5["input/EQUIL_CONTROL/$key"] = val
@@ -206,7 +219,7 @@ function write_outputs_to_HDF5(ctrl::DconControl, equil::Equilibrium.PlasmaEquil
         # TODO: assuming EQUIL_OUTPUT is going to be deprecated
         # TODO: should we store the equilibrium? difficult since it could be a gfile, sol.in, etc.
         # TODO: if we do one input file, can just pass that in instead and loop easily since its parsed
-        # as a dict already (for (k, v) in inputs["DCON_CONTROL"]...). We have to do this since custom structs
+        # as a dict already (for (k, v) in inputs["ForceFreeStates_CONTROL"]...). We have to do this since custom structs
         # don't inherently have an iterator by default
 
         # Write derived run parameters
diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
new file mode 100644
index 00000000..9886d81e
--- /dev/null
+++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
@@ -0,0 +1,3 @@
+module PerturbedEquilibrium
+# Module for perturbed equilibrium calculations
+end
\ No newline at end of file
diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl
index 081821df..79085691 100644
--- a/src/Vacuum/Vacuum.jl
+++ b/src/Vacuum/Vacuum.jl
@@ -5,7 +5,7 @@ using TOML, Interpolations, SpecialFunctions, LinearAlgebra, Printf
 include("VacuumStructs.jl")
 include("VacuumInternals.jl")
 
-export mscvac, set_dcon_params, VacuumInput, compute_vacuum_response
+export mscvac, set_surface_params, VacuumInput, compute_vacuum_response
 export compute_vacuum_field
 export kernel!
 export WallShapeSettings
@@ -18,7 +18,7 @@ const libdir = joinpath(@__DIR__, "..", "..", "deps")
 const libvac = joinpath(libdir, "libvac")
 
 """
-    set_dcon_params(mtheta, lmin, lmax, nnin, qa1in, xin, zin, deltain)
+    set_surface_params(mtheta, lmin, lmax, nnin, qa1in, xin, zin, deltain)
 
 Initialize DCON parameters for vacuum field calculations.
 
@@ -48,14 +48,14 @@ xin = rand(Float64, n_modes)
 zin = rand(Float64, n_modes)
 deltain = rand(Float64, n_modes)
 
-set_dcon_params(mtheta, lmin, lmax, nnin, qa1in, xin, zin, deltain)
+set_surface_params(mtheta, lmin, lmax, nnin, qa1in, xin, zin, deltain)
 ```
 """
-function set_dcon_params(mtheta::Integer, lmin::Integer, lmax::Integer, nnin::Integer,
+function set_surface_params(mtheta::Integer, lmin::Integer, lmax::Integer, nnin::Integer,
     qa1in::Float64,
     xin::Vector{Float64}, zin::Vector{Float64}, deltain::Vector{Float64})
 
-    return ccall((:set_dcon_params_, libvac),
+    return ccall((:set_surface_params_, libvac),
         Nothing,
         (Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint},
             Ref{Cdouble},
@@ -66,30 +66,30 @@ function set_dcon_params(mtheta::Integer, lmin::Integer, lmax::Integer, nnin::In
 end
 
 """
-    unset_dcon_params()
+    unset_surface_params()
 
-Unset DCON parameters previously set by `set_dcon_params`.
+Unset DCON parameters previously set by `set_surface_params`.
 
 This subroutine deallocates in-memory arrays (`x_dcon`, `z_dcon`, and `delta_dcon`)
 and resets the internal DCON state for future vacuum calculations.
 
 # Notes
 
-- Must be called after `set_dcon_params` if you want to reset the DCON memory.
+- Must be called after `set_surface_params` if you want to reset the surface data memory.
 - No arguments are required.
 
 # Example
 
 ```julia
 # Set parameters
-set_dcon_params(mtheta, lmin, lmax, nnin, qa1in, xin, zin, deltain)
+set_surface_params(mtheta, lmin, lmax, nnin, qa1in, xin, zin, deltain)
 
-# Reset DCON parameters
-unset_dcon_params()
+# Reset surface parameters
+unset_surface_params()
 ```
 """
-function unset_dcon_params()
-    ccall((:unset_dcon_params_, libvac), Nothing, ())
+function unset_surface_params()
+    ccall((:unset_surface_params_, libvac), Nothing, ())
 end
 
 """
@@ -109,7 +109,7 @@ Compute the vacuum response matrix for magnetostatic perturbations.
 - `farwall_flag`: Whether to use far-wall approximation (Bool)
 - `grrio`: Green's function data (Array{Float64,2})
 - `xzptso`: Source point coordinates (Array{Float64,2})
-- `op_ahgfile`: Optional communication file for when set_dcon_params is not called (String or Nothing)
+- `op_ahgfile`: Optional communication file for when set_surface_params is not called (String or Nothing)
 
 # Returns
 
@@ -118,13 +118,13 @@ Compute the vacuum response matrix for magnetostatic perturbations.
 
 # Note
 
-Requires prior initialization with `set_dcon_params()` before calling this function.
+Requires prior initialization with `set_surface_params()` before calling this function.
 
 # Examples
 
 ```julia
 # Initialize parameters first
-set_dcon_params(mtheta, lmin, lmax, nnin, qa1in, xin, zin, deltain)
+set_surface_params(mtheta, lmin, lmax, nnin, qa1in, xin, zin, deltain)
 
 # Set up vacuum calculation
 mpert = 5
diff --git a/src/Vacuum/fortran/vacuum_io.f b/src/Vacuum/fortran/vacuum_io.f
index 756ad7fa..ed756e81 100644
--- a/src/Vacuum/fortran/vacuum_io.f
+++ b/src/Vacuum/fortran/vacuum_io.f
@@ -508,14 +508,14 @@ subroutine setahgdir ( directory )
       end
 
 c-----------------------------------------------------------------------
-c     subprogram 7. set_dcon_params.
+c     subprogram 7. set_surface_params.
 c     In memory io for DCON
 c     Speedup over original ascii io (readahg) for multiple calls
 c-----------------------------------------------------------------------
 c-----------------------------------------------------------------------
 c     declarations.
 c-----------------------------------------------------------------------
-      subroutine set_dcon_params ( mthin,lmin,lmax,nnin,qa1in,xin,
+      subroutine set_surface_params ( mthin,lmin,lmax,nnin,qa1in,xin,
      $     zin, deltain )
       USE vglobal_mod, ONLY: r8
       use vglobal_mod, only: dcon_set, mth_dcon, lmin_dcon,lmax_dcon,
@@ -528,7 +528,7 @@ subroutine set_dcon_params ( mthin,lmin,lmax,nnin,qa1in,xin,
 c-----------------------------------------------------------------------
       mthin1 = mthin + 1
       mthin5 = mthin1 + 4
-      if(dcon_set) call unset_dcon_params
+      if(dcon_set) call unset_surface_params
       allocate(x_dcon(mthin5), z_dcon(mthin5), delta_dcon(mthin5))
       mth_dcon = mthin
       lmin_dcon = lmin
@@ -553,13 +553,13 @@ subroutine set_dcon_params ( mthin,lmin,lmax,nnin,qa1in,xin,
       end
 
 c-----------------------------------------------------------------------
-c     subprogram 8. unset_dcon_params.
+c     subprogram 8. unset_surface_params.
 c     Tell vacuum to ignore in-memory dcon params and read ascii values
 c-----------------------------------------------------------------------
 c-----------------------------------------------------------------------
 c     declarations.
 c-----------------------------------------------------------------------
-      subroutine unset_dcon_params
+      subroutine unset_surface_params
       use vglobal_mod, only: dcon_set, x_dcon, z_dcon, delta_dcon
 c-----------------------------------------------------------------------
 c     computations.
diff --git a/test/runtests_fullruns.jl b/test/runtests_fullruns.jl
index feca80fe..68003ba1 100644
--- a/test/runtests_fullruns.jl
+++ b/test/runtests_fullruns.jl
@@ -1,10 +1,10 @@
-# Run DCON.Main on the provided example directories and assert it completes without throwing.
-@testset "Full DCON runs" begin
+# Run ForceFreeStates.Main on the provided example directories and assert it completes without throwing.
+@testset "Full ForceFreeStates runs" begin
     ex1 = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example")
     @info "Running Solovev ideal example"
-    @test isnothing(DCON.Main(ex1))
+    @test isnothing(ForceFreeStates.Main(ex1))
 
     ex2 = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example_multi_n")
     @info "Running Solovev ideal multi-n example"
-    @test isnothing(DCON.Main(ex2))
+    @test isnothing(ForceFreeStates.Main(ex2))
 end
\ No newline at end of file
diff --git a/test/runtests_ode.jl b/test/runtests_ode.jl
index d497f083..258e3158 100644
--- a/test/runtests_ode.jl
+++ b/test/runtests_ode.jl
@@ -1,7 +1,7 @@
 using LinearAlgebra
 
 # TODO: perhaps this isn't the best place for this function?
-# Should I do include("../dcon/utils.jl") instead? or maybe save these functions in a separate file?
+# Should I do include("../ForceFreeStates/utils.jl") instead? or maybe save these functions in a separate file?
 # associated TODO: come up with Gaussian reduction test that doesn't rely on external data
 
 function load_u_matrix(filename)
@@ -30,7 +30,7 @@ end
         # Test that resize_storage! doubles the size of storage arrays
         mpert = 3
         numsteps_init = 10
-        odet = JPEC.DCON.OdeState(mpert, numsteps_init, 10, 5)
+        odet = JPEC.ForceFreeStates.OdeState(mpert, numsteps_init, 10, 5)
         
         # Fill some data
         odet.step = 8
@@ -42,7 +42,7 @@ end
         end
         
         # Resize storage
-        JPEC.DCON.resize_storage!(odet)
+        JPEC.ForceFreeStates.resize_storage!(odet)
         
         # Check new size is doubled
         @test length(odet.psi_store) == 2 * numsteps_init
@@ -59,7 +59,7 @@ end
         end
 
         # Check that you can resize again
-        JPEC.DCON.resize_storage!(odet)
+        JPEC.ForceFreeStates.resize_storage!(odet)
         @test length(odet.psi_store) == 4 * numsteps_init
         @test length(odet.q_store) == 4 * numsteps_init
         @test size(odet.u_store, 4) == 4 * numsteps_init
@@ -70,7 +70,7 @@ end
         # Test that trim_storage! resizes arrays to actual step count
         mpert = 3
         numsteps_init = 20
-        odet = JPEC.DCON.OdeState(mpert, numsteps_init, 10, 5)
+        odet = JPEC.ForceFreeStates.OdeState(mpert, numsteps_init, 10, 5)
         
         # Set step to less than initial size
         odet.step = 12
@@ -82,7 +82,7 @@ end
         end
         
         # Trim storage
-        JPEC.DCON.trim_storage!(odet)
+        JPEC.ForceFreeStates.trim_storage!(odet)
         
         # Check sizes match step count
         @test length(odet.psi_store) == odet.step
@@ -98,70 +98,70 @@ end
     @testset "compute_tols" begin
         # Test tolerance computation
         mpert = 3
-        ctrl = JPEC.DCON.DconControl()
+        ctrl = JPEC.ForceFreeStates.ForceFreeStatesControl()
         ctrl.tol_r = 1e-6
         ctrl.tol_nr = 1e-4
         ctrl.crossover = 0.01
         
-        intr = JPEC.DCON.DconInternal(; mpert=mpert)
+        intr = JPEC.ForceFreeStates.ForceFreeStatesInternal(; mpert=mpert)
         intr.msing = 2
-        intr.sing = [JPEC.DCON.SingType(), JPEC.DCON.SingType()]
+        intr.sing = [JPEC.ForceFreeStates.SingType(), JPEC.ForceFreeStates.SingType()]
         intr.sing[1].q = 2.0
         intr.sing[1].n = [1]
         intr.sing[2].q = 3.0
         intr.sing[2].n = [1]
         
-        odet = JPEC.DCON.OdeState(mpert, 10, 10, 2)
+        odet = JPEC.ForceFreeStates.OdeState(mpert, 10, 10, 2)
         
         # Test 1: Far from singular surface (singfac > crossover)
         odet.ising = 1
         odet.q = 1.5  # Far from q=2.0
-        rtol = JPEC.DCON.compute_tols(ctrl, intr, odet)
+        rtol = JPEC.ForceFreeStates.compute_tols(ctrl, intr, odet)
         @test rtol == ctrl.tol_nr  # Should use non-resonant tolerance
         
         # Test 2: Close to singular surface (singfac < crossover)
         odet.ising = 1
         odet.q = 1.999  # Very close to q=2.0
-        rtol = JPEC.DCON.compute_tols(ctrl, intr, odet)
+        rtol = JPEC.ForceFreeStates.compute_tols(ctrl, intr, odet)
         @test rtol == ctrl.tol_r  # Should use resonant tolerance
         
         # Test 3: Between two singular surfaces
         odet.ising = 2
         odet.q = 2.5  # Between q=2.0 and q=3.0
-        rtol = JPEC.DCON.compute_tols(ctrl, intr, odet)
+        rtol = JPEC.ForceFreeStates.compute_tols(ctrl, intr, odet)
         @test rtol == ctrl.tol_nr  # Should use min distance to either surface
         
         # Test 4: Beyond all singular surfaces
         odet.ising = 3
         odet.q = 4.0
-        rtol = JPEC.DCON.compute_tols(ctrl, intr, odet)
+        rtol = JPEC.ForceFreeStates.compute_tols(ctrl, intr, odet)
         @test rtol == ctrl.tol_nr
 
         # Edge case - no singular surfaces
         mpert = 2
-        ctrl = JPEC.DCON.DconControl()
+        ctrl = JPEC.ForceFreeStates.ForceFreeStatesControl()
         ctrl.tol_r = 1e-6
         ctrl.tol_nr = 1e-4
         ctrl.crossover = 0.01
         
-        intr = JPEC.DCON.DconInternal(; mpert=mpert)
+        intr = JPEC.ForceFreeStates.ForceFreeStatesInternal(; mpert=mpert)
         intr.msing = 0
         intr.sing = []
         
-        odet = JPEC.DCON.OdeState(mpert, 10, 10, 0)
+        odet = JPEC.ForceFreeStates.OdeState(mpert, 10, 10, 0)
         odet.ising = 1
         odet.q = 2.0
         
         # Should return non-resonant tolerance when no singular surfaces
-        rtol = JPEC.DCON.compute_tols(ctrl, intr, odet)
+        rtol = JPEC.ForceFreeStates.compute_tols(ctrl, intr, odet)
         @test rtol == ctrl.tol_nr
     end
 
     @testset "transform_u!" begin
         # Test transformation of solution vectors
         mpert = 2
-        intr = JPEC.DCON.DconInternal(; mpert=mpert, numpert_total=mpert)
-        odet = JPEC.DCON.OdeState(mpert, 10, 5, 2)
+        intr = JPEC.ForceFreeStates.ForceFreeStatesInternal(; mpert=mpert, numpert_total=mpert)
+        odet = JPEC.ForceFreeStates.OdeState(mpert, 10, 5, 2)
         
         # Set up a simple fixup scenario
         odet.ifix = 1
@@ -190,7 +190,7 @@ end
         u_orig = copy(odet.u_store)
         
         # Apply transformation
-        JPEC.DCON.transform_u!(odet, intr)
+        JPEC.ForceFreeStates.transform_u!(odet, intr)
         
         # Check that u_store was modified (transformation applied)
         @test !all(odet.u_store .== u_orig)
@@ -205,16 +205,16 @@ end
         # Initialize to random u
         mpert = 5
         ifix = 1
-        odet = JPEC.DCON.OdeState(mpert, 10, 10, 10)
+        odet = JPEC.ForceFreeStates.OdeState(mpert, 10, 10, 10)
         odet.u = randn(ComplexF64, mpert, mpert, 2)
         odet.unorm = [norm(odet.u[:, i, 1]) for i in 1:mpert]
         odet.ifix = ifix
         odet.fixfac = zeros(ComplexF64, mpert, mpert, ifix)
-        intr = JPEC.DCON.DconInternal(; numpert_total=mpert)
+        intr = JPEC.ForceFreeStates.ForceFreeStatesInternal(; numpert_total=mpert)
 
         # Save copy of original u and run
         u_orig = copy(odet.u)
-        JPEC.DCON.ode_fixup!(odet.u, odet, intr, false)
+        JPEC.ForceFreeStates.ode_fixup!(odet.u, odet, intr, false)
 
         # Very simple tests
         @test !all(odet.u .== u_orig)  # u should have changed
@@ -224,7 +224,7 @@ end
 
         # --- Real Fortran data check ---
         mpert = 31
-        odet = JPEC.DCON.OdeState(mpert, 10, 10, 10)
+        odet = JPEC.ForceFreeStates.OdeState(mpert, 10, 10, 10)
         # We'll load in Fortran data for u pulled before and after a fixup
         # Note that this was generated by manually setting
         # unorm = [norm(odet.u[:,i,1]) for i in 1:msol] in the Fortran to avoid
@@ -233,9 +233,9 @@ end
         odet.unorm = [norm(odet.u[:, i, 1]) for i in 1:mpert]
         odet.ifix = ifix
         odet.fixfac = zeros(ComplexF64, mpert, mpert, ifix)
-        intr = JPEC.DCON.DconInternal(; numpert_total=mpert)
+        intr = JPEC.ForceFreeStates.ForceFreeStatesInternal(; numpert_total=mpert)
 
-        JPEC.DCON.ode_fixup!(odet.u, odet, intr, false)
+        JPEC.ForceFreeStates.ode_fixup!(odet.u, odet, intr, false)
 
         u_fortran = load_u_matrix(joinpath(@__DIR__, "test_data", "u_postfixup.dat"))
         # test that the outputs are approximately equivalent (1e-3 seems ok to account for loading differences)
@@ -243,7 +243,7 @@ end
 
         # Test with a simple 2x2 case where we can predict the result
         mpert = 2
-        odet = JPEC.DCON.OdeState(mpert, 10, 10, 10)
+        odet = JPEC.ForceFreeStates.OdeState(mpert, 10, 10, 10)
         
         # Set up a simple u matrix where first column has larger norm
         # u[:, 1, 1] = [3, 4] (norm = 5)
@@ -255,11 +255,11 @@ end
         odet.unorm = [norm(odet.u[:, i, 1]) for i in 1:mpert]
         odet.ifix = 1
         odet.fixfac = zeros(ComplexF64, mpert, mpert, 1)
-        intr = JPEC.DCON.DconInternal(; numpert_total=mpert)
+        intr = JPEC.ForceFreeStates.ForceFreeStatesInternal(; numpert_total=mpert)
         
         u_before = copy(odet.u)
         
-        JPEC.DCON.ode_fixup!(odet.u, odet, intr, false)
+        JPEC.ForceFreeStates.ode_fixup!(odet.u, odet, intr, false)
         
         # After fixup:
         # - index should sort by norm: [1, 2] (largest first)
@@ -277,9 +277,9 @@ end
 
     @testset "ode_unorm!" begin
         mpert = 2
-        odet = JPEC.DCON.OdeState(mpert, 10, 10, 10)
-        intr = JPEC.DCON.DconInternal(; mpert=mpert)
-        ctrl = JPEC.DCON.DconControl()
+        odet = JPEC.ForceFreeStates.OdeState(mpert, 10, 10, 10)
+        intr = JPEC.ForceFreeStates.ForceFreeStatesInternal(; mpert=mpert)
+        ctrl = JPEC.ForceFreeStates.ForceFreeStatesControl()
         ctrl.ucrit = 10.0
 
         # Case 1: Basic norm computation
@@ -287,7 +287,7 @@ end
         odet.u[:, 1, 1] .= [3, 4]          # norm = 5
         odet.u[:, 2, 1] .= [0, 2]          # norm = 2
 
-        JPEC.DCON.ode_unorm!(odet.u, odet, ctrl, intr, false)
+        JPEC.ForceFreeStates.ode_unorm!(odet.u, odet, ctrl, intr, false)
         # After the first run with new=True (default), unorm0 should be set to unorm
         # and new should be false
         @test odet.unorm[1:intr.mpert] ≈ [5, 2]
@@ -297,27 +297,27 @@ end
         # Case 2: Error on zero norm
         odet.u[:, 1, 1] .= 0
         odet.new = true
-        @test_throws ErrorException JPEC.DCON.ode_unorm!(odet.u, odet, ctrl, intr, false)
+        @test_throws ErrorException JPEC.ForceFreeStates.ode_unorm!(odet.u, odet, ctrl, intr, false)
 
         # Case 3: Normalization on second call
         odet.u[:, 1, 1] .= [3, 4]   # norm = 5
         odet.u[:, 2, 1] .= [0, 2]   # norm = 2
         odet.new = false
-        JPEC.DCON.ode_unorm!(odet.u, odet, ctrl, intr, false)
+        JPEC.ForceFreeStates.ode_unorm!(odet.u, odet, ctrl, intr, false)
         @test odet.unorm[1:intr.mpert] ≈ [1, 1]
 
         # Case 4: Trigger fixup via ucrit
         odet.unorm0 = ones(intr.mpert)
         odet.u[:, 1, 1] .= [1000, 0]   # large norm
         odet.u[:, 2, 1] .= [1, 0]      # small norm
-        JPEC.DCON.ode_unorm!(odet.u, odet, ctrl, intr, false)
+        JPEC.ForceFreeStates.ode_unorm!(odet.u, odet, ctrl, intr, false)
         @test odet.new == true  # implies fixup ran
 
         # Case 5: Trigger fixup via sing_flag
         odet.new = false
         odet.u[:, 1, 1] .= [1, 0]
         odet.u[:, 2, 1] .= [1, 0]
-        JPEC.DCON.ode_unorm!(odet.u, odet, ctrl, intr, true)
+        JPEC.ForceFreeStates.ode_unorm!(odet.u, odet, ctrl, intr, true)
         @test odet.new == true  # fixup triggered
     end
 
@@ -328,7 +328,7 @@ end
         numunorms_init = 20
         msing = 10
         
-        odet = JPEC.DCON.OdeState(numpert_total, numsteps_init, numunorms_init, msing)
+        odet = JPEC.ForceFreeStates.OdeState(numpert_total, numsteps_init, numunorms_init, msing)
         
         # Check fields are initialized correctly
         @test odet.numpert_total == numpert_total
diff --git a/test/runtests_vacuum_fortran.jl b/test/runtests_vacuum_fortran.jl
index f3299bd0..6d6cda9b 100644
--- a/test/runtests_vacuum_fortran.jl
+++ b/test/runtests_vacuum_fortran.jl
@@ -6,9 +6,9 @@
         zin = rand(Float64, lmax - lmin + 1)
         deltain = rand(Float64, lmax - lmin + 1)
 
-        @info "Testing set_dcon_params…"
-        JPEC.Vacuum.set_dcon_params(mtheta_in, lmin, lmax, nnin, qa1in, xin, zin, deltain)
-        @info "set_dcon_params OK!"
+        @info "Testing set_surface_params…"
+        JPEC.Vacuum.set_surface_params(mtheta_in, lmin, lmax, nnin, qa1in, xin, zin, deltain)
+        @info "set_surface_params OK!"
 
         @info "Testing mscvac…"
         mpert = 5
diff --git a/test/test_data/regression_solovev_ideal_example/dcon.toml b/test/test_data/regression_solovev_ideal_example/ForceFreeStates.toml
similarity index 99%
rename from test/test_data/regression_solovev_ideal_example/dcon.toml
rename to test/test_data/regression_solovev_ideal_example/ForceFreeStates.toml
index 118a6aff..0a45b560 100644
--- a/test/test_data/regression_solovev_ideal_example/dcon.toml
+++ b/test/test_data/regression_solovev_ideal_example/ForceFreeStates.toml
@@ -1,4 +1,4 @@
-[DCON_CONTROL]
+[ForceFreeStates_CONTROL]
 bal_flag = false              # Ideal MHD ballooning criterion for short wavelengths
 mat_flag = true               # Construct coefficient matrices for diagnostic purposes
 ode_flag = true               # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
diff --git a/test/test_data/regression_solovev_ideal_example_multi_n/dcon.toml b/test/test_data/regression_solovev_ideal_example_multi_n/ForceFreeStates.toml
similarity index 99%
rename from test/test_data/regression_solovev_ideal_example_multi_n/dcon.toml
rename to test/test_data/regression_solovev_ideal_example_multi_n/ForceFreeStates.toml
index 3396fc58..367770d8 100644
--- a/test/test_data/regression_solovev_ideal_example_multi_n/dcon.toml
+++ b/test/test_data/regression_solovev_ideal_example_multi_n/ForceFreeStates.toml
@@ -1,4 +1,4 @@
-[DCON_CONTROL]
+[ForceFreeStates_CONTROL]
 bal_flag = false              # Ideal MHD ballooning criterion for short wavelengths
 mat_flag = true               # Construct coefficient matrices for diagnostic purposes
 ode_flag = true               # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)