From f589bd7f120af2033c681d5044717a1320d6cb5b Mon Sep 17 00:00:00 2001
From: james <81617086+je-cook@users.noreply.github.com>
Date: Mon, 8 Jun 2026 10:08:55 +0100
Subject: [PATCH 1/7] scan cleanup

---
 process/core/scan.py | 170 +++++++++++++++----------------------------
 1 file changed, 57 insertions(+), 113 deletions(-)

diff --git a/process/core/scan.py b/process/core/scan.py
index d62a84a769..6956b3e59d 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -18,7 +18,7 @@
 from process.core.log import logging_model_handler, show_errors
 from process.core.solver import constraints
 from process.core.solver.solver_handler import SolverHandler
-from process.data_structure.scan_variables import IPNSCNS, NOUTVARS, ScanData
+from process.data_structure.scan_variables import IPNSCNS, NOUTVARS
 from process.models.availability import AvailabilityModel
 
 if TYPE_CHECKING:
@@ -325,20 +325,18 @@ def scan_1d(self):
             max_sweep_value_length - len(str(sweep_val).replace(".", ""))
             for sweep_val in sweep_values
         ]
-        for iscan in range(1, self.data.scan.isweep + 1):
-            if scan_1d_ifail_dict[iscan] == 1:
+        for iscan in range(self.data.scan.isweep):
+            pstring = (
+                f"Scan {iscan:02d}: {nsweep_var.name} = {sweep_values[iscan]} "
+                + " " * offsets[iscan]
+                + "\u001b[32m{}CONVERGED \u001b[0m"
+            )
+            if scan_1d_ifail_dict[iscan + 1] == 1:
                 converged_count += 1
-                print(
-                    f"Scan {iscan:02d}: {nsweep_var.fname} = {sweep_values[iscan - 1]} "
-                    + " " * offsets[iscan - 1]
-                    + "\u001b[32mCONVERGED \u001b[0m"
-                )
+                pstring.format("")
             else:
-                print(
-                    f"Scan {iscan:02d}: {nsweep_var.fname} = {sweep_values[iscan - 1]} "
-                    + " " * offsets[iscan - 1]
-                    + "\u001b[31mUNCONVERGED \u001b[0m"
-                )
+                pstring.format("UN")
+            print(pstring)
         converged_percentage = converged_count / self.data.scan.isweep * 100
         print(f"\nConvergence Percentage: {converged_percentage:.2f}%")
 
@@ -399,95 +397,59 @@ def scan_2d(self):
 
         for iscan_1 in range(1, self.data.scan.isweep + 1):
             for iscan_2 in range(1, self.data.scan.isweep_2 + 1):
+                string = (
+                    f"Scan {scan_point:02d}: ({nsweep_var.name} = "
+                    f"{sweep_1_values[iscan_1 - 1]}, {nsweep_2_var.name} "
+                    f"= {sweep_2_values[iscan_2 - 1]}) "
+                    + " " * offsets[iscan_1 - 1][iscan_2 - 1]
+                    + "\u001b[32m{}CONVERGED \u001b[0m"
+                )
                 if scan_2d_ifail_list[iscan_1][iscan_2] == 1:
                     converged_count += 1
-                    print(
-                        (
-                            f"Scan {scan_point:02d}: ({nsweep_var.fname} = "
-                            f"{sweep_1_values[iscan_1 - 1]}, {nsweep_2_var.fname} "
-                            f"= {sweep_2_values[iscan_2 - 1]}) "
-                        )
-                        + " " * offsets[iscan_1 - 1][iscan_2 - 1]
-                        + "\u001b[32mCONVERGED \u001b[0m"
-                    )
-                    scan_point += 1
+                    print(string.format())
                 else:
-                    print(
-                        (
-                            f"Scan {scan_point:02d}: ({nsweep_var.fname} = "
-                            f"{sweep_1_values[iscan_1 - 1]}, {nsweep_2_var.fname} = "
-                            f"{sweep_2_values[iscan_2 - 1]}) "
-                        )
-                        + " " * offsets[iscan_1 - 1][iscan_2 - 1]
-                        + "\u001b[31mUNCONVERGED \u001b[0m"
-                    )
-                    scan_point += 1
+                    print(string.format("UN"))
+                scan_point += 1
         converged_percentage = (
             converged_count / (self.data.scan.isweep * self.data.scan.isweep_2) * 100
         )
         print(f"\nConvergence Percentage: {converged_percentage:.2f}%")
 
-    @staticmethod
-    def scan_2d_init(scan_data: ScanData):
-        """Scan 2d initialisation"""
-        process_output.ovarre(
-            constants.MFILE,
-            "Number of first variable scan points",
-            "(isweep)",
-            scan_data.isweep,
-        )
-        process_output.ovarre(
-            constants.MFILE,
-            "Number of second variable scan points",
-            "(isweep_2)",
-            scan_data.isweep_2,
-        )
-        process_output.ovarre(
-            constants.MFILE,
-            "Scanning first variable number",
-            "(nsweep)",
-            scan_data.nsweep,
-        )
-        process_output.ovarre(
-            constants.MFILE,
-            "Scanning second variable number",
-            "(nsweep_2)",
-            scan_data.nsweep_2,
-        )
-        process_output.ovarre(
-            constants.MFILE,
-            "Scanning second variable number",
-            "(nsweep_2)",
-            scan_data.nsweep_2,
-        )
-        process_output.ovarre(
-            constants.MFILE,
-            "Scanning second variable number",
-            "(nsweep_2)",
-            scan_data.nsweep_2,
-        )
+    def scan_2d_init(self):
+        sv = self.data.scan
+        for d, n, v in (
+            ("Number of first variable scan points", "(isweep)", sv.isweep),
+            ("Number of second variable scan points", "(isweep_2)", sv.isweep_2),
+            ("Scanning first variable number", "(nsweep)", sv.nsweep),
+            ("Scanning second variable number", "(nsweep_2)", sv.nsweep_2),
+            ("Scanning second variable number", "(nsweep_2)", sv.nsweep_2),
+            ("Scanning second variable number", "(nsweep_2)", sv.nsweep_2),
+        ):
+            process_output.ovarin(constants.MFILE, d, n, v)
+
+    def _set_v_x_label(self, iscan, twod=False):
+        if twod:
+            sv = self.scan_select(self.data.scan.nsweep_2, self.data.scan.sweep_2, iscan)
+        else:
+            sv = self.scan_select(self.data.scan.nsweep, self.data.scan.sweep, iscan)
+        self.data.globals.vlabel = sv.fname
+        self.data.globals.xlabel = sv.description
 
     def scan_1d_write_point_header(self, iscan: int):
         """Scan 1d header"""
         self.data.globals.iscan_global = iscan
-        sv = self.scan_select(self.data.scan.nsweep, self.data.scan.sweep, iscan)
-
-        self.data.globals.vlabel = sv.fname
-        self.data.globals.xlabel = sv.description
+        self._set_v_x_label(iscan)
 
         process_output.oblnkl(constants.NOUT)
-        process_output.ostars(constants.NOUT, 110)
+        process_output.oblnkl(constants.MFILE)
 
         process_output.write(
             constants.NOUT,
-            f"***** Scan point {iscan} of {self.data.scan.isweep} : "
-            f"{self.data.globals.xlabel}"
+            f"Scan point {iscan} of {self.data.scan.isweep} : {self.data.globals.xlabel}"
             f", {self.data.globals.vlabel} = {self.data.scan.sweep[iscan - 1]} "
-            "*****",
+            "",
         )
-        process_output.ostars(constants.NOUT, 110)
-        process_output.oblnkl(constants.MFILE)
-        process_output.ovarre(constants.MFILE, "Scan point number", "(iscan)", iscan)
+        process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan)
 
         print(
             f"Starting scan point {iscan} of {self.data.scan.isweep} : "
@@ -501,18 +463,12 @@ def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2):
 
         # Makes iscan available globally (read-only)
         self.data.globals.iscan_global = iscan
-        sv_1 = self.scan_select(self.data.scan.nsweep, self.data.scan.sweep, iscan_1)
-
-        self.data.globals.vlabel = sv_1.fname
-        self.data.globals.xlabel = sv_1.data.description
 
-        sv_2 = self.scan_select(self.data.scan.nsweep_2, self.data.scan.sweep_2, iscan_r)
-
-        self.data.globals.vlabel_2 = sv_2.fname
-        self.data.globals.xlabel_2 = sv_2.data.description
+        self._set_v_x_label(iscan_1)
+        self._set_v_x_label(iscan_r)
 
         process_output.oblnkl(constants.NOUT)
-        process_output.ostars(constants.NOUT, 110)
+        process_output.oblnkl(constants.MFILE)
 
         process_output.write(
             constants.NOUT,
@@ -522,9 +478,7 @@ def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2):
             f" {self.data.globals.vlabel_2} = {self.data.scan.sweep_2[iscan_r - 1]} "
             "*****",
         )
-        process_output.ostars(constants.NOUT, 110)
-        process_output.oblnkl(constants.MFILE)
-        process_output.ovarre(constants.MFILE, "Scan point number", "(iscan)", iscan)
+        process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan)
 
         print(
             f"Starting scan point {iscan}:  {self.data.globals.xlabel}, "
@@ -535,27 +489,17 @@ def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2):
 
         return iscan_r
 
-    @staticmethod
-    def scan_1d_write_plot(scan_data: ScanData):
-        """Scan 1d plotter"""
-        if scan_data.first_call_1d:
-            process_output.ovarre(
-                constants.MFILE,
-                "Number of scan points",
-                "(isweep)",
-                scan_data.isweep,
-            )
-            process_output.ovarre(
-                constants.MFILE,
-                "Scanning variable number",
-                "(nsweep)",
-                scan_data.nsweep,
-            )
+    def scan_1d_write_plot(self):
+        if self.data.scan.first_call_1d:
+            for d, n, v in (
+                ("Number of scan points", "(isweep)", self.data.scan.isweep),
+                ("Scanning variable number", "(nsweep)", self.data.scan.nsweep),
+            ):
+                process_output.ovarin(constants.MFILE, d, n, v)
 
-            scan_data.first_call_1d = False
+            self.data.scan.first_call_1d = False
 
     def scan_select(self, nsweep, sweep, iscan) -> ScanVariables:
         """Select a scan"""
         sv = ScanVariables(nsweep)
         sv.set(self.data, sweep[iscan - 1])
-        return sv

From 1ba1c052dc38504cbf7ed37a7d3a37f56571685b Mon Sep 17 00:00:00 2001
From: james <81617086+je-cook@users.noreply.github.com>
Date: Sat, 13 Jun 2026 10:13:21 +0100
Subject: [PATCH 2/7] scan variable cleanup

---
 process/core/scan.py | 57 ++++++++++++++++++++++++--------------------
 1 file changed, 31 insertions(+), 26 deletions(-)

diff --git a/process/core/scan.py b/process/core/scan.py
index 6956b3e59d..5876ad0f0e 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -329,13 +329,13 @@ def scan_1d(self):
             pstring = (
                 f"Scan {iscan:02d}: {nsweep_var.name} = {sweep_values[iscan]} "
                 + " " * offsets[iscan]
-                + "\u001b[32m{}CONVERGED \u001b[0m"
+                + "\u001b[3{}CONVERGED \u001b[0m"
             )
             if scan_1d_ifail_dict[iscan + 1] == 1:
                 converged_count += 1
-                pstring.format("")
+                pstring.format("2m")
             else:
-                pstring.format("UN")
+                pstring.format("1mUN")
             print(pstring)
         converged_percentage = converged_count / self.data.scan.isweep * 100
         print(f"\nConvergence Percentage: {converged_percentage:.2f}%")
@@ -369,7 +369,24 @@ def scan_2d(self):
                 scan_2d_ifail_list[iscan_1][iscan_2] = ifail
                 iscan += 1
 
-        print("Scan Convergence Summary\n")
+        self.output_2d_summary(scan_2d_ifail_list)
+
+    def scan_2d_init(self):
+        sv = self.data.scan
+        for d, n, v in (
+            ("Number of first variable scan points", "(isweep)", sv.isweep),
+            ("Number of second variable scan points", "(isweep_2)", sv.isweep_2),
+            ("Scanning first variable number", "(nsweep)", sv.nsweep),
+            ("Scanning second variable number", "(nsweep_2)", sv.nsweep_2),
+            ("Scanning second variable number", "(nsweep_2)", sv.nsweep_2),
+            ("Scanning second variable number", "(nsweep_2)", sv.nsweep_2),
+        ):
+            process_output.ovarin(constants.MFILE, d, n, v)
+
+    def output_2d_summary(self, scan_2d_ifail_list):
+        print(
+            " ****************************************** Scan Convergence Summary ****************************************** \n"
+        )
         sweep_1_values = self.data.scan.sweep[: self.data.scan.isweep]
         sweep_2_values = self.data.scan.sweep_2[: self.data.scan.isweep_2]
         nsweep_var = self.scan_select(
@@ -402,36 +419,25 @@ def scan_2d(self):
                     f"{sweep_1_values[iscan_1 - 1]}, {nsweep_2_var.name} "
                     f"= {sweep_2_values[iscan_2 - 1]}) "
                     + " " * offsets[iscan_1 - 1][iscan_2 - 1]
-                    + "\u001b[32m{}CONVERGED \u001b[0m"
+                    + "\u001b[3{}CONVERGED \u001b[0m"
                 )
                 if scan_2d_ifail_list[iscan_1][iscan_2] == 1:
                     converged_count += 1
-                    print(string.format())
+                    print(string.format("2m"))
                 else:
-                    print(string.format("UN"))
+                    print(string.format("1mUN"))
                 scan_point += 1
         converged_percentage = (
             converged_count / (self.data.scan.isweep * self.data.scan.isweep_2) * 100
         )
         print(f"\nConvergence Percentage: {converged_percentage:.2f}%")
 
-    def scan_2d_init(self):
-        sv = self.data.scan
-        for d, n, v in (
-            ("Number of first variable scan points", "(isweep)", sv.isweep),
-            ("Number of second variable scan points", "(isweep_2)", sv.isweep_2),
-            ("Scanning first variable number", "(nsweep)", sv.nsweep),
-            ("Scanning second variable number", "(nsweep_2)", sv.nsweep_2),
-            ("Scanning second variable number", "(nsweep_2)", sv.nsweep_2),
-            ("Scanning second variable number", "(nsweep_2)", sv.nsweep_2),
-        ):
-            process_output.ovarin(constants.MFILE, d, n, v)
-
     def _set_v_x_label(self, iscan, twod=False):
-        if twod:
-            sv = self.scan_select(self.data.scan.nsweep_2, self.data.scan.sweep_2, iscan)
-        else:
-            sv = self.scan_select(self.data.scan.nsweep, self.data.scan.sweep, iscan)
+        sv = (
+            self.scan_select(self.data.scan.nsweep_2, self.data.scan.sweep_2, iscan)
+            if twod
+            else self.scan_select(self.data.scan.nsweep, self.data.scan.sweep, iscan)
+        )
         self.data.globals.vlabel = sv.fname
         self.data.globals.xlabel = sv.description
 
@@ -472,11 +478,10 @@ def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2):
 
         process_output.write(
             constants.NOUT,
-            f"***** 2D Scan point {iscan} of "
-            f"{self.data.scan.isweep * self.data.scan.isweep_2} : "
+            f"2D Scan point {iscan} of {self.data.scan.isweep * self.data.scan.isweep_2} : "
             f"{self.data.globals.vlabel} = {self.data.scan.sweep[iscan_1 - 1]} and"
             f" {self.data.globals.vlabel_2} = {self.data.scan.sweep_2[iscan_r - 1]} "
-            "*****",
+            "",
         )
         process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan)
 

From 2f1d211b4f9d44743a5a355e95ee90fa3bc134dc Mon Sep 17 00:00:00 2001
From: james <81617086+je-cook@users.noreply.github.com>
Date: Sat, 13 Jun 2026 10:14:07 +0100
Subject: [PATCH 3/7] solver_handler cleanup

---
 process/core/solver/solver_handler.py | 55 +++++++++++++--------------
 1 file changed, 26 insertions(+), 29 deletions(-)

diff --git a/process/core/solver/solver_handler.py b/process/core/solver/solver_handler.py
index 4901d5983a..6c4cb30c87 100644
--- a/process/core/solver/solver_handler.py
+++ b/process/core/solver/solver_handler.py
@@ -1,6 +1,7 @@
 """Module containing solver handler routines"""
 
 import logging
+from contextlib import contextmanager
 
 from tabulate import tabulate
 
@@ -50,14 +51,9 @@ def run(self):
         # Initialise iteration variables and bounds in Python: relies on Fortran
         # iteration variables being defined above
         # Trim maximum size arrays down to actually used size
-        n = self.data.numerics.nvar
-        x = self.data.numerics.xcm[:n]
-        bndl = self.data.numerics.itv_scaled_lower_bounds[:n]
-        bndu = self.data.numerics.itv_scaled_upper_bounds[:n]
-
-        # Define total number of constraints and equality constraints
-        m = self.data.numerics.neqns + self.data.numerics.nineqns
-        meq = self.data.numerics.neqns
+        x = self.data.numerics.xcm[: self.data.numerics.nvar]
+        bndl = self.data.numerics.itv_scaled_lower_bounds[: self.data.numerics.nvar]
+        bndu = self.data.numerics.itv_scaled_upper_bounds[: self.data.numerics.nvar]
 
         # Evaluators() calculates the objective and constraint functions and
         # their gradients for a given vector x
@@ -68,31 +64,18 @@ def run(self):
         self.solver.set_evaluators(evaluators)
         self.solver.set_bounds(bndl, bndu)
         self.solver.set_opt_params(x)
-        self.solver.set_constraints(m, meq)
+        # Define total number of constraints and equality constraints
+        self.solver.set_constraints(
+            m=self.data.numerics.neqns + self.data.numerics.nineqns,
+            meq=self.data.numerics.neqns,
+        )
         ifail = self.solver.solve()
 
         # If VMCON optimisation has failed then try altering value of epsfcn
         if self.solver_name == "vmcon":
-            if ifail != SolverOutputCondition.CONVERGED:
-                print("Trying again with new epsfcn")
-                # epsfcn is only used in evaluators.Evaluators()
-                # TODO epsfcn could be set in Evaluators instance now, don't need to
-                # set/unset in self.data.numerics module
-                self.data.numerics.epsfcn *= 10  # try new larger value
-                print("new epsfcn = ", self.data.numerics.epsfcn)
-
-                ifail = self.solver.solve()
-                # First solution attempt failed
-                # (ifail != SolverOutputCondition.CONVERGED): supply ifail value
-                # to next attempt
-                self.data.numerics.epsfcn /= 10  # reset value
-
-            if ifail != SolverOutputCondition.CONVERGED:
-                print("Trying again with new epsfcn")
-                self.data.numerics.epsfcn /= 10  # try new smaller value
-                print("new epsfcn = ", self.data.numerics.epsfcn)
-                ifail = self.solver.solve()
-                self.data.numerics.epsfcn *= 10  # reset value
+            if ifail != 1:
+                with epsfcn_context(self.data.numerics):
+                    self.solver.solve()
 
             # If VMCON has exited with error code 5
             # (ifail = SolverOutputCondition.NO_SOLUTION) try another run using a
@@ -346,3 +329,17 @@ def _optimisation_parameters_output(self):
                 numalign="left",
             ),
         )
+
+
+@contextmanager
+def epsfcn_context(numerics):
+    print("Trying again with new epsfcn")
+    # epsfcn is only used in evaluators.Evaluators()
+    # TODO epsfcn could be set in Evaluators instance now, don't need to
+    # set/unset in numerics module
+    numerics.epsfcn *= 10  # try new larger value
+    print("new epsfcn = ", numerics.epsfcn)
+    try:
+        yield
+    finally:
+        numerics.epsfcn /= 10  # reset value

From 87291ef1531bf2213f3bcfb97ece331f51abbdd1 Mon Sep 17 00:00:00 2001
From: james <81617086+je-cook@users.noreply.github.com>
Date: Sat, 13 Jun 2026 10:14:53 +0100
Subject: [PATCH 4/7] OH MY GOD THE DUPLICATION IT BURNS

---
 process/core/io/plot/scans.py | 1 +
 process/core/scan.py          | 1 +
 2 files changed, 2 insertions(+)

diff --git a/process/core/io/plot/scans.py b/process/core/io/plot/scans.py
index d1325f59f4..00dba3ba85 100644
--- a/process/core/io/plot/scans.py
+++ b/process/core/io/plot/scans.py
@@ -36,6 +36,7 @@
 from matplotlib.ticker import MultipleLocator, PercentFormatter
 
 from process.core.io.mfile import MFile
+from process.core.io.mfile.cli import mfile
 from process.core.io.variable_metadata import var_dicts as meta
 from process.core.scan import ScanVariables
 from process.data_structure.numerics import SolverOutputCondition
diff --git a/process/core/scan.py b/process/core/scan.py
index 5876ad0f0e..4f407f7ce9 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -508,3 +508,4 @@ def scan_select(self, nsweep, sweep, iscan) -> ScanVariables:
         """Select a scan"""
         sv = ScanVariables(nsweep)
         sv.set(self.data, sweep[iscan - 1])
+        return sv

From 114d5f9334e4c7ecccb86a4b54669b7067a6a7e6 Mon Sep 17 00:00:00 2001
From: james <81617086+je-cook@users.noreply.github.com>
Date: Mon, 15 Jun 2026 09:21:55 +0100
Subject: [PATCH 5/7] WIP reorganise

---
 process/core/scan.py | 152 +++++++++++++++++++++----------------------
 1 file changed, 75 insertions(+), 77 deletions(-)

diff --git a/process/core/scan.py b/process/core/scan.py
index 4f407f7ce9..16930b2194 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -225,25 +225,23 @@ def get_val(self, mfile, scan):
 
 
 class Scan:
-    """Perform a parameter scan using the Fortran scan module."""
+    """Perform a parameter scan
+
+    Parameters
+    ----------
+    models :
+        Physics and engineering model objects
+    solver :
+        Which solver to use, as specified in solver.py
+    data :
+        Data structure object
+    """
 
     def __init__(self, models: Model, solver: str, data: DataStructure):
-        """Immediately run the run_scan() method.
-
-        Parameters
-        ----------
-        models :
-            Physics and engineering model objects
-        solver :
-            Which solver to use, as specified in solver.py
-        data :
-            Data structure object
-        """
         self.models = models
         self.solver = solver
         self.data = data
         self.solver_handler = SolverHandler(models, solver, data)
-        self.run_scan()
 
     def run_scan(self):
         """Call a solver over a range of values of one of the variables.
@@ -258,44 +256,42 @@ def run_scan(self):
         ProcessValueError
             isweep value greater than IPNSCNS
         """
-        if self.data.scan.isweep == 0:
-            # Solve single problem, rather than an array of problems (scan)
-            # doopt() can also run just an evaluation
-            start_time = time.time()
-            ifail = self.doopt()
-            write_output_files(
-                models=self.models,
-                data=self.data,
-                ifail=ifail,
-                runtime=time.time() - start_time,
-            )
-            show_errors(constants.NOUT)
-            return
-
-        if self.data.scan.isweep > IPNSCNS:
+        if self.data.scan.isweep > 0 and self.data.scan.isweep > IPNSCNS:
             raise ProcessValueError(
                 "Illegal value of isweep",
                 isweep=self.data.scan.isweep,
                 IPNSCNS=IPNSCNS,
             )
+        isw = self.data.scan.isweep or 1
+        self.scan_array = np.arange(isw * (self.data.scan.isweep_2 or 1)).reshape(
+            isw, -1
+        )
+        if self.data.scan.isweep == 0:
+            # Solve single problem, rather than an array of problems (scan)
+            # doopt() can also run just an evaluation
+            self._start_time = time.time()
+            self._ifail = self.solver_handler.run()
+            self._finish_time = time.time()
+            return
 
         if self.data.scan.scan_dim == 2:
             self.scan_2d()
         else:
             self.scan_1d()
 
-    def doopt(self):
-        """Run the optimiser or solver."""
-        ifail = self.solver_handler.run()
+    def write_outputs(self):
+        write_output_files(
+            models=self.models,
+            data=self.data,
+            ifail=self._ifail,
+            runtime=self._finish_time - self._start_time,
+        )
         constraints.constraints_output(self.data, self.solver)
+        show_errors(constants.NOUT)
 
-        return ifail
+    def setup_scan(self, scan_index): ...
 
     def scan_1d(self):
-        """Run a 1-D scan."""
-        # initialise dict which will contain ifail values for each scan point
-        scan_1d_ifail_dict = {}
-
         for iscan in range(1, self.data.scan.isweep + 1):
             self.scan_1d_write_point_header(iscan)
             start_time = time.time()
@@ -311,6 +307,43 @@ def scan_1d(self):
             show_errors(constants.NOUT)
             logging_model_handler.clear_logs()
 
+    def scan_2d(self):
+        iscan = 1
+
+        # initialise array which will contain ifail values for each scan point
+        scan_2d_ifail_list = np.zeros(
+            (NOUTVARS, IPNSCNS),
+            dtype=np.float64,
+            order="F",
+        )
+        for iscan_1 in range(1, self.data.scan.isweep + 1):
+            for iscan_2 in range(1, self.data.scan.isweep_2 + 1):
+                self.scan_2d_write_point_header(iscan, iscan_1, iscan_2)
+                start_time = time.time()
+                ifail = self.doopt()
+                write_output_files(
+                    models=self.models,
+                    data=self.data,
+                    ifail=ifail,
+                    runtime=time.time() - start_time,
+                )
+
+                show_errors(constants.NOUT)
+                logging_model_handler.clear_logs()
+                scan_2d_ifail_list[iscan_1][iscan_2] = ifail
+                iscan += 1
+
+
+def scan_summary(scan):
+
+    def write_outputs(self):
+        self.post_optimise(self._ifail)
+
+    def scan_1d(self):
+        """Run a 1-D scan."""
+        # initialise dict which will contain ifail values for each scan point
+        scan_1d_ifail_dict = {}
+
         # outvar now contains results
         self.scan_1d_write_plot(self.data.scan)
         print("Scan Convergence Summary \n")
@@ -327,7 +360,7 @@ def scan_1d(self):
         ]
         for iscan in range(self.data.scan.isweep):
             pstring = (
-                f"Scan {iscan:02d}: {nsweep_var.name} = {sweep_values[iscan]} "
+                f"Scan {iscan:02d}: {nsweep_var.fname} = {sweep_values[iscan]} "
                 + " " * offsets[iscan]
                 + "\u001b[3{}CONVERGED \u001b[0m"
             )
@@ -340,37 +373,6 @@ def scan_1d(self):
         converged_percentage = converged_count / self.data.scan.isweep * 100
         print(f"\nConvergence Percentage: {converged_percentage:.2f}%")
 
-    def scan_2d(self):
-        """Run a 2-D scan."""
-        # Initialise intent(out) arrays
-        self.scan_2d_init(self.data.scan)
-        iscan = 1
-
-        # initialise array which will contain ifail values for each scan point
-        scan_2d_ifail_list = np.zeros(
-            (NOUTVARS, IPNSCNS),
-            dtype=np.float64,
-            order="F",
-        )
-        for iscan_1 in range(1, self.data.scan.isweep + 1):
-            for iscan_2 in range(1, self.data.scan.isweep_2 + 1):
-                self.scan_2d_write_point_header(iscan, iscan_1, iscan_2)
-                start_time = time.time()
-                ifail = self.doopt()
-                write_output_files(
-                    models=self.models,
-                    data=self.data,
-                    ifail=ifail,
-                    runtime=time.time() - start_time,
-                )
-
-                show_errors(constants.NOUT)
-                logging_model_handler.clear_logs()
-                scan_2d_ifail_list[iscan_1][iscan_2] = ifail
-                iscan += 1
-
-        self.output_2d_summary(scan_2d_ifail_list)
-
     def scan_2d_init(self):
         sv = self.data.scan
         for d, n, v in (
@@ -384,9 +386,7 @@ def scan_2d_init(self):
             process_output.ovarin(constants.MFILE, d, n, v)
 
     def output_2d_summary(self, scan_2d_ifail_list):
-        print(
-            " ****************************************** Scan Convergence Summary ****************************************** \n"
-        )
+        print("Scan Convergence Summary\n")
         sweep_1_values = self.data.scan.sweep[: self.data.scan.isweep]
         sweep_2_values = self.data.scan.sweep_2[: self.data.scan.isweep_2]
         nsweep_var = self.scan_select(
@@ -415,8 +415,8 @@ def output_2d_summary(self, scan_2d_ifail_list):
         for iscan_1 in range(1, self.data.scan.isweep + 1):
             for iscan_2 in range(1, self.data.scan.isweep_2 + 1):
                 string = (
-                    f"Scan {scan_point:02d}: ({nsweep_var.name} = "
-                    f"{sweep_1_values[iscan_1 - 1]}, {nsweep_2_var.name} "
+                    f"Scan {scan_point:02d}: ({nsweep_var.fname} = "
+                    f"{sweep_1_values[iscan_1 - 1]}, {nsweep_2_var.fname} "
                     f"= {sweep_2_values[iscan_2 - 1]}) "
                     + " " * offsets[iscan_1 - 1][iscan_2 - 1]
                     + "\u001b[3{}CONVERGED \u001b[0m"
@@ -452,8 +452,7 @@ def scan_1d_write_point_header(self, iscan: int):
         process_output.write(
             constants.NOUT,
             f"Scan point {iscan} of {self.data.scan.isweep} : {self.data.globals.xlabel}"
-            f", {self.data.globals.vlabel} = {self.data.scan.sweep[iscan - 1]} "
-            "",
+            f", {self.data.globals.vlabel} = {self.data.scan.sweep[iscan - 1]} ",
         )
         process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan)
 
@@ -480,8 +479,7 @@ def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2):
             constants.NOUT,
             f"2D Scan point {iscan} of {self.data.scan.isweep * self.data.scan.isweep_2} : "
             f"{self.data.globals.vlabel} = {self.data.scan.sweep[iscan_1 - 1]} and"
-            f" {self.data.globals.vlabel_2} = {self.data.scan.sweep_2[iscan_r - 1]} "
-            "",
+            f" {self.data.globals.vlabel_2} = {self.data.scan.sweep_2[iscan_r - 1]} ",
         )
         process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan)
 

From 20b5d17011c62f4646bbe09fb3b44f02a01ab3c2 Mon Sep 17 00:00:00 2001
From: james <81617086+je-cook@users.noreply.github.com>
Date: Mon, 15 Jun 2026 09:28:37 +0100
Subject: [PATCH 6/7] remove unnecessary variable

---
 process/core/scan.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/process/core/scan.py b/process/core/scan.py
index 16930b2194..1e175cd548 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -373,6 +373,13 @@ def scan_1d(self):
         converged_percentage = converged_count / self.data.scan.isweep * 100
         print(f"\nConvergence Percentage: {converged_percentage:.2f}%")
 
+    def scan_2d(self):
+        """Run a 2-D scan."""
+        # Initialise intent(out) arrays
+        self.scan_2d_init()
+
+        self.output_2d_summary(scan_2d_ifail_list)
+
     def scan_2d_init(self):
         sv = self.data.scan
         for d, n, v in (

From b449ccd31ebbc6bc19fa2d4eb8b7a42009116d37 Mon Sep 17 00:00:00 2001
From: james <81617086+je-cook@users.noreply.github.com>
Date: Mon, 29 Jun 2026 09:12:50 +0100
Subject: [PATCH 7/7] sudo code mvp

---
 process/core/input.py                      |  13 +-
 process/core/io/plot/scans.py              |   1 -
 process/core/scan.py                       | 398 +++++++++------------
 process/core/solver/iteration_variables.py |   5 +-
 process/core/solver/solver.py              |  12 -
 process/data_structure/global_variables.py |  18 +-
 process/data_structure/scan_variables.py   |  67 ++--
 process/main.py                            |   1 +
 8 files changed, 221 insertions(+), 294 deletions(-)

diff --git a/process/core/input.py b/process/core/input.py
index 93886b8a16..af1f485909 100644
--- a/process/core/input.py
+++ b/process/core/input.py
@@ -1133,24 +1133,15 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
         "scan",
         int,
         choices=range(IPNSCNS + 1),
+        array=True,
     ),
     "nsweep": InputVariable(
         "scan",
         int,
         choices=range(1, IPNSCNV + 1),
-    ),
-    "isweep_2": InputVariable(
-        "scan",
-        int,
-        choices=range(IPNSCNS + 1),
-    ),
-    "nsweep_2": InputVariable(
-        "scan",
-        int,
-        choices=range(1, IPNSCNV + 1),
+        array=True,
     ),
     "sweep": InputVariable("scan", float, array=True),
-    "sweep_2": InputVariable("scan", float, array=True),
     "impvardiv": InputVariable(
         "reinke",
         int,
diff --git a/process/core/io/plot/scans.py b/process/core/io/plot/scans.py
index 00dba3ba85..d1325f59f4 100644
--- a/process/core/io/plot/scans.py
+++ b/process/core/io/plot/scans.py
@@ -36,7 +36,6 @@
 from matplotlib.ticker import MultipleLocator, PercentFormatter
 
 from process.core.io.mfile import MFile
-from process.core.io.mfile.cli import mfile
 from process.core.io.variable_metadata import var_dicts as meta
 from process.core.scan import ScanVariables
 from process.data_structure.numerics import SolverOutputCondition
diff --git a/process/core/scan.py b/process/core/scan.py
index 1e175cd548..cca39075de 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -18,11 +18,11 @@
 from process.core.log import logging_model_handler, show_errors
 from process.core.solver import constraints
 from process.core.solver.solver_handler import SolverHandler
-from process.data_structure.scan_variables import IPNSCNS, NOUTVARS
 from process.models.availability import AvailabilityModel
 
 if TYPE_CHECKING:
-    from process.core.model import DataStructure, Model
+    from process.core.model import DataStructure
+    from process.main import Models
 
 
 logger = logging.getLogger(__name__)
@@ -224,6 +224,13 @@ def get_val(self, mfile, scan):
     n_tf_coil_turns = (81, Area.T)
 
 
+@dataclass
+class ScanRes:
+    iscan: int
+    ifail: int
+    solver: SolverHandler
+
+
 class Scan:
     """Perform a parameter scan
 
@@ -237,15 +244,13 @@ class Scan:
         Data structure object
     """
 
-    def __init__(self, models: Model, solver: str, data: DataStructure):
+    def __init__(self, models: Models, solver: str, data: DataStructure):
         self.models = models
         self.solver = solver
         self.data = data
-        self.solver_handler = SolverHandler(models, solver, data)
-
-    def run_scan(self):
-        """Call a solver over a range of values of one of the variables.
 
+    def _run(self, iscan, nsweep, sweep, data):
+        """
         This method calls the optimisation routine VMCON a number of times, by
         performing a sweep over a range of values of a particular variable. A
         number of output variable values are written to the MFILE.DAT file at
@@ -256,23 +261,16 @@ def run_scan(self):
         ProcessValueError
             isweep value greater than IPNSCNS
         """
-        if self.data.scan.isweep > 0 and self.data.scan.isweep > IPNSCNS:
-            raise ProcessValueError(
-                "Illegal value of isweep",
-                isweep=self.data.scan.isweep,
-                IPNSCNS=IPNSCNS,
-            )
-        isw = self.data.scan.isweep or 1
-        self.scan_array = np.arange(isw * (self.data.scan.isweep_2 or 1)).reshape(
-            isw, -1
-        )
-        if self.data.scan.isweep == 0:
-            # Solve single problem, rather than an array of problems (scan)
-            # doopt() can also run just an evaluation
-            self._start_time = time.time()
-            self._ifail = self.solver_handler.run()
-            self._finish_time = time.time()
-            return
+        sh = SolverHandler(self.models, self.solver, data)
+        # TODO queue the output to avoid race condition (?)
+        if data.scan.nsweep is not None:
+            self.write_point_header(iscan)
+        start_time = time.time()
+        ifail = sh.run()
+        end_time = time.time() - start_time
+        write_output_files(models=self.models, data=data, ifail=ifail, runtime=end_time)
+        nums = data.numerics
+        nums.sqsumsq = sum(r**2 for r in nums.rcm[: nums.neqns]) ** 0.5
 
         if self.data.scan.scan_dim == 2:
             self.scan_2d()
@@ -286,171 +284,23 @@ def write_outputs(self):
             ifail=self._ifail,
             runtime=self._finish_time - self._start_time,
         )
-        constraints.constraints_output(self.data, self.solver)
         show_errors(constants.NOUT)
 
-    def setup_scan(self, scan_index): ...
-
-    def scan_1d(self):
-        for iscan in range(1, self.data.scan.isweep + 1):
-            self.scan_1d_write_point_header(iscan)
-            start_time = time.time()
-            ifail = self.doopt()
-            scan_1d_ifail_dict[iscan] = ifail
-            write_output_files(
-                models=self.models,
-                data=self.data,
-                ifail=ifail,
-                runtime=time.time() - start_time,
-            )
-
-            show_errors(constants.NOUT)
-            logging_model_handler.clear_logs()
+        logging_model_handler.clear_logs()
+        optimisation_output(data)
+        constraints.constraints_output(data, self.solver)
 
-    def scan_2d(self):
-        iscan = 1
+        return ScanRes(iscan, ifail, sh)
 
-        # initialise array which will contain ifail values for each scan point
-        scan_2d_ifail_list = np.zeros(
-            (NOUTVARS, IPNSCNS),
-            dtype=np.float64,
-            order="F",
-        )
-        for iscan_1 in range(1, self.data.scan.isweep + 1):
-            for iscan_2 in range(1, self.data.scan.isweep_2 + 1):
-                self.scan_2d_write_point_header(iscan, iscan_1, iscan_2)
-                start_time = time.time()
-                ifail = self.doopt()
-                write_output_files(
-                    models=self.models,
-                    data=self.data,
-                    ifail=ifail,
-                    runtime=time.time() - start_time,
-                )
-
-                show_errors(constants.NOUT)
-                logging_model_handler.clear_logs()
-                scan_2d_ifail_list[iscan_1][iscan_2] = ifail
-                iscan += 1
-
-
-def scan_summary(scan):
-
-    def write_outputs(self):
-        self.post_optimise(self._ifail)
-
-    def scan_1d(self):
-        """Run a 1-D scan."""
-        # initialise dict which will contain ifail values for each scan point
-        scan_1d_ifail_dict = {}
-
-        # outvar now contains results
-        self.scan_1d_write_plot(self.data.scan)
-        print("Scan Convergence Summary \n")
-        sweep_values = self.data.scan.sweep[: self.data.scan.isweep]
-        nsweep_var = self.scan_select(
-            self.data.scan.nsweep, self.data.scan.sweep, self.data.scan.isweep
-        )
-        converged_count = 0
-        # offsets for aligning the converged/unconverged column
-        max_sweep_value_length = len(str(np.max(sweep_values)).replace(".", ""))
-        offsets = [
-            max_sweep_value_length - len(str(sweep_val).replace(".", ""))
-            for sweep_val in sweep_values
+    def _set_v_x_label(self, iscan: list[int]):
+        sv = [
+            self.scan_select(self.data.scan.nsweep, self.data.scan.sweep, isc)
+            for isc in iscan
         ]
-        for iscan in range(self.data.scan.isweep):
-            pstring = (
-                f"Scan {iscan:02d}: {nsweep_var.fname} = {sweep_values[iscan]} "
-                + " " * offsets[iscan]
-                + "\u001b[3{}CONVERGED \u001b[0m"
-            )
-            if scan_1d_ifail_dict[iscan + 1] == 1:
-                converged_count += 1
-                pstring.format("2m")
-            else:
-                pstring.format("1mUN")
-            print(pstring)
-        converged_percentage = converged_count / self.data.scan.isweep * 100
-        print(f"\nConvergence Percentage: {converged_percentage:.2f}%")
-
-    def scan_2d(self):
-        """Run a 2-D scan."""
-        # Initialise intent(out) arrays
-        self.scan_2d_init()
-
-        self.output_2d_summary(scan_2d_ifail_list)
-
-    def scan_2d_init(self):
-        sv = self.data.scan
-        for d, n, v in (
-            ("Number of first variable scan points", "(isweep)", sv.isweep),
-            ("Number of second variable scan points", "(isweep_2)", sv.isweep_2),
-            ("Scanning first variable number", "(nsweep)", sv.nsweep),
-            ("Scanning second variable number", "(nsweep_2)", sv.nsweep_2),
-            ("Scanning second variable number", "(nsweep_2)", sv.nsweep_2),
-            ("Scanning second variable number", "(nsweep_2)", sv.nsweep_2),
-        ):
-            process_output.ovarin(constants.MFILE, d, n, v)
-
-    def output_2d_summary(self, scan_2d_ifail_list):
-        print("Scan Convergence Summary\n")
-        sweep_1_values = self.data.scan.sweep[: self.data.scan.isweep]
-        sweep_2_values = self.data.scan.sweep_2[: self.data.scan.isweep_2]
-        nsweep_var = self.scan_select(
-            self.data.scan.nsweep, self.data.scan.sweep, self.data.scan.isweep
-        )
-        nsweep_2_var = self.scan_select(
-            self.data.scan.nsweep_2, self.data.scan.sweep_2, self.data.scan.isweep_2
-        )
-        converged_count = 0
-        scan_point = 1
-        # offsets for aligning the converged/unconverged column
-        max_sweep1_value_length = len(str(np.max(sweep_1_values)).replace(".", ""))
-        max_sweep2_value_length = len(str(np.max(sweep_2_values)).replace(".", ""))
-        offsets = np.zeros(
-            (self.data.scan.isweep, self.data.scan.isweep_2), dtype=int, order="F"
-        )
-        for count1, sweep1 in enumerate(sweep_1_values):
-            for count2, sweep2 in enumerate(sweep_2_values):
-                offsets[count1][count2] = (
-                    max_sweep1_value_length
-                    - len(str(sweep1).replace(".", ""))
-                    + max_sweep2_value_length
-                    - len(str(sweep2).replace(".", ""))
-                )
-
-        for iscan_1 in range(1, self.data.scan.isweep + 1):
-            for iscan_2 in range(1, self.data.scan.isweep_2 + 1):
-                string = (
-                    f"Scan {scan_point:02d}: ({nsweep_var.fname} = "
-                    f"{sweep_1_values[iscan_1 - 1]}, {nsweep_2_var.fname} "
-                    f"= {sweep_2_values[iscan_2 - 1]}) "
-                    + " " * offsets[iscan_1 - 1][iscan_2 - 1]
-                    + "\u001b[3{}CONVERGED \u001b[0m"
-                )
-                if scan_2d_ifail_list[iscan_1][iscan_2] == 1:
-                    converged_count += 1
-                    print(string.format("2m"))
-                else:
-                    print(string.format("1mUN"))
-                scan_point += 1
-        converged_percentage = (
-            converged_count / (self.data.scan.isweep * self.data.scan.isweep_2) * 100
-        )
-        print(f"\nConvergence Percentage: {converged_percentage:.2f}%")
+        self.data.globals.vlabel = [s.fname for s in sv]
+        self.data.globals.xlabel = [s.data.description for s in sv]
 
-    def _set_v_x_label(self, iscan, twod=False):
-        sv = (
-            self.scan_select(self.data.scan.nsweep_2, self.data.scan.sweep_2, iscan)
-            if twod
-            else self.scan_select(self.data.scan.nsweep, self.data.scan.sweep, iscan)
-        )
-        self.data.globals.vlabel = sv.fname
-        self.data.globals.xlabel = sv.description
-
-    def scan_1d_write_point_header(self, iscan: int):
-        """Scan 1d header"""
-        self.data.globals.iscan_global = iscan
+    def write_point_header(self, iscan):
         self._set_v_x_label(iscan)
 
         process_output.oblnkl(constants.NOUT)
@@ -458,59 +308,163 @@ def scan_1d_write_point_header(self, iscan: int):
 
         process_output.write(
             constants.NOUT,
-            f"Scan point {iscan} of {self.data.scan.isweep} : {self.data.globals.xlabel}"
-            f", {self.data.globals.vlabel} = {self.data.scan.sweep[iscan - 1]} ",
+            f"Scan point {iscan} of {np.prod(self.data.scan.isweep)} : \n".join(
+                f"{v} = {self.data.scan.sweep[iscan[no] - 1]}"
+                for no, v in enumerate(self.data.globals.vlabel)
+            ),
         )
         process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan)
 
         print(
-            f"Starting scan point {iscan} of {self.data.scan.isweep} : "
-            f"{self.data.globals.xlabel} , {self.data.globals.vlabel}"
-            f" = {self.data.scan.sweep[iscan - 1]}"
+            f"Starting scan point {iscan}:  {self.data.globals.xlabel}, \n".join(
+                f"{v} = {self.data.scan.sweep[iscan[no] - 1]}"
+                for no, v in enumerate(self.data.globals.vlabel)
+            )
         )
 
-    def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2):
-        """Scan 2d header"""
-        iscan_r = self.data.scan.isweep_2 - iscan_2 + 1 if iscan_1 % 2 == 0 else iscan_2
-
-        # Makes iscan available globally (read-only)
-        self.data.globals.iscan_global = iscan
+    def scan_select(self, nsweep, sweep, iscan):
+        sv = ScanVariables(nsweep)
+        sv.set(self.data, sweep[iscan - 1])
+        return sv
 
-        self._set_v_x_label(iscan_1)
-        self._set_v_x_label(iscan_r)
+    def run(self):
+        """Call a solver over a range of values of one of the variables.
 
-        process_output.oblnkl(constants.NOUT)
-        process_output.oblnkl(constants.MFILE)
+        This method calls the optimisation routine VMCON a number of times, by
+        performing a sweep over a range of values of a particular variable. A
+        number of output variable values are written to the MFILE.DAT file at
+        each scan point, for plotting or other post-processing purposes.
+        """
+        # vectorise running of self._run
+        if self.data.scan.nsweep is not None:
+            for d, n, v in (
+                ("Number of scan points", "(isweep)", self.data.scan.isweep),
+                ("Scanning variable number", "(nsweep)", self.data.scan.nsweep),
+            ):
+                process_output.ovarin(constants.MFILE, d, n, v)
 
-        process_output.write(
-            constants.NOUT,
-            f"2D Scan point {iscan} of {self.data.scan.isweep * self.data.scan.isweep_2} : "
-            f"{self.data.globals.vlabel} = {self.data.scan.sweep[iscan_1 - 1]} and"
-            f" {self.data.globals.vlabel_2} = {self.data.scan.sweep_2[iscan_r - 1]} ",
+        # TODO copy of self.data for each vectorised run (?)
+        scan_res = np.vectorise(self._run)(
+            self.data.scan.isweep, self.data.scan.nsweep, self.data.scan.sweep, self.data
         )
-        process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan)
+
+        if self.data.scan.nsweep is not None:
+            self.summary(scan_res)
+
+    def summary(self, scan_res):
+        print("Scan Convergence Summary\n")
+        sweep_values = self.data.scan.sweep
+        nsweep_var = [ScanVariables(nsw) for nsw in self.data.scan.nsweep]
+
+        conv_list = []
+        converged_count = 0
+        conv_str = "\u001b[3{}CONVERGED \u001b[0m"
+        for no, sr in enumerate(scan_res):
+            if sr.ifail == 1:
+                converged_count += 1
+                conv = conv_str.format("2m")
+            else:
+                conv = conv_str.format("1mUN")
+            conv_list.append([
+                "{sr.iscan:02d}",
+                nsweep_var[no].fname,
+                sweep_values[sr.iscan],
+                conv,
+            ])
 
         print(
-            f"Starting scan point {iscan}:  {self.data.globals.xlabel}, "
-            f"{self.data.globals.vlabel} = {self.data.scan.sweep[iscan_1 - 1]}"
-            f" and {self.data.globals.xlabel_2}, "
-            f"{self.data.globals.vlabel_2} = {self.data.scan.sweep_2[iscan_r - 1]} "
+            tabulate(conv_list, headers=["Iscan", "Sweep Var", "Sweep Val", "Converged"])
         )
 
-        return iscan_r
+        converged_percentage = converged_count / np.prod(self.data.scan.isweep) * 100
+        print(f"\nConvergence Percentage: {converged_percentage:.2f}%")
 
-    def scan_1d_write_plot(self):
-        if self.data.scan.first_call_1d:
-            for d, n, v in (
-                ("Number of scan points", "(isweep)", self.data.scan.isweep),
-                ("Scanning variable number", "(nsweep)", self.data.scan.nsweep),
-            ):
-                process_output.ovarin(constants.MFILE, d, n, v)
 
-            self.data.scan.first_call_1d = False
+def optimisation_output(data: DataStructure):
+    nums = data.numerics
 
-    def scan_select(self, nsweep, sweep, iscan) -> ScanVariables:
-        """Select a scan"""
-        sv = ScanVariables(nsweep)
-        sv.set(self.data, sweep[iscan - 1])
-        return sv
+    written_warning = False
+
+    # Output optimisation parameters
+    solution_vector_table = []
+    for i in range(nums.nvar):
+        nums.xcs[i] = nums.xcm[i] * nums.scafc[i]
+
+        name = nums.lablxc[nums.ixc[i] - 1]
+        solution_vector_table.append([name, nums.xcs[i], nums.xcm[i]])
+
+        xminn = 1.01 * nums.itv_scaled_lower_bounds[i]
+        xmaxx = 0.99 * nums.itv_scaled_upper_bounds[i]
+
+        # Write to output file if close to optimisation parameter bounds
+        if nums.xcm[i] < xminn or nums.xcm[i] > xmaxx:
+            if not written_warning:
+                written_warning = True
+                process_output.ocmmnt(
+                    constants.NOUT,
+                    (
+                        "Certain operating limits have been reached,"
+                        "\n as shown by the following optimisation parameters that are"
+                        "\n at or near to the edge of their prescribed range :\n"
+                    ),
+                )
+
+            xcval = nums.xcm[i] * nums.scafc[i]
+
+            if nums.xcm[i] < xminn:
+                location, bound = "below", "lower"
+                bounds = nums.itv_scaled_lower_bounds
+            else:
+                location, bound = "above", "upper"
+                bounds = nums.itv_scaled_upper_bounds
+            process_output.write(
+                constants.NOUT,
+                f"   {name:<30}= {xcval} is at or {location} its {bound} bound:"
+                f" {bounds[i] * nums.scafc[i]}",
+            )
+
+        xnorm = (
+            1.0
+            if nums.boundu[i] == nums.boundl[i]
+            else min(
+                max(
+                    (nums.xcm[i] - nums.itv_scaled_lower_bounds[i])
+                    / (
+                        nums.itv_scaled_upper_bounds[i] - nums.itv_scaled_lower_bounds[i]
+                    ),
+                    0.0,
+                ),
+                1.0,
+            )
+        )
+
+        # Write optimisation parameters to mfile
+        for d, var, v in (
+            (nums.lablxc[nums.ixc[i] - 1], f"(itvar{i + 1:03d})", nums.xcs[i]),
+            (f"{name} (final value/initial value)", f"(xcm{i + 1:03d})", nums.xcm[i]),
+            (f"{name} (range normalised)", f"(nitvar{i + 1:03d})", xnorm),
+            (
+                f"{name} (upper bound)",
+                f"(boundu{i + 1:03d})",
+                nums.itv_scaled_upper_bounds[i] * nums.scafc[i],
+            ),
+            (
+                f"{name} (lower bound)",
+                f"(boundl{i + 1:03d})",
+                nums.itv_scaled_lower_bounds[i] * nums.scafc[i],
+            ),
+        ):
+            process_output.ovarre(constants.MFILE, d, var, v)
+
+    # Write optimisation parameter headings to output file
+    process_output.osubhd(
+        constants.NOUT, "The solution vector is comprised as follows :"
+    )
+    process_output.write(
+        constants.NOUT,
+        tabulate(
+            solution_vector_table,
+            headers=["", "Final value", "Final / initial"],
+            numalign="left",
+        ),
+    )
diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py
index c40f932ec1..85b0c07d4c 100644
--- a/process/core/solver/iteration_variables.py
+++ b/process/core/solver/iteration_variables.py
@@ -326,10 +326,7 @@ def load_iteration_variables(data):
 
         # warn of the iteration variable is also a scan variable because this will cause
         # the optimiser and scan to overwrite the same variable and conflict
-        if iteration_variable.name in {
-            data.globals.vlabel,
-            data.globals.vlabel_2,
-        }:
+        if iteration_variable.name in data.globals.vlabel:
             logger.critical(
                 (
                     "The sweep variable is also an iteration variable and will be "
diff --git a/process/core/solver/solver.py b/process/core/solver/solver.py
index 269177e338..513bebeafa 100644
--- a/process/core/solver/solver.py
+++ b/process/core/solver/solver.py
@@ -325,11 +325,6 @@ def verror(self):
                     "A feasible solution may be difficult to achieve.",
                     "Try changing or adding variables to IXC.",
                 ),
-                4: (
-                    "An uphill search direction was found.",
-                    "Try changing the equations in ICC, or",
-                    "adding new variables to IXC.",
-                ),
                 SolverOutputCondition.NO_SOLUTION: (
                     "The quadratic programming technique was unable to",
                     "find a feasible point.\n",
@@ -337,13 +332,6 @@ def verror(self):
                     "their initial values (especially if only 1 optimisation",
                     "iteration was performed).",
                 ),
-                6: (
-                    "The quadratic programming technique was restricted",
-                    "by an artificial bound, or failed due to a singular",
-                    "matrix.",
-                    "Try changing the equations in ICC, or",
-                    "adding new variables to IXC.",
-                ),
             }.get(self.info, "Unknown Error code")
         )
 
diff --git a/process/data_structure/global_variables.py b/process/data_structure/global_variables.py
index f00c9441ab..a84275727e 100644
--- a/process/data_structure/global_variables.py
+++ b/process/data_structure/global_variables.py
@@ -1,6 +1,6 @@
 """Module containing variables necessary for running PROCESS"""
 
-from dataclasses import dataclass
+from dataclasses import dataclass, field
 
 
 @dataclass(slots=True)
@@ -22,19 +22,11 @@ class GlobalData:
     output_prefix: str = ""
     """Output file path prefix"""
 
-    xlabel: str = ""
-    """Scan parameter description label"""
+    xlabel: list[str] = field(default_factory=lambda: [""])
+    """Scan parameters description label"""
 
-    vlabel: str = ""
-    """Scan value name label"""
-
-    xlabel_2: str = ""
-    """Scan parameter description label (2nd dimension)"""
-
-    vlabel_2: str = ""
-    """Scan value name label (2nd dimension)"""
-
-    iscan_global: int = 0
+    vlabel: list[str] = field(default_factory=lambda: [""])
+    """Scan values name label"""
 
     convergence_parameter: float = 0.0
     """VMCON convergence parameter 'sum'"""
diff --git a/process/data_structure/scan_variables.py b/process/data_structure/scan_variables.py
index 0a268e81e1..18d06ccf4e 100644
--- a/process/data_structure/scan_variables.py
+++ b/process/data_structure/scan_variables.py
@@ -4,17 +4,12 @@
 
 import numpy as np
 
+from process.core.exceptions import ProcessValueError
+
 IPNSCNS = 1000
 """Maximum number of scan points"""
 
 
-IPNSCNV = 81
-"""Number of available scan variables"""
-
-
-NOUTVARS = 84
-
-
 @dataclass(slots=True)
 class ScanData:
     """Dataclass holding scan variables"""
@@ -22,35 +17,45 @@ class ScanData:
     scan_dim: int = 1
     """1-D or 2-D scan switch (1=1D, 2=2D)"""
 
-    isweep: int = 0
+    isweep: list[int] | int = 1
     """Number of scan points to calculate"""
 
-    isweep_2: int = 0
-    """Number of 2D scan points to calculate"""
-
-    nsweep: int = 1
+    nsweep: list[int] | int | None = None
     """Switch denoting quantity to scan
+    see `process.core.scan.ScanVariables` for available options
     """
 
-    nsweep_2: int = 3
-    """nsweep_2 /3/ : switch denoting quantity to scan for 2D scan:"""
-
-    sweep: list[float] = field(
-        default_factory=lambda: np.zeros(IPNSCNS, dtype=np.float64)
-    )
-    """sweep(IPNSCNS) /../: actual values to use in scan"""
-
-    sweep_2: list[float] = field(
-        default_factory=lambda: np.zeros(IPNSCNS, dtype=np.float64)
-    )
-    """sweep_2(IPNSCNS) /../: actual values to use in 2D scan"""
-
-    # Vars in subroutines scan_1d and scan_2d requiring re-initialising before
-    # each new run
-
-    first_call_1d: bool = True
-
-    first_call_2d: bool = True
+    sweep: np.ndarray = field(default_factory=lambda: np.zeros(1, dtype=np.float64))
+    """Actual values to use in scan"""
+
+    def __post_init__(self):
+        if isinstance(self.isweep, int):
+            # avoid old 0 default
+            self.isweep = [self.isweep or 1]
+
+        if len(self.isweep) > 2 or len(self.sweep.shape) > 2:
+            raise NotImplementedError("N-D Scans not currently supported")
+
+        if max(self.isweep) > IPNSCNS:
+            raise ProcessValueError(
+                "Illegal value of isweep",
+                isweep=self.isweep,
+                IPNSCNS=IPNSCNS,
+            )
+        if self.nsweep != len(self.isweep):
+            raise ValueError(
+                "Number of sweep variables not equal to scan point dimensions"
+            )
+        if self.sweep.shape != self.isweep:
+            if self.isweep != 1:
+                self.isweep = list(self.sweep.shape)
+            else:
+                print("Unset sweep values set to zero")
+                # TODO append to size instead of resetting
+                self.sweep = np.zeros(self.isweep, dtype=np.float64)
+
+        self.nsweep = np.asarray(self.nsweep, dtype=int)
+        self.isweep = np.asarray(self.isweep, dtype=int)
 
 
 CREATE_DICTS_FROM_DATACLASS = ScanData
diff --git a/process/main.py b/process/main.py
index 2501505685..2f326b924c 100644
--- a/process/main.py
+++ b/process/main.py
@@ -459,6 +459,7 @@ def run_scan(self):
                 "select either 1 (optimise) or -2 (no optimisation)."
             )
         self.scan = Scan(self.models, self.solver, self.data)
+        self.scan.run()
 
     @staticmethod
     def show_errors():