From 495e7f1603c127fc7fa2b5d54859ca53118327b7 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 | 378 ++++++++++++-------------------------------
 1 file changed, 104 insertions(+), 274 deletions(-)

diff --git a/process/core/scan.py b/process/core/scan.py
index e74a960c0d..276ef1cbdc 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -308,20 +308,17 @@ def post_optimise(self, ifail: int):
         )
 
         process_output.oheadr(constants.NOUT, "Numerics")
-        if self.solver == "fsolve":
-            process_output.ocmmnt(
-                constants.NOUT, "PROCESS has performed an fsolve (evaluation) run."
-            )
-        else:
-            process_output.ocmmnt(
-                constants.NOUT, "PROCESS has performed a VMCON (optimisation) run."
-            )
+        process_output.ocmmnt(
+            constants.NOUT,
+            f"PROCESS has performed a {'fsolve' if self.solver == 'fsolve' else 'VMCON'}"
+            " (optimisation) run.",
+        )
         if ifail != 1:
             process_output.ovarre(constants.NOUT, "Error flag", "(ifail)", ifail)
             process_output.oheadr(
                 constants.IOTTY, "PROCESS COULD NOT FIND A FEASIBLE SOLUTION"
             )
-            process_output.oblnkl(constants.IOTTY)
+            print()
 
             logger.critical("Solver returns with ifail /= 1. %s", ifail)
 
@@ -329,7 +326,7 @@ def post_optimise(self, ifail: int):
             if self.solver == "vmcon":
                 self.verror(ifail)
             process_output.oblnkl(constants.NOUT)
-            process_output.oblnkl(constants.IOTTY)
+            print()
         else:
             # Solution found
             if self.solver != "fsolve":
@@ -350,33 +347,13 @@ def post_optimise(self, ifail: int):
             process_output.ovarre(constants.NOUT, "Error flag", "(ifail)", ifail)
 
             if self.data.numerics.sqsumsq >= 1.0e-2:
-                process_output.oblnkl(constants.NOUT)
-                process_output.ocmmnt(
-                    constants.NOUT,
-                    "WARNING: Constraint residues are HIGH; consider re-running",
-                )
-                process_output.ocmmnt(
-                    constants.NOUT,
-                    "   with lower values of EPSVMC to confirm convergence...",
-                )
-                process_output.ocmmnt(
-                    constants.NOUT,
-                    "   (should be able to get down to about 1.0E-8 okay)",
-                )
-                process_output.oblnkl(constants.NOUT)
-                process_output.ocmmnt(
-                    constants.IOTTY,
-                    "WARNING: Constraint residues are HIGH; consider re-running",
-                )
-                process_output.ocmmnt(
-                    constants.IOTTY,
-                    "   with lower values of EPSVMC to confirm convergence...",
-                )
-                process_output.ocmmnt(
-                    constants.IOTTY,
-                    "   (should be able to get down to about 1.0E-8 okay)",
+                string = (
+                    "WARNING: Constraint residues are HIGH; consider re-running\n"
+                    "   with lower values of EPSVMC to confirm convergence...\n"
+                    "   (should be able to get down to about 1.0E-8 okay)\n"
                 )
-                process_output.oblnkl(constants.IOTTY)
+                process_output.ocmmnt(constants.NOUT, ("\n" + string))
+                print(string)
 
                 logger.warning(
                     f"High final constraint residues. {self.data.numerics.sqsumsq=}"
@@ -753,149 +730,58 @@ def verror(ifail: int):
 
         """
         if ifail == -1:
-            process_output.ocmmnt(constants.NOUT, "User-terminated execution of VMCON.")
-            process_output.ocmmnt(constants.IOTTY, "User-terminated execution of VMCON.")
+            strings = ("User-terminated execution of VMCON.",)
         elif ifail == 0:
-            process_output.ocmmnt(
-                constants.NOUT, "Improper input parameters to the VMCON routine."
+            strings = (
+                "Improper input parameters to the VMCON routine.",
+                "PROCESS coding must be checked.",
             )
-            process_output.ocmmnt(constants.NOUT, "PROCESS coding must be checked.")
-
-            process_output.ocmmnt(
-                constants.IOTTY, "Improper input parameters to the VMCON routine."
-            )
-            process_output.ocmmnt(constants.IOTTY, "PROCESS coding must be checked.")
         elif ifail == 2:
-            process_output.ocmmnt(
-                constants.NOUT,
-                "The maximum number of calls has been reached without solution.",
-            )
-            process_output.ocmmnt(
-                constants.NOUT,
-                "The code may be stuck in a minimum in the residual space that is "
-                "significantly above zero.",
-            )
-            process_output.oblnkl(constants.NOUT)
-            process_output.ocmmnt(
-                constants.NOUT, "There is either no solution possible, or the code"
-            )
-            process_output.ocmmnt(
-                constants.NOUT, "is failing to escape from a deep local minimum."
-            )
-            process_output.ocmmnt(
-                constants.NOUT,
-                "Try changing the variables in IXC, or modify their initial values.",
-            )
-
-            process_output.ocmmnt(
-                constants.IOTTY,
+            strings = (
                 "The maximum number of calls has been reached without solution.",
-            )
-            process_output.ocmmnt(
-                constants.IOTTY,
-                "The code may be stuck in a minimum in the residual space that is "
-                "significantly above zero.",
-            )
-            process_output.oblnkl(constants.NOUT)
-            process_output.oblnkl(constants.IOTTY)
-            process_output.ocmmnt(
-                constants.IOTTY, "There is either no solution possible, or the code"
-            )
-            process_output.ocmmnt(
-                constants.IOTTY, "is failing to escape from a deep local minimum."
-            )
-            process_output.ocmmnt(
-                constants.IOTTY,
+                (
+                    "The code may be stuck in a minimum in the residual space that is "
+                    "significantly above zero."
+                ),
+                "",
+                "There is either no solution possible, or the code",
+                "is failing to escape from a deep local minimum.",
                 "Try changing the variables in IXC, or modify their initial values.",
             )
         elif ifail == 3:
-            process_output.ocmmnt(
-                constants.NOUT, "The line search required the maximum of 10 calls."
-            )
-            process_output.ocmmnt(
-                constants.NOUT, "A feasible solution may be difficult to achieve."
-            )
-            process_output.ocmmnt(
-                constants.NOUT, "Try changing or adding variables to IXC."
-            )
-
-            process_output.ocmmnt(
-                constants.IOTTY, "The line search required the maximum of 10 calls."
-            )
-            process_output.ocmmnt(
-                constants.IOTTY, "A feasible solution may be difficult to achieve."
-            )
-            process_output.ocmmnt(
-                constants.IOTTY, "Try changing or adding variables to IXC."
+            strings = (
+                "The line search required the maximum of 10 calls.",
+                "A feasible solution may be difficult to achieve.",
+                "Try changing or adding variables to IXC.",
             )
         elif ifail == 4:
-            process_output.ocmmnt(
-                constants.NOUT, "An uphill search direction was found."
-            )
-            process_output.ocmmnt(
-                constants.NOUT, "Try changing the equations in ICC, or"
-            )
-            process_output.ocmmnt(constants.NOUT, "adding new variables to IXC.")
-
-            process_output.ocmmnt(
-                constants.IOTTY, "An uphill search direction was found."
-            )
-            process_output.ocmmnt(
-                constants.IOTTY, "Try changing the equations in ICC, or"
+            strings = (
+                "An uphill search direction was found.",
+                "Try changing the equations in ICC, or",
+                "adding new variables to IXC.",
             )
-            process_output.ocmmnt(constants.IOTTY, "adding new variables to IXC.")
         elif ifail == 5:
-            process_output.ocmmnt(
-                constants.NOUT, "The quadratic programming technique was unable to"
-            )
-            process_output.ocmmnt(constants.NOUT, "find a feasible point.")
-            process_output.oblnkl(constants.NOUT)
-            process_output.ocmmnt(
-                constants.NOUT, "Try changing or adding variables to IXC, or modify"
-            )
-            process_output.ocmmnt(
-                constants.NOUT,
+            strings = (
+                "The quadratic programming technique was unable to",
+                "find a feasible point.",
+                "",
+                "Try changing or adding variables to IXC, or modify",
                 "their initial values (especially if only 1 optimisation",
+                "iteration was performed).",
             )
-            process_output.ocmmnt(constants.NOUT, "iteration was performed).")
 
-            process_output.ocmmnt(
-                constants.IOTTY, "The quadratic programming technique was unable to"
-            )
-            process_output.ocmmnt(constants.IOTTY, "find a feasible point.")
-            process_output.oblnkl(constants.IOTTY)
-            process_output.ocmmnt(
-                constants.IOTTY, "Try changing or adding variables to IXC, or modify"
-            )
-            process_output.ocmmnt(
-                constants.IOTTY,
-                "their initial values (especially if only 1 optimisation",
-            )
-            process_output.ocmmnt(constants.IOTTY, "iteration was performed).")
         elif ifail == 6:
-            process_output.ocmmnt(
-                constants.NOUT, "The quadratic programming technique was restricted"
-            )
-            process_output.ocmmnt(
-                constants.NOUT, "by an artificial bound, or failed due to a singular"
+            strings = (
+                "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.",
             )
-            process_output.ocmmnt(constants.NOUT, "matrix.")
-            process_output.ocmmnt(
-                constants.NOUT, "Try changing the equations in ICC, or"
-            )
-            process_output.ocmmnt(constants.NOUT, "adding new variables to IXC.")
 
-            process_output.ocmmnt(
-                constants.IOTTY, "The quadratic programming technique was restricted"
-            )
-            process_output.ocmmnt(
-                constants.IOTTY, "by an artificial bound, or failed due to a singular"
-            )
-            process_output.ocmmnt(constants.IOTTY, "matrix.")
-            process_output.ocmmnt(
-                constants.IOTTY, "Try changing the equations in ICC, or"
-            )
-            process_output.ocmmnt(constants.IOTTY, "adding new variables to IXC.")
+        strings = "\n".join(strings)
+        process_output.ocmmnt(constants.NOUT, strings)
+        print(strings)
 
     def scan_1d(self):
         """Run a 1-D scan."""
@@ -931,20 +817,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}%")
 
@@ -1005,95 +889,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} : "
@@ -1107,18 +955,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,
@@ -1128,9 +970,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}, "
@@ -1141,27 +981,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 bdb2fb9eead309d4374f38440f70c4c63042646d 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 | 216 ++++++++++++++++++++-----------------------
 1 file changed, 102 insertions(+), 114 deletions(-)

diff --git a/process/core/scan.py b/process/core/scan.py
index 276ef1cbdc..c5dd955afd 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -329,20 +329,11 @@ def post_optimise(self, ifail: int):
             print()
         else:
             # Solution found
-            if self.solver != "fsolve":
-                process_output.ocmmnt(
-                    constants.NOUT, "and found a feasible set of parameters."
-                )
-                process_output.oheadr(
-                    constants.IOTTY, "PROCESS found a feasible solution"
-                )
-            else:
-                process_output.ocmmnt(
-                    constants.NOUT, "and found a consistent set of parameters."
-                )
-                process_output.oheadr(
-                    constants.IOTTY, "PROCESS found a consistent solution"
-                )
+            descr = "consistent" if self.solver == "fsolve" else "feasible"
+            process_output.ocmmnt(
+                constants.NOUT, f"and found a {descr} set of parameters."
+            )
+            process_output.oheadr(constants.IOTTY, f"PROCESS found a {descr} solution")
             process_output.oblnkl(constants.NOUT)
             process_output.ovarre(constants.NOUT, "Error flag", "(ifail)", ifail)
 
@@ -359,28 +350,22 @@ def post_optimise(self, ifail: int):
                     f"High final constraint residues. {self.data.numerics.sqsumsq=}"
                 )
 
-        process_output.ovarre(
-            constants.NOUT,
-            "Number of iteration variables",
-            "(nvar)",
-            self.data.numerics.nvar,
-        )
-        process_output.ovarre(
-            constants.NOUT,
-            "Number of constraints (total)",
-            "(neqns+nineqns)",
-            self.data.numerics.neqns + self.data.numerics.nineqns,
-        )
-        process_output.ovarre(
-            constants.NOUT,
-            "Optimisation switch",
-            "(ioptimz)",
-            self.data.numerics.ioptimz,
-        )
+        for d, var, v in (
+            ("Number of iteration variables", "(nvar)", self.data.numerics.nvar),
+            (
+                "Number of constraints (total)",
+                "(neqns+nineqns)",
+                self.data.numerics.neqns + self.data.numerics.nineqns,
+            ),
+            ("Optimisation switch", "(ioptimz)", self.data.numerics.ioptimz),
+        ):
+            process_output.ovarin(constants.NOUT, d, var, v)
+
         process_output.ocmmnt(
             constants.NOUT,
             f"     {PROCESSRunMode(self.data.numerics.ioptimz).description}",
         )
+
         # Objective function output: none for fsolve
         if self.solver != "fsolve":
             process_output.ovarre(
@@ -390,69 +375,66 @@ def post_optimise(self, ifail: int):
                 self.data.numerics.minmax,
             )
 
-            objf_name = f'"{FiguresOfMerit(abs(self.data.numerics.minmax)).description}"'
-
-            self.data.numerics.objf_name = objf_name
-
-            process_output.ovarre(
-                constants.NOUT,
-                "Objective function name",
-                "(objf_name)",
-                self.data.numerics.objf_name,
-            )
-            process_output.ovarre(
-                constants.NOUT,
-                "Normalised objective function",
-                "(norm_objf)",
-                self.data.numerics.norm_objf,
-                "OP ",
+            self.data.numerics.objf_name = (
+                f'"{FiguresOfMerit(abs(self.data.numerics.minmax)).description}"'
             )
 
-        process_output.ovarre(
-            constants.NOUT,
-            "Square root of the sum of squares of the constraint residuals",
-            "(sqsumsq)",
-            self.data.numerics.sqsumsq,
-            "OP ",
-        )
-        if self.solver != "fsolve":
-            process_output.ovarre(
-                constants.NOUT,
-                "VMCON convergence parameter",
-                "(convergence_parameter)",
-                self.data.globals.convergence_parameter,
-                "OP ",
-            )
-            process_output.ovarre(
-                constants.NOUT,
-                "Number of optimising solver iterations",
-                "(nviter)",
-                self.data.numerics.nviter,
-                "OP ",
-            )
+            for d, var, v, o in (
+                (
+                    "Objective function name",
+                    "(objf_name)",
+                    self.data.numerics.objf_name,
+                    "",
+                ),
+                (
+                    "Normalised objective function",
+                    "(norm_objf)",
+                    self.data.numerics.norm_objf,
+                    "OP ",
+                ),
+                (
+                    "VMCON convergence parameter",
+                    "(convergence_parameter)",
+                    self.data.globals.convergence_parameter,
+                    "OP ",
+                ),
+                (
+                    "Number of optimising solver iterations",
+                    "(nviter)",
+                    self.data.numerics.nviter,
+                    "OP ",
+                ),
+                (
+                    "Square root of the sum of squares of the constraint residuals",
+                    "(sqsumsq)",
+                    self.data.numerics.sqsumsq,
+                    "OP ",
+                ),
+            ):
+                process_output.ovarre(constants.NOUT, d, var, v, o)
+
         process_output.oblnkl(constants.NOUT)
 
         if self.solver == "fsolve":
-            if ifail == 1:
-                msg = "PROCESS has solved using fsolve."
-            else:
-                msg = "PROCESS failed to solve using fsolve."
             process_output.write(
                 constants.NOUT,
-                f"{msg}\n",
+                "PROCESS has solved using fsolve.\n"
+                if ifail == 1
+                else "PROCESS failed to solve using fsolve.\n",
             )
         else:
-            if ifail == 1:
-                string1 = "PROCESS has successfully optimised"
-            else:
-                string1 = "PROCESS has failed to optimise"
-
-            string2 = "minimise" if self.data.numerics.minmax > 0 else "maximise"
-
             process_output.write(
                 constants.NOUT,
-                f"{string1} the optimisation parameters to {string2} "
-                f"the objective function: {objf_name}\n",
+                (
+                    (
+                        "PROCESS has successfully optimised"
+                        if ifail == 1
+                        else "PROCESS has failed to optimise"
+                    )
+                    + " the optimisation parameters to"
+                    + ("minimise" if self.data.numerics.minmax > 0 else "maximise")
+                    + f" the objective function: {self.data.numerics.objf_name}\n"
+                ),
             )
 
         written_warning = False
@@ -510,10 +492,10 @@ def post_optimise(self, ifail: int):
                 self.data.numerics.xcs[i],
             )
 
-            if self.data.numerics.boundu[i] == self.data.numerics.boundl[i]:
-                xnorm = 1.0
-            else:
-                xnorm = min(
+            xnorm = (
+                1.0
+                if self.data.numerics.boundu[i] == self.data.numerics.boundl[i]
+                else min(
                     max(
                         (
                             self.data.numerics.xcm[i]
@@ -527,6 +509,7 @@ def post_optimise(self, ifail: int):
                     ),
                     1.0,
                 )
+            )
 
             process_output.ovarre(
                 constants.MFILE,
@@ -821,13 +804,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}%")
@@ -861,7 +844,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(
@@ -894,36 +894,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
 
@@ -964,11 +953,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 12019505a662a73c50fad2afb9cca9a8d371f8b4 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 | 53 +++++++++++++--------------
 1 file changed, 26 insertions(+), 27 deletions(-)

diff --git a/process/core/solver/solver_handler.py b/process/core/solver/solver_handler.py
index c7c14c3af5..9a117bd864 100644
--- a/process/core/solver/solver_handler.py
+++ b/process/core/solver/solver_handler.py
@@ -1,5 +1,7 @@
 """Module containing solver handler routines"""
 
+from contextlib import contextmanager
+
 from process.core.solver.evaluators import Evaluators
 from process.core.solver.iteration_variables import (
     load_iteration_variables,
@@ -38,14 +40,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
@@ -56,30 +53,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 != 1:
-                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 != 1): supply ifail value
-                # to next attempt
-                self.data.numerics.epsfcn /= 10  # reset value
-
-            if ifail != 1:
-                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
+                with epsfcn_context(self.data.numerics):
+                    self.solver.solve()
 
             # If VMCON has exited with error code 5 try another run using a multiple
             # of the identity matrix as input for the Hessian b(n,n)
@@ -106,3 +91,17 @@ def output(self):
         # than required, size
         self.data.numerics.xcm[: self.solver.x.shape[0]] = self.solver.x
         self.data.numerics.rcm[: self.solver.conf.shape[0]] = self.solver.conf
+
+
+@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 568d81eb3e3eb37ecc4a3c8115c5e593de7a3463 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 | 5 +++--
 process/core/scan.py          | 1 +
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/process/core/io/plot/scans.py b/process/core/io/plot/scans.py
index 43aa121f11..3e73cdc16d 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
 
@@ -532,9 +533,9 @@ def plot_1d_scan(
         # This uses exclusively the last output_name2 defined in an earlier loop ignoring
         # all other output_name2s...
         if output_name == "plasma_current_MA":
-            extra_str = f"plasma_current{f'_vs_{output_name2}' if len(output_names2) > 0 else ''}"  # ruff:ignore[line-too-long]
+            extra_str = f"plasma_current{f'_vs_{output_name2}' if len(output_names2) > 0 else ''}"  # noqa: E501
         elif stack_plots and output_names[-1] == output_name:
-            extra_str = f"{output_name}{f'_vs_{output_name2}' if len(output_names2) > 0 else '_vs_'.join(output_names)}"  # ruff:ignore[line-too-long]
+            extra_str = f"{output_name}{f'_vs_{output_name2}' if len(output_names2) > 0 else '_vs_'.join(output_names)}"  # noqa: E501
         else:
             extra_str = (
                 f"{output_name}{f'_vs_{output_name2}' if len(output_names2) > 0 else ''}"
diff --git a/process/core/scan.py b/process/core/scan.py
index c5dd955afd..5f2a8ed165 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -983,3 +983,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 5c7190a9817df19cd2a37b108f4bb83cd34810eb 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 | 336 +++++++++++++++++++++----------------------
 1 file changed, 160 insertions(+), 176 deletions(-)

diff --git a/process/core/scan.py b/process/core/scan.py
index 5f2a8ed165..f66b20f2c9 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -227,25 +227,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.
@@ -260,38 +258,87 @@ 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
+
+        if self.data.scan.scan_dim == 2:
+            self.scan_2d()
+        else:
+            self.scan_1d()
+
+    def write_outputs(self):
+        write_output_files(
+            models=self.models,
+            data=self.data,
+            ifail=self._ifail,
+            runtime=self._finish_time - self._start_time,
+        )
+        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)
-            return
+            logging_model_handler.clear_logs()
 
-        if self.data.scan.isweep > IPNSCNS:
-            raise ProcessValueError(
-                "Illegal value of isweep",
-                isweep=self.data.scan.isweep,
-                IPNSCNS=IPNSCNS,
-            )
+    def scan_2d(self):
+        iscan = 1
 
-        if self.data.scan.scan_dim == 2:
-            self.scan_2d()
-        else:
-            self.scan_1d()
+        # 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 doopt(self):
-        """Run the optimiser or solver."""
-        ifail = self.solver_handler.run()
-        self.post_optimise(ifail)
 
-        return ifail
+def scan_summary(scan):
+
+    def write_outputs(self):
+        self.post_optimise(self._ifail)
 
     def post_optimise(self, ifail: int):
         """Called after calling the optimising equation solver from Python.
@@ -484,14 +531,6 @@ def post_optimise(self, ifail: int):
                     f" {bounds[i] * self.data.numerics.scafc[i]}",
                 )
 
-            # Write optimisation parameters to mfile
-            process_output.ovarre(
-                constants.MFILE,
-                self.data.numerics.lablxc[self.data.numerics.ixc[i] - 1],
-                f"(itvar{i + 1:03d})",
-                self.data.numerics.xcs[i],
-            )
-
             xnorm = (
                 1.0
                 if self.data.numerics.boundu[i] == self.data.numerics.boundl[i]
@@ -511,32 +550,33 @@ def post_optimise(self, ifail: int):
                 )
             )
 
-            process_output.ovarre(
-                constants.MFILE,
-                f"{name} (final value/initial value)",
-                f"(xcm{i + 1:03d})",
-                self.data.numerics.xcm[i],
-            )
-            process_output.ovarre(
-                constants.MFILE,
-                f"{name} (range normalised)",
-                f"(nitvar{i + 1:03d})",
-                xnorm,
-            )
-            process_output.ovarre(
-                constants.MFILE,
-                f"{name} (upper bound)",
-                f"(boundu{i + 1:03d})",
-                self.data.numerics.itv_scaled_upper_bounds[i]
-                * self.data.numerics.scafc[i],
-            )
-            process_output.ovarre(
-                constants.MFILE,
-                f"{name} (lower bound)",
-                f"(boundl{i + 1:03d})",
-                self.data.numerics.itv_scaled_lower_bounds[i]
-                * self.data.numerics.scafc[i],
-            )
+            # Write optimisation parameters to mfile
+            for d, var, v in (
+                (
+                    self.data.numerics.lablxc[self.data.numerics.ixc[i] - 1],
+                    f"(itvar{i + 1:03d})",
+                    self.data.numerics.xcs[i],
+                ),
+                (
+                    f"{name} (final value/initial value)",
+                    f"(xcm{i + 1:03d})",
+                    self.data.numerics.xcm[i],
+                ),
+                (f"{name} (range normalised)", f"(nitvar{i + 1:03d})", xnorm),
+                (
+                    f"{name} (upper bound)",
+                    f"(boundu{i + 1:03d})",
+                    self.data.numerics.itv_scaled_upper_bounds[i]
+                    * self.data.numerics.scafc[i],
+                ),
+                (
+                    f"{name} (lower bound)",
+                    f"(boundl{i + 1:03d})",
+                    self.data.numerics.itv_scaled_lower_bounds[i]
+                    * self.data.numerics.scafc[i],
+                ),
+            ):
+                process_output.ovarre(constants.MFILE, d, var, v)
 
         # Write optimisation parameter headings to output file
         process_output.osubhd(
@@ -572,32 +612,30 @@ def post_optimise(self, ifail: int):
                 f"{err[i]} {lab[i]}",
                 con1[i],
             ])
-            process_output.ovarre(
-                constants.MFILE,
-                f"{name:<33} normalised residue",
-                f"(eq_con{self.data.numerics.icc[i]:03d})",
-                con1[i],
-            )
 
-            process_output.ovarre(
-                constants.MFILE,
-                f"{name:<33} residual",
-                f"(res_eq_con{self.data.numerics.icc[i]:03d})",
-                err[i],
-            )
-            process_output.ovarre(
-                constants.MFILE,
-                f"{name} constraint value",
-                f"(val_eq_con{self.data.numerics.icc[i]:03d})",
-                con2[i],
-            )
-
-            process_output.ovarre(
-                constants.MFILE,
-                f"{name} units",
-                f"(eq_units_con{self.data.numerics.icc[i]:03d})",
-                f"'{lab[i]}'",
-            )
+            for d, var, v in (
+                (
+                    f"{name:<33} normalised residue",
+                    f"(eq_con{self.data.numerics.icc[i]:03d})",
+                    con1[i],
+                ),
+                (
+                    f"{name:<33} residual",
+                    f"(res_eq_con{self.data.numerics.icc[i]:03d})",
+                    err[i],
+                ),
+                (
+                    f"{name} constraint value",
+                    f"(val_eq_con{self.data.numerics.icc[i]:03d})",
+                    con2[i],
+                ),
+                (
+                    f"{name} units",
+                    f"(eq_units_con{self.data.numerics.icc[i]:03d})",
+                    f"'{lab[i]}'",
+                ),
+            ):
+                process_output.ovarre(constants.MFILE, d, var, v)
 
         # Write equality constraints to output file
         process_output.write(
@@ -648,39 +686,35 @@ def post_optimise(self, ifail: int):
                     f"{constraint.residual} {constraint.units}",
                     f"{constraint.normalised_residual}",
                 ])
-                process_output.ovarre(
-                    constants.MFILE,
-                    f"{name} normalised residue",
-                    f"(ineq_con{self.data.numerics.icc[i]:03d})",
-                    -constraint.normalised_residual,
-                )
-                process_output.ovarre(
-                    constants.MFILE,
-                    f"{name} physical value",
-                    f"(ineq_value_con{self.data.numerics.icc[i]:03d})",
-                    constraint.constraint_value,
-                )
 
-                process_output.ovarre(
-                    constants.MFILE,
-                    f"{name} symbol",
-                    f"(ineq_symbol_con{self.data.numerics.icc[i]:03d})",
-                    f"'{constraint.symbol}'",
-                )
-
-                process_output.ovarre(
-                    constants.MFILE,
-                    f"{name} units",
-                    f"(ineq_units_con{self.data.numerics.icc[i]:03d})",
-                    f"'{constraint.units}'",
-                )
-
-                process_output.ovarre(
-                    constants.MFILE,
-                    f"{name} physical bound",
-                    f"(ineq_bound_con{self.data.numerics.icc[i]:03d})",
-                    constraint.constraint_bound,
-                )
+                for d, var, v in (
+                    (
+                        "normalised residue",
+                        f"(ineq_con{self.data.numerics.icc[i]:03d})",
+                        -constraint.normalised_residual,
+                    ),
+                    (
+                        "physical value",
+                        f"(ineq_value_con{self.data.numerics.icc[i]:03d})",
+                        constraint.constraint_value,
+                    ),
+                    (
+                        "symbol",
+                        f"(ineq_symbol_con{self.data.numerics.icc[i]:03d})",
+                        f"'{constraint.symbol}'",
+                    ),
+                    (
+                        "units",
+                        f"(ineq_units_con{self.data.numerics.icc[i]:03d})",
+                        f"'{constraint.units}'",
+                    ),
+                    (
+                        "physical bound",
+                        f"(ineq_bound_con{self.data.numerics.icc[i]:03d})",
+                        constraint.constraint_bound,
+                    ),
+                ):
+                    process_output.ovarre(constants.MFILE, f"{name} {d}", var, v)
 
             process_output.write(
                 constants.NOUT,
@@ -771,21 +805,6 @@ def scan_1d(self):
         # 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()
-            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()
-
         # outvar now contains results
         self.scan_1d_write_plot(self.data.scan)
         print("Scan Convergence Summary \n")
@@ -802,7 +821,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"
             )
@@ -815,37 +834,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 (
@@ -859,9 +847,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(
@@ -890,8 +876,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"
@@ -927,8 +913,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)
 
@@ -955,8 +940,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 ab3402c4726c9887c295b0917fd2f33e13dd0784 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 f66b20f2c9..8fc3f0d5ad 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -834,6 +834,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 257f97a013e1525b83509dc2988ac5dd3b2c89be 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/caller.py                     | 100 ++-
 process/core/final.py                      | 105 ---
 process/core/input.py                      |  13 +-
 process/core/output.py                     | 153 ----
 process/core/scan.py                       | 846 ++++-----------------
 process/core/solver/constraints.py         | 123 ++-
 process/core/solver/iteration_variables.py |   5 +-
 process/core/solver/solver.py              |  60 ++
 process/core/solver/solver_handler.py      | 123 +++
 process/data_structure/global_variables.py |  18 +-
 process/data_structure/scan_variables.py   |  70 +-
 process/main.py                            | 145 +++-
 12 files changed, 756 insertions(+), 1005 deletions(-)

diff --git a/process/core/caller.py b/process/core/caller.py
index e618e51262..1634777e53 100644
--- a/process/core/caller.py
+++ b/process/core/caller.py
@@ -9,13 +9,14 @@
 from tabulate import tabulate
 
 from process.core import constants
-from process.core.final import finalise
+from process.core import process_output as po
 from process.core.io.mfile import MFile
 from process.core.process_output import OutputFileManager, ovarre
 from process.core.solver import constraints
 from process.core.solver.iteration_variables import set_scaled_iteration_variable
 from process.core.solver.objectives import objective_function
 from process.data_structure.blanket_variables import BlktModelTypes
+from process.data_structure.numerics import PROCESSRunMode
 from process.models.tfcoil.base import TFConductorModel
 from process.models.tfcoil.superconducting import SuperconductingTFTurnType
 
@@ -396,6 +397,103 @@ def _call_models_once(self, xc: np.ndarray):
         # FISPACT and LOCA model (not used)- removed
 
 
+def finalise(models, data, ifail: int, non_idempotent_msg: str | None = None):
+    """Routine to print out the final point in the scan.
+
+    Writes to OUT.DAT and MFILE.DAT.
+
+    Parameters
+    ----------
+    models : process.main.Models
+        physics and engineering model objects
+    data: DataStructure
+        data structure object to provide data to evaluate the constraints
+    ifail : int
+        error flag
+    non_idempotent_msg : None | str, optional
+        warning about non-idempotent variables, defaults to None
+    """
+    if ifail == 1:
+        po.oheadr(constants.NOUT, "Final Feasible Point")
+    else:
+        po.oheadr(constants.NOUT, "Final UNFEASIBLE Point")
+
+    # Output relevant to no optimisation
+    if data.numerics.ioptimz == PROCESSRunMode.EVALUATION:
+        output_evaluation(data)
+
+    # Print non-idempotence warning to OUT.DAT only
+    if non_idempotent_msg:
+        po.oheadr(constants.NOUT, "NON-IDEMPOTENT VARIABLES")
+        po.ocmmnt(constants.NOUT, non_idempotent_msg)
+
+    # Write output to OUT.DAT and MFILE.DAT
+    models.write(data, constants.NOUT)
+
+
+def output_evaluation(data):
+    """Write output for an evaluation run of PROCESS
+
+    Parameters
+    ----------
+    data: DataStructure
+        data structure object to provide data to evaluate the constraints
+    """
+    po.oheadr(constants.NOUT, "Numerics")
+    po.ocmmnt(constants.NOUT, "PROCESS has performed an evaluation run.")
+    po.oblnkl(constants.NOUT)
+
+    # Evaluate objective function
+    norm_objf = objective_function(data.numerics.minmax, data)
+    po.ovarre(constants.MFILE, "Normalised objective function", "(norm_objf)", norm_objf)
+
+    # Print the residuals of the constraint equations
+
+    residual_error, value, residual, symbols, units = constraints.constraint_eqns(
+        data.numerics.neqns + data.numerics.nineqns, -1, data
+    )
+
+    labels = [
+        data.numerics.lablcc[j - 1]
+        for j in data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns]
+    ]
+
+    def _fmt(a, units):
+        return [f"{c} {u}" for c, u in zip(a, units, strict=False)]
+
+    po.write(
+        constants.NOUT,
+        tabulate(
+            {
+                "Constraint Name": labels,
+                "Constraint Type": symbols,
+                "Physical constraint": _fmt(value, units),
+                "Constraint residual": _fmt(residual, units),
+                "Normalised residual": residual_error,
+            },
+            headers="keys",
+        ),
+    )
+
+    for i in range(data.numerics.neqns):
+        constraint_id = data.numerics.icc[i]
+        po.ovarre(
+            constants.MFILE,
+            f"{labels[i]} normalised residue",
+            f"(eq_con{constraint_id:03d})",
+            residual_error[i],
+        )
+
+    for i in range(data.numerics.nineqns):
+        constraint_id = data.numerics.icc[data.numerics.neqns + i]
+        po.ovarre(
+            constants.MFILE,
+            f"{labels[data.numerics.neqns + i]}",
+            f"(ineq_con{constraint_id:03d})",
+            residual_error[data.numerics.neqns + i],
+        )
+
+
 def write_output_files(
     models: Models, data: DataStructure, ifail: int, *, runtime: float | None = None
 ):
diff --git a/process/core/final.py b/process/core/final.py
index d5da7593bc..e69de29bb2 100644
--- a/process/core/final.py
+++ b/process/core/final.py
@@ -1,105 +0,0 @@
-"""Final output at the end of a scan."""
-
-from tabulate import tabulate
-
-from process.core import constants
-from process.core import output as op
-from process.core import process_output as po
-from process.core.solver import constraints
-from process.core.solver.objectives import objective_function
-from process.data_structure.numerics import PROCESSRunMode
-
-
-def finalise(models, data, ifail: int, non_idempotent_msg: str | None = None):
-    """Routine to print out the final point in the scan.
-
-    Writes to OUT.DAT and MFILE.DAT.
-
-    Parameters
-    ----------
-    models : process.main.Models
-        physics and engineering model objects
-    data: DataStructure
-        data structure object to provide data to evaluate the constraints
-    ifail : int
-        error flag
-    non_idempotent_msg : None | str, optional
-        warning about non-idempotent variables, defaults to None
-    """
-    if ifail == 1:
-        po.oheadr(constants.NOUT, "Final Feasible Point")
-    else:
-        po.oheadr(constants.NOUT, "Final UNFEASIBLE Point")
-
-    # Output relevant to no optimisation
-    if data.numerics.ioptimz == PROCESSRunMode.EVALUATION:
-        output_evaluation(data)
-
-    # Print non-idempotence warning to OUT.DAT only
-    if non_idempotent_msg:
-        po.oheadr(constants.NOUT, "NON-IDEMPOTENT VARIABLES")
-        po.ocmmnt(constants.NOUT, non_idempotent_msg)
-
-    # Write output to OUT.DAT and MFILE.DAT
-    op.write(models, data, constants.NOUT)
-
-
-def output_evaluation(data):
-    """Write output for an evaluation run of PROCESS
-
-    Parameters
-    ----------
-    data: DataStructure
-        data structure object to provide data to evaluate the constraints
-    """
-    po.oheadr(constants.NOUT, "Numerics")
-    po.ocmmnt(constants.NOUT, "PROCESS has performed an evaluation run.")
-    po.oblnkl(constants.NOUT)
-
-    # Evaluate objective function
-    norm_objf = objective_function(data.numerics.minmax, data)
-    po.ovarre(constants.MFILE, "Normalised objective function", "(norm_objf)", norm_objf)
-
-    # Print the residuals of the constraint equations
-
-    residual_error, value, residual, symbols, units = constraints.constraint_eqns(
-        data.numerics.neqns + data.numerics.nineqns, -1, data
-    )
-
-    labels = [
-        data.numerics.lablcc[j]
-        for j in [
-            i - 1
-            for i in data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns]
-        ]
-    ]
-    physical_constraint = [f"{c} {u}" for c, u in zip(value, units, strict=False)]
-    physical_residual = [f"{c} {u}" for c, u in zip(residual, units, strict=False)]
-
-    table_data = {
-        "Constraint Name": labels,
-        "Constraint Type": symbols,
-        "Physical constraint": physical_constraint,
-        "Constraint residual": physical_residual,
-        "Normalised residual": residual_error,
-    }
-
-    po.write(constants.NOUT, tabulate(table_data, headers="keys"))
-
-    for i in range(data.numerics.neqns):
-        constraint_id = data.numerics.icc[i]
-        po.ovarre(
-            constants.MFILE,
-            f"{labels[i]} normalised residue",
-            f"(eq_con{constraint_id:03d})",
-            residual_error[i],
-        )
-
-    for i in range(data.numerics.nineqns):
-        constraint_id = data.numerics.icc[data.numerics.neqns + i]
-        po.ovarre(
-            constants.MFILE,
-            f"{labels[data.numerics.neqns + i]}",
-            f"(ineq_con{constraint_id:03d})",
-            residual_error[data.numerics.neqns + i],
-        )
diff --git a/process/core/input.py b/process/core/input.py
index 1d7c132b6a..57a0f2757b 100644
--- a/process/core/input.py
+++ b/process/core/input.py
@@ -1138,24 +1138,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/output.py b/process/core/output.py
index 770c1f486d..e69de29bb2 100644
--- a/process/core/output.py
+++ b/process/core/output.py
@@ -1,153 +0,0 @@
-"""Module containing routine to write the results to the main output file (OUT.DAT)"""
-
-from process.core.log import logging_model_handler
-from process.data_structure.blanket_variables import BlktModelTypes
-from process.models.tfcoil.base import TFConductorModel
-from process.models.tfcoil.superconducting import (
-    SuperconductingTFTurnType,
-)
-
-
-def write(models, data, _outfile):
-    """Write the results to the main output file (OUT.DAT).
-
-    Write the program results to a file, in a tidy format.
-
-    Parameters
-    ----------
-    models : process.main.Models
-        physics and engineering model objects
-    _outfile : int
-        Fortran output unit identifier
-
-    Raises
-    ------
-    ValueError
-        If unsupported superconducting TF turn type is used
-    """
-    # ensure we are capturing warnings that occur in the 'output' stage
-    # as these are warnings that occur at our solution point.
-    # So we clear existing warnings
-    logging_model_handler.start_capturing()
-    logging_model_handler.clear_logs()
-
-    # Call stellarator output routine instead if relevant
-    if data.stellarator.istell != 0:
-        models.stellarator.output()
-        return
-
-    #  Call IFE output routine instead if relevant
-    if data.ife.ife != 0:
-        models.ife.output()
-        return
-
-    # Costs model
-    # Cost switch values
-    # No.  |  model
-    # ---- | ------
-    # 0    |  1990 costs model
-    # 1    |  2015 Kovari model
-    # 2    |  Custom model
-    models.costs.output()
-
-    # Availability model
-    models.availability.output()
-
-    # Physics model
-    models.physics.output()
-
-    # Detailed physics, currently only done at final point as values are not used
-    # by any other functions
-    models.physics_detailed.output()
-
-    # TODO what is this? Not in caller.py?
-    models.current_drive.output()
-
-    # Pulsed reactor model
-    models.pulse.output()
-
-    models.divertor.output()
-
-    # Machine Build Model
-    models.build.output()
-
-    # Cryostat build
-    models.cryostat.output()
-
-    # Toroidal field coil copper model
-    if data.tfcoil.i_tf_sup == TFConductorModel.WATER_COOLED_COPPER:
-        models.copper_tf_coil.output()
-
-    # Toroidal field coil superconductor model
-    if data.tfcoil.i_tf_sup == TFConductorModel.SUPERCONDUCTING:
-        tf_turn_type = SuperconductingTFTurnType(
-            data.superconducting_tfcoil.i_tf_turn_type
-        )
-        if tf_turn_type == SuperconductingTFTurnType.CABLE_IN_CONDUIT:
-            models.cicc_sctfcoil.output()
-        elif tf_turn_type == SuperconductingTFTurnType.CROSS_CONDUCTOR:
-            models.croco_sctfcoil.output()
-        else:
-            raise ValueError(
-                "Unsupported superconducting TF turn type: "
-                f"{data.superconducting_tfcoil.i_tf_turn_type}"
-            )
-
-    # Toroidal field coil aluminium model
-    if data.tfcoil.i_tf_sup == TFConductorModel.HELIUM_COOLED_ALUMINIUM:
-        models.aluminium_tf_coil.output()
-
-    # Tight aspect ratio machine model
-    if (
-        data.physics.itart == 1
-        and data.tfcoil.i_tf_sup != TFConductorModel.SUPERCONDUCTING
-    ):
-        models.tfcoil.output()
-
-    # Poloidal field coil model
-    models.pfcoil.output()
-
-    # Structure Model
-    models.structure.output()
-
-    # Blanket model
-    # Blanket switch values
-    # No.  |  model
-    # ---- | ------
-    # 1    |  CCFE HCPB model
-    # 2    |  KIT HCPB model
-    # 3    |  CCFE HCPB model with Tritium Breeding Ratio calculation
-    # 4    |  KIT HCLL model
-    # 5    |  DCLL model
-
-    models.shield.output()
-    models.vacuum_vessel.output()
-
-    # First wall geometry
-    models.fw.output()
-
-    if data.fwbs.i_blanket_type == BlktModelTypes.CCFE_HCPB:
-        # CCFE HCPB model
-        models.ccfe_hcpb.output()
-
-    elif data.fwbs.i_blanket_type == BlktModelTypes.DCLL:
-        # DCLL model
-        models.dcll.output()
-
-    # FISPACT and LOCA model (not used)- removed
-
-    # Power model
-    models.power.output()
-
-    # Vacuum model
-    models.vacuum.output()
-
-    # Buildings model
-    models.buildings.output()
-
-    # Water usage in secondary cooling system
-    models.water_use.output()
-
-    # stop capturing warnings so that Outfile does not end up with
-    # a lot of non-model logs
-    logging_model_handler.stop_capturing()
diff --git a/process/core/scan.py b/process/core/scan.py
index 8fc3f0d5ad..0458b7f50a 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -20,11 +20,12 @@
 from process.core.solver import constraints
 from process.core.solver.solver_handler import SolverHandler
 from process.data_structure.numerics import FiguresOfMerit, PROCESSRunMode
-from process.data_structure.scan_variables import IPNSCNS, NOUTVARS, ScanData
+from process.data_structure.scan_variables import IPNSCNS, ScanData
 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__)
@@ -226,6 +227,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
 
@@ -239,15 +247,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
@@ -258,720 +264,198 @@ 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
-
-        if self.data.scan.scan_dim == 2:
-            self.scan_2d()
-        else:
-            self.scan_1d()
-
-    def write_outputs(self):
-        write_output_files(
-            models=self.models,
-            data=self.data,
-            ifail=self._ifail,
-            runtime=self._finish_time - self._start_time,
-        )
-        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,
-            )
+        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
 
-            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):
+        show_errors(constants.NOUT)
 
-    def write_outputs(self):
-        self.post_optimise(self._ifail)
+        logging_model_handler.clear_logs()
+        optimisation_output(data)
+        constraints.constraints_output(data, self.solver)
 
-    def post_optimise(self, ifail: int):
-        """Called after calling the optimising equation solver from Python.
+        return ScanRes(iscan, ifail, sh)
 
-        ifail   : input integer : error flag
+    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
+        ]
+        self.data.globals.vlabel = [s.fname for s in sv]
+        self.data.globals.xlabel = [s.data.description for s in sv]
 
-        Parameters
-        ----------
-        ifail: int :
+    def write_point_header(self, iscan):
+        self._set_v_x_label(iscan)
 
-        """
-        self.data.numerics.sqsumsq = (
-            sum(r**2 for r in self.data.numerics.rcm[: self.data.numerics.neqns]) ** 0.5
-        )
+        process_output.oblnkl(constants.NOUT)
+        process_output.oblnkl(constants.MFILE)
 
-        process_output.oheadr(constants.NOUT, "Numerics")
-        process_output.ocmmnt(
+        process_output.write(
             constants.NOUT,
-            f"PROCESS has performed a {'fsolve' if self.solver == 'fsolve' else 'VMCON'}"
-            " (optimisation) run.",
+            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)
+            ),
         )
-        if ifail != 1:
-            process_output.ovarre(constants.NOUT, "Error flag", "(ifail)", ifail)
-            process_output.oheadr(
-                constants.IOTTY, "PROCESS COULD NOT FIND A FEASIBLE SOLUTION"
-            )
-            print()
-
-            logger.critical("Solver returns with ifail /= 1. %s", ifail)
+        process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan)
 
-            # Error code handler for VMCON
-            if self.solver == "vmcon":
-                self.verror(ifail)
-            process_output.oblnkl(constants.NOUT)
-            print()
-        else:
-            # Solution found
-            descr = "consistent" if self.solver == "fsolve" else "feasible"
-            process_output.ocmmnt(
-                constants.NOUT, f"and found a {descr} set of parameters."
+        print(
+            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)
             )
-            process_output.oheadr(constants.IOTTY, f"PROCESS found a {descr} solution")
-            process_output.oblnkl(constants.NOUT)
-            process_output.ovarre(constants.NOUT, "Error flag", "(ifail)", ifail)
-
-            if self.data.numerics.sqsumsq >= 1.0e-2:
-                string = (
-                    "WARNING: Constraint residues are HIGH; consider re-running\n"
-                    "   with lower values of EPSVMC to confirm convergence...\n"
-                    "   (should be able to get down to about 1.0E-8 okay)\n"
-                )
-                process_output.ocmmnt(constants.NOUT, ("\n" + string))
-                print(string)
-
-                logger.warning(
-                    f"High final constraint residues. {self.data.numerics.sqsumsq=}"
-                )
-
-        for d, var, v in (
-            ("Number of iteration variables", "(nvar)", self.data.numerics.nvar),
-            (
-                "Number of constraints (total)",
-                "(neqns+nineqns)",
-                self.data.numerics.neqns + self.data.numerics.nineqns,
-            ),
-            ("Optimisation switch", "(ioptimz)", self.data.numerics.ioptimz),
-        ):
-            process_output.ovarin(constants.NOUT, d, var, v)
-
-        process_output.ocmmnt(
-            constants.NOUT,
-            f"     {PROCESSRunMode(self.data.numerics.ioptimz).description}",
         )
 
-        # Objective function output: none for fsolve
-        if self.solver != "fsolve":
-            process_output.ovarre(
-                constants.NOUT,
-                "Figure of merit switch",
-                "(minmax)",
-                self.data.numerics.minmax,
-            )
+    def scan_select(self, nsweep, sweep, iscan):
+        sv = ScanVariables(nsweep)
+        sv.set(self.data, sweep[iscan - 1])
+        return sv
 
-            self.data.numerics.objf_name = (
-                f'"{FiguresOfMerit(abs(self.data.numerics.minmax)).description}"'
-            )
+    def run(self):
+        """Call a solver over a range of values of one of the variables.
 
-            for d, var, v, o in (
-                (
-                    "Objective function name",
-                    "(objf_name)",
-                    self.data.numerics.objf_name,
-                    "",
-                ),
-                (
-                    "Normalised objective function",
-                    "(norm_objf)",
-                    self.data.numerics.norm_objf,
-                    "OP ",
-                ),
-                (
-                    "VMCON convergence parameter",
-                    "(convergence_parameter)",
-                    self.data.globals.convergence_parameter,
-                    "OP ",
-                ),
-                (
-                    "Number of optimising solver iterations",
-                    "(nviter)",
-                    self.data.numerics.nviter,
-                    "OP ",
-                ),
-                (
-                    "Square root of the sum of squares of the constraint residuals",
-                    "(sqsumsq)",
-                    self.data.numerics.sqsumsq,
-                    "OP ",
-                ),
+        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.ovarre(constants.NOUT, d, var, v, o)
-
-        process_output.oblnkl(constants.NOUT)
+                process_output.ovarin(constants.MFILE, d, n, v)
 
-        if self.solver == "fsolve":
-            process_output.write(
-                constants.NOUT,
-                "PROCESS has solved using fsolve.\n"
-                if ifail == 1
-                else "PROCESS failed to solve using fsolve.\n",
-            )
-        else:
-            process_output.write(
-                constants.NOUT,
-                (
-                    (
-                        "PROCESS has successfully optimised"
-                        if ifail == 1
-                        else "PROCESS has failed to optimise"
-                    )
-                    + " the optimisation parameters to"
-                    + ("minimise" if self.data.numerics.minmax > 0 else "maximise")
-                    + f" the objective function: {self.data.numerics.objf_name}\n"
-                ),
-            )
+        # 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
+        )
 
-        written_warning = False
+        if self.data.scan.nsweep is not None:
+            self.summary(scan_res)
 
-        # Output optimisation parameters
-        solution_vector_table = []
-        for i in range(self.data.numerics.nvar):
-            self.data.numerics.xcs[i] = (
-                self.data.numerics.xcm[i] * self.data.numerics.scafc[i]
-            )
+    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]
 
-            name = self.data.numerics.lablxc[self.data.numerics.ixc[i] - 1]
-            solution_vector_table.append([
-                name,
-                self.data.numerics.xcs[i],
-                self.data.numerics.xcm[i],
+        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,
             ])
 
-            xminn = 1.01 * self.data.numerics.itv_scaled_lower_bounds[i]
-            xmaxx = 0.99 * self.data.numerics.itv_scaled_upper_bounds[i]
-
-            # Write to output file if close to optimisation parameter bounds
-            if self.data.numerics.xcm[i] < xminn or self.data.numerics.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 = self.data.numerics.xcm[i] * self.data.numerics.scafc[i]
-
-                if self.data.numerics.xcm[i] < xminn:
-                    location, bound = "below", "lower"
-                    bounds = self.data.numerics.itv_scaled_lower_bounds
-                else:
-                    location, bound = "above", "upper"
-                    bounds = self.data.numerics.itv_scaled_upper_bounds
-                process_output.write(
-                    constants.NOUT,
-                    f"   {name:<30}= {xcval} is at or {location} its {bound} bound:"
-                    f" {bounds[i] * self.data.numerics.scafc[i]}",
-                )
-
-            xnorm = (
-                1.0
-                if self.data.numerics.boundu[i] == self.data.numerics.boundl[i]
-                else min(
-                    max(
-                        (
-                            self.data.numerics.xcm[i]
-                            - self.data.numerics.itv_scaled_lower_bounds[i]
-                        )
-                        / (
-                            self.data.numerics.itv_scaled_upper_bounds[i]
-                            - self.data.numerics.itv_scaled_lower_bounds[i]
-                        ),
-                        0.0,
-                    ),
-                    1.0,
-                )
-            )
-
-            # Write optimisation parameters to mfile
-            for d, var, v in (
-                (
-                    self.data.numerics.lablxc[self.data.numerics.ixc[i] - 1],
-                    f"(itvar{i + 1:03d})",
-                    self.data.numerics.xcs[i],
-                ),
-                (
-                    f"{name} (final value/initial value)",
-                    f"(xcm{i + 1:03d})",
-                    self.data.numerics.xcm[i],
-                ),
-                (f"{name} (range normalised)", f"(nitvar{i + 1:03d})", xnorm),
-                (
-                    f"{name} (upper bound)",
-                    f"(boundu{i + 1:03d})",
-                    self.data.numerics.itv_scaled_upper_bounds[i]
-                    * self.data.numerics.scafc[i],
-                ),
-                (
-                    f"{name} (lower bound)",
-                    f"(boundl{i + 1:03d})",
-                    self.data.numerics.itv_scaled_lower_bounds[i]
-                    * self.data.numerics.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",
-            ),
-        )
-
-        process_output.osubhd(
-            constants.NOUT,
-            "The following equality constraint residues should be close to zero:",
+        print(
+            tabulate(conv_list, headers=["Iscan", "Sweep Var", "Sweep Val", "Converged"])
         )
 
-        con1, con2, err, _, lab = constraints.constraint_eqns(
-            self.data.numerics.neqns + self.data.numerics.nineqns, -1, self.data
-        )
+        converged_percentage = converged_count / np.prod(self.data.scan.isweep) * 100
+        print(f"\nConvergence Percentage: {converged_percentage:.2f}%")
 
-        # Write equality constraints to mfile
-        equality_constraint_table = []
-        for i in range(self.data.numerics.neqns):
-            name = self.data.numerics.lablcc[self.data.numerics.icc[i] - 1]
-
-            equality_constraint_table.append([
-                name,
-                "=",
-                f"{con2[i]} {lab[i]}",
-                f"{err[i]} {lab[i]}",
-                con1[i],
-            ])
 
-            for d, var, v in (
-                (
-                    f"{name:<33} normalised residue",
-                    f"(eq_con{self.data.numerics.icc[i]:03d})",
-                    con1[i],
-                ),
-                (
-                    f"{name:<33} residual",
-                    f"(res_eq_con{self.data.numerics.icc[i]:03d})",
-                    err[i],
-                ),
-                (
-                    f"{name} constraint value",
-                    f"(val_eq_con{self.data.numerics.icc[i]:03d})",
-                    con2[i],
-                ),
-                (
-                    f"{name} units",
-                    f"(eq_units_con{self.data.numerics.icc[i]:03d})",
-                    f"'{lab[i]}'",
-                ),
-            ):
-                process_output.ovarre(constants.MFILE, d, var, v)
+def optimisation_output(data: DataStructure):
+    nums = data.numerics
 
-        # Write equality constraints to output file
-        process_output.write(
-            constants.NOUT,
-            tabulate(
-                equality_constraint_table,
-                headers=[
-                    "",
-                    "",
-                    "Physical constraint",
-                    "Constraint residue",
-                    "Normalised residue",
-                ],
-                numalign="left",
-            ),
-        )
+    written_warning = False
 
-        # Write inequality constraints
-        if self.data.numerics.nineqns > 0:
-            inequality_constraint_table = []
-            # Inequalities not necessarily satisfied when evaluating
-            process_output.osubhd(
-                constants.NOUT,
-                "Negative inequality constraint (normalised) residuals "
-                "indicate a constraint is satisfied.",
-            )
-            if self.solver == "fsolve":
-                process_output.osubhd(
-                    constants.NOUT,
-                    "This MFile was produced via an evaluation, not an optimisation, "
-                    "and so the constraints might be violated.",
-                )
+    # Output optimisation parameters
+    solution_vector_table = []
+    for i in range(nums.nvar):
+        nums.xcs[i] = nums.xcm[i] * nums.scafc[i]
 
-            for i in range(
-                self.data.numerics.neqns,
-                self.data.numerics.neqns + self.data.numerics.nineqns,
-            ):
-                name = self.data.numerics.lablcc[self.data.numerics.icc[i] - 1]
-                constraint = constraints.ConstraintManager.evaluate_constraint(
-                    int(self.data.numerics.icc[i]), self.data
-                )
+        name = nums.lablxc[nums.ixc[i] - 1]
+        solution_vector_table.append([name, nums.xcs[i], nums.xcm[i]])
 
-                inequality_constraint_table.append([
-                    name,
-                    f"{constraint.constraint_value} {constraint.units}",
-                    constraint.symbol,
-                    f"{constraint.constraint_bound} {constraint.units}",
-                    f"{constraint.residual} {constraint.units}",
-                    f"{constraint.normalised_residual}",
-                ])
+        xminn = 1.01 * nums.itv_scaled_lower_bounds[i]
+        xmaxx = 0.99 * nums.itv_scaled_upper_bounds[i]
 
-                for d, var, v in (
-                    (
-                        "normalised residue",
-                        f"(ineq_con{self.data.numerics.icc[i]:03d})",
-                        -constraint.normalised_residual,
-                    ),
-                    (
-                        "physical value",
-                        f"(ineq_value_con{self.data.numerics.icc[i]:03d})",
-                        constraint.constraint_value,
-                    ),
-                    (
-                        "symbol",
-                        f"(ineq_symbol_con{self.data.numerics.icc[i]:03d})",
-                        f"'{constraint.symbol}'",
-                    ),
-                    (
-                        "units",
-                        f"(ineq_units_con{self.data.numerics.icc[i]:03d})",
-                        f"'{constraint.units}'",
-                    ),
+        # 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,
                     (
-                        "physical bound",
-                        f"(ineq_bound_con{self.data.numerics.icc[i]:03d})",
-                        constraint.constraint_bound,
+                        "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"
                     ),
-                ):
-                    process_output.ovarre(constants.MFILE, f"{name} {d}", var, v)
+                )
 
+            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,
-                tabulate(
-                    inequality_constraint_table,
-                    headers=[
-                        "",
-                        "Physical constraint",
-                        "",
-                        "Physical constraint bound",
-                        "Constraint residue",
-                        "Normalised residue",
-                    ],
-                    numalign="left",
-                ),
+                f"   {name:<30}= {xcval} is at or {location} its {bound} bound:"
+                f" {bounds[i] * nums.scafc[i]}",
             )
 
-    @staticmethod
-    def verror(ifail: int):
-        """Routine to print out relevant messages in the case of an
-        unfeasible result from a VMCON (optimisation) run
-
-        ifail  : input integer : error flag
-        This routine prints out relevant messages in the case of
-        an unfeasible result from a VMCON (optimisation) run.
-
-        Parameters
-        ----------
-        ifail: int :
-
-        """
-        if ifail == -1:
-            strings = ("User-terminated execution of VMCON.",)
-        elif ifail == 0:
-            strings = (
-                "Improper input parameters to the VMCON routine.",
-                "PROCESS coding must be checked.",
-            )
-        elif ifail == 2:
-            strings = (
-                "The maximum number of calls has been reached without solution.",
-                (
-                    "The code may be stuck in a minimum in the residual space that is "
-                    "significantly above zero."
+        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,
                 ),
-                "",
-                "There is either no solution possible, or the code",
-                "is failing to escape from a deep local minimum.",
-                "Try changing the variables in IXC, or modify their initial values.",
-            )
-        elif ifail == 3:
-            strings = (
-                "The line search required the maximum of 10 calls.",
-                "A feasible solution may be difficult to achieve.",
-                "Try changing or adding variables to IXC.",
-            )
-        elif ifail == 4:
-            strings = (
-                "An uphill search direction was found.",
-                "Try changing the equations in ICC, or",
-                "adding new variables to IXC.",
-            )
-        elif ifail == 5:
-            strings = (
-                "The quadratic programming technique was unable to",
-                "find a feasible point.",
-                "",
-                "Try changing or adding variables to IXC, or modify",
-                "their initial values (especially if only 1 optimisation",
-                "iteration was performed).",
-            )
-
-        elif ifail == 6:
-            strings = (
-                "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.",
+                1.0,
             )
-
-        strings = "\n".join(strings)
-        process_output.ocmmnt(constants.NOUT, strings)
-        print(strings)
-
-    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
-        ]
-        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),
+        # 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.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}%")
-
-    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
-        self._set_v_x_label(iscan)
-
-        process_output.oblnkl(constants.NOUT)
-        process_output.oblnkl(constants.MFILE)
-
-        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]} ",
-        )
-        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]}"
-        )
-
-    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
-
-        self._set_v_x_label(iscan_1)
-        self._set_v_x_label(iscan_r)
-
-        process_output.oblnkl(constants.NOUT)
-        process_output.oblnkl(constants.MFILE)
-
-        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]} ",
-        )
-        process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan)
-
-        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]} "
-        )
-
-        return iscan_r
-
-    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 scan_select(self, nsweep, sweep, iscan) -> ScanVariables:
-        """Select a scan"""
-        sv = ScanVariables(nsweep)
-        sv.set(self.data, sweep[iscan - 1])
-        return sv
+            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/constraints.py b/process/core/solver/constraints.py
index e8def90ef6..dbef355604 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -6,8 +6,9 @@
 from typing import ClassVar, Literal
 
 import numpy as np
+from tabulate import tabulate
 
-from process.core import constants
+from process.core import constants, process_output
 from process.core.exceptions import ProcessError, ProcessValueError
 from process.core.model import DataStructure
 from process.data_structure.build_variables import TFCSRadialConfiguration
@@ -2010,3 +2011,123 @@ def constraint_eqns(m: int, ieqn: int, data: DataStructure):
         units.append(tmp_units)
 
     return np.array(cc), np.array(con), np.array(err), symbol, units
+
+
+def constraints_output(data: DataStructure, solver_name: str):
+    nums = data.numerics
+
+    process_output.osubhd(
+        constants.NOUT,
+        "The following equality constraint residues should be close to zero:",
+    )
+
+    con1, con2, err, _, lab = constraint_eqns(nums.neqns + nums.nineqns, -1, data)
+
+    # Write equality constraints to mfile
+    equality_constraint_table = []
+    for i in range(nums.neqns):
+        name = nums.lablcc[nums.icc[i] - 1]
+
+        equality_constraint_table.append([
+            name,
+            "=",
+            f"{con2[i]} {lab[i]}",
+            f"{err[i]} {lab[i]}",
+            con1[i],
+        ])
+
+        for d, var, v in (
+            (f"{name:<33} normalised residue", f"(eq_con{nums.icc[i]:03d})", con1[i]),
+            (f"{name:<33} residual", f"(res_eq_con{nums.icc[i]:03d})", err[i]),
+            (f"{name} constraint value", f"(val_eq_con{nums.icc[i]:03d})", con2[i]),
+            (f"{name} units", f"(eq_units_con{nums.icc[i]:03d})", f"'{lab[i]}'"),
+        ):
+            process_output.ovarre(constants.MFILE, d, var, v)
+
+    # Write equality constraints to output file
+    process_output.write(
+        constants.NOUT,
+        tabulate(
+            equality_constraint_table,
+            headers=[
+                "",
+                "",
+                "Physical constraint",
+                "Constraint residue",
+                "Normalised residue",
+            ],
+            numalign="left",
+        ),
+    )
+
+    # Write inequality constraints
+    if nums.nineqns > 0:
+        inequality_constraint_table = []
+        # Inequalities not necessarily satisfied when evaluating
+        process_output.osubhd(
+            constants.NOUT,
+            "Negative inequality constraint (normalised) residuals indicate a constraint is satisfied.",
+        )
+        if solver_name == "fsolve":
+            process_output.osubhd(
+                constants.NOUT,
+                "This MFile was produced via an evaluation, not an optimisation, and so the constraints "
+                "might be violated.",
+            )
+
+        for i in range(
+            nums.neqns,
+            nums.neqns + nums.nineqns,
+        ):
+            name = nums.lablcc[nums.icc[i] - 1]
+            constraint = ConstraintManager.evaluate_constraint(int(nums.icc[i]), data)
+
+            inequality_constraint_table.append([
+                name,
+                f"{constraint.constraint_value} {constraint.units}",
+                constraint.symbol,
+                f"{constraint.constraint_bound} {constraint.units}",
+                f"{constraint.residual} {constraint.units}",
+                f"{constraint.normalised_residual}",
+            ])
+
+            for d, var, v in (
+                (
+                    "normalised residue",
+                    f"(ineq_con{nums.icc[i]:03d})",
+                    -constraint.normalised_residual,
+                ),
+                (
+                    "physical value",
+                    f"(ineq_value_con{nums.icc[i]:03d})",
+                    constraint.constraint_value,
+                ),
+                (
+                    "symbol",
+                    f"(ineq_symbol_con{nums.icc[i]:03d})",
+                    f"'{constraint.symbol}'",
+                ),
+                ("units", f"(ineq_units_con{nums.icc[i]:03d})", f"'{constraint.units}'"),
+                (
+                    "physical bound",
+                    f"(ineq_bound_con{nums.icc[i]:03d})",
+                    constraint.constraint_bound,
+                ),
+            ):
+                process_output.ovarre(constants.MFILE, f"{name} {d}", var, v)
+
+        process_output.write(
+            constants.NOUT,
+            tabulate(
+                inequality_constraint_table,
+                headers=[
+                    "",
+                    "Physical constraint",
+                    "",
+                    "Physical constraint bound",
+                    "Constraint residue",
+                    "Normalised residue",
+                ],
+                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 1c564e07dd..38e5770147 100644
--- a/process/core/solver/solver.py
+++ b/process/core/solver/solver.py
@@ -16,6 +16,7 @@
 )
 from scipy.optimize import fsolve
 
+from process.core import constants, process_output
 from process.core.exceptions import ProcessValueError
 from process.core.model import DataStructure
 from process.core.solver.evaluators import Evaluators
@@ -293,6 +294,65 @@ def _ineq_cons_satisfied(
 
         return self.info
 
+    def verror(self):
+        """Routine to print out relevant messages in the case of an
+        unfeasible result from a VMCON (optimisation) run
+
+        This routine prints out relevant messages in the case of
+        an unfeasible result from a VMCON (optimisation) run.
+
+        Parameters
+        ----------
+        ifail: int :
+
+        """
+        strings = "\n".join(
+            {
+                -1: ("User-terminated execution of VMCON.",),
+                0: (
+                    "Improper input parameters to the VMCON routine.",
+                    "PROCESS coding must be checked.",
+                ),
+                2: (
+                    "The maximum number of calls has been reached without solution.",
+                    (
+                        "The code may be stuck in a minimum in the residual space that"
+                        " is significantly above zero.\n"
+                    ),
+                    "There is either no solution possible, or the code",
+                    "is failing to escape from a deep local minimum.",
+                    "Try changing the variables in IXC, or modify their initial values.",
+                ),
+                3: (
+                    "The line search required the maximum of 10 calls.",
+                    "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.",
+                ),
+                5: (
+                    "The quadratic programming technique was unable to",
+                    "find a feasible point.\n",
+                    "Try changing or adding variables to IXC, or modify",
+                    "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")
+        )
+
+        process_output.ocmmnt(constants.NOUT, strings)
+        print(strings)
+
 
 class VmconBounded(Vmcon):
     """A solver that uses VMCON but checks x is in bounds before running"""
diff --git a/process/core/solver/solver_handler.py b/process/core/solver/solver_handler.py
index 9a117bd864..251928c610 100644
--- a/process/core/solver/solver_handler.py
+++ b/process/core/solver/solver_handler.py
@@ -1,13 +1,18 @@
 """Module containing solver handler routines"""
 
+import logging
 from contextlib import contextmanager
 
+from process.core import constants, process_output
 from process.core.solver.evaluators import Evaluators
 from process.core.solver.iteration_variables import (
     load_iteration_variables,
     load_scaled_bounds,
 )
 from process.core.solver.solver import get_solver
+from process.data_structure.numerics import FiguresOfMerit, PROCESSRunMode
+
+logger = logging.getLogger(__name__)
 
 
 class SolverHandler:
@@ -92,6 +97,124 @@ def output(self):
         self.data.numerics.xcm[: self.solver.x.shape[0]] = self.solver.x
         self.data.numerics.rcm[: self.solver.conf.shape[0]] = self.solver.conf
 
+        nums = self.data.numerics
+
+        process_output.oheadr(constants.NOUT, "Numerics")
+        process_output.ocmmnt(
+            constants.NOUT,
+            f"PROCESS has performed a {'fsolve' if self.solver == 'fsolve' else 'VMCON'} (optimisation) run.",
+        )
+        ifail = self.solver.info
+        if ifail != 1:
+            process_output.ovarin(constants.NOUT, "Error flag", "(ifail)", ifail)
+            process_output.oheadr(
+                constants.IOTTY, "PROCESS COULD NOT FIND A FEASIBLE SOLUTION"
+            )
+            print()
+
+            logger.critical("Solver returns with ifail /= 1. %s", ifail)
+
+            if self.solver_name == "vmcon":
+                self.solver.verror()
+
+            process_output.oblnkl(constants.NOUT)
+            print()
+        else:
+            # Solution found
+            descr = "consistent" if self.solver == "fsolve" else "feasible"
+            process_output.ocmmnt(
+                constants.NOUT, f"and found a {descr} set of parameters."
+            )
+            process_output.oheadr(constants.IOTTY, f"PROCESS found a {descr} solution")
+            process_output.oblnkl(constants.NOUT)
+            process_output.ovarin(constants.NOUT, "Error flag", "(ifail)", ifail)
+
+            if nums.sqsumsq >= 1.0e-2:
+                string = (
+                    "WARNING: Constraint residues are HIGH; consider re-running\n"
+                    "   with lower values of EPSVMC to confirm convergence...\n"
+                    "   (should be able to get down to about 1.0E-8 okay)\n"
+                )
+                process_output.ocmmnt(constants.NOUT, ("\n" + string))
+                print(string)
+
+                logger.warning(f"High final constraint residues. {nums.sqsumsq=}")
+
+        for d, var, v in (
+            ("Number of iteration variables", "(nvar)", nums.nvar),
+            (
+                "Number of constraints (total)",
+                "(neqns+nineqns)",
+                nums.neqns + nums.nineqns,
+            ),
+            ("Optimisation switch", "(ioptimz)", nums.ioptimz),
+        ):
+            process_output.ovarin(constants.NOUT, d, var, v)
+
+        process_output.ocmmnt(
+            constants.NOUT,
+            f"     {PROCESSRunMode(nums.ioptimz).description}",
+        )
+
+        # Objective function output: none for fsolve
+        if self.solver_name != "fsolve":
+            process_output.ovarin(
+                constants.NOUT,
+                "Figure of merit switch",
+                "(minmax)",
+                nums.minmax,
+            )
+
+            nums.objf_name = f'"{FiguresOfMerit(abs(nums.minmax)).description}"'
+
+            for d, var, v, o in (
+                ("Objective function name", "(objf_name)", nums.objf_name, ""),
+                ("Normalised objective function", "(norm_objf)", nums.norm_objf, "OP "),
+                (
+                    "VMCON convergence parameter",
+                    "(convergence_parameter)",
+                    self.data.globals.convergence_parameter,
+                    "OP ",
+                ),
+                (
+                    "Number of optimising solver iterations",
+                    "(nviter)",
+                    nums.nviter,
+                    "OP ",
+                ),
+                (
+                    "Square root of the sum of squares of the constraint residuals",
+                    "(sqsumsq)",
+                    nums.sqsumsq,
+                    "OP ",
+                ),
+            ):
+                process_output.ovarre(constants.NOUT, d, var, v, o)
+
+        process_output.oblnkl(constants.NOUT)
+
+        if self.solver_name == "fsolve":
+            process_output.write(
+                constants.NOUT,
+                "PROCESS has solved using fsolve.\n"
+                if ifail == 1
+                else "PROCESS failed to solve using fsolve.\n",
+            )
+        else:
+            process_output.write(
+                constants.NOUT,
+                (
+                    (
+                        "PROCESS has successfully optimised"
+                        if ifail == 1
+                        else "PROCESS has failed to optimise"
+                    )
+                    + " the optimisation parameters to"
+                    + ("minimise" if nums.minmax > 0 else "maximise")
+                    + f" the objective function: {nums.objf_name}\n"
+                ),
+            )
+
 
 @contextmanager
 def epsfcn_context(numerics):
diff --git a/process/data_structure/global_variables.py b/process/data_structure/global_variables.py
index f3e9849c3e..100198ee33 100644
--- a/process/data_structure/global_variables.py
+++ b/process/data_structure/global_variables.py
@@ -1,4 +1,4 @@
-from dataclasses import dataclass
+from dataclasses import dataclass, field
 
 
 @dataclass(slots=True)
@@ -18,19 +18,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 5e0180b1b6..d381e1bbe8 100644
--- a/process/data_structure/scan_variables.py
+++ b/process/data_structure/scan_variables.py
@@ -9,51 +9,53 @@
 
 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:
-    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 009317ab68..9a66cf2175 100644
--- a/process/main.py
+++ b/process/main.py
@@ -55,6 +55,7 @@
 from process.core.model import DataStructure, Model
 from process.core.process_output import OutputFileManager, oheadr
 from process.core.scan import Scan
+from process.data_structure.blanket_variables import BlktModelTypes
 from process.data_structure.cost_variables import CostModels
 from process.data_structure.numerics import PROCESSRunMode
 from process.models.availability import Availability
@@ -110,7 +111,7 @@
 from process.models.stellarator.neoclassics import Neoclassics
 from process.models.stellarator.stellarator import Stellarator
 from process.models.structure import Structure
-from process.models.tfcoil.base import TFCoil
+from process.models.tfcoil.base import TFCoil, TFConductorModel
 from process.models.tfcoil.resistive import (
     AluminiumTFCoil,
     CopperTFCoil,
@@ -120,6 +121,7 @@
     CICCSuperconductingTFCoil,
     CROCOSuperconductingTFCoil,
     SuperconductingTFCoil,
+    SuperconductingTFTurnType,
 )
 from process.models.vacuum import Vacuum, VacuumVessel
 from process.models.water_use import WaterUse
@@ -453,7 +455,6 @@ def run_scan(self):
         # ioptimz == 1: optimisation
         if self.data.numerics.ioptimz == PROCESSRunMode.OPTIMISATION:
             pass
-        # ioptimz == -2: evaluation
         elif self.data.numerics.ioptimz == PROCESSRunMode.EVALUATION:
             # No optimisation:
             # solve equality (consistency) constraints only using fsolve (HYBRD)
@@ -464,6 +465,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():
@@ -837,6 +839,145 @@ def setup_data_structure(self):
         for model in self.models:
             model.data = self.data
 
+    def write(self, data, _outfile):
+        """Write the results to the main output file (OUT.DAT).
+
+        Write the program results to a file, in a tidy format.
+
+        Parameters
+        ----------
+        self : process.main.Models
+            physics and engineering model objects
+        _outfile : int
+            Fortran output unit identifier
+
+        """
+        # ensure we are capturing warnings that occur in the 'output' stage as these are warnings
+        # that occur at our solution point. So we clear existing warnings
+        logging_model_handler.start_capturing()
+        logging_model_handler.clear_logs()
+
+        # Call stellarator output routine instead if relevant
+        if data.stellarator.istell != 0:
+            self.stellarator.output()
+            return
+
+        #  Call IFE output routine instead if relevant
+        if data.ife.ife != 0:
+            self.ife.output()
+            return
+
+        # Costs model
+        # Cost switch values
+        # No.  |  model
+        # ---- | ------
+        # 0    |  1990 costs model
+        # 1    |  2015 Kovari model
+        # 2    |  Custom model
+        self.costs.output()
+
+        # Availability model
+        self.availability.output()
+
+        # Physics model
+        self.physics.output()
+
+        # Detailed physics, currently only done at final point as values are not used
+        # by any other functions
+        self.physics_detailed.output()
+
+        # TODO what is this? Not in caller.py?
+        self.current_drive.output()
+
+        # Pulsed reactor model
+        self.pulse.output()
+
+        self.divertor.output()
+
+        # Machine Build Model
+        self.build.output()
+
+        # Cryostat build
+        self.cryostat.output()
+
+        # Toroidal field coil copper model
+        if data.tfcoil.i_tf_sup == TFConductorModel.WATER_COOLED_COPPER:
+            self.copper_tf_coil.output()
+
+        # Toroidal field coil superconductor model
+        if data.tfcoil.i_tf_sup == TFConductorModel.SUPERCONDUCTING:
+            tf_turn_type = SuperconductingTFTurnType(
+                data.superconducting_tfcoil.i_tf_turn_type
+            )
+            if tf_turn_type == SuperconductingTFTurnType.CABLE_IN_CONDUIT:
+                self.cicc_sctfcoil.output()
+            elif tf_turn_type == SuperconductingTFTurnType.CROSS_CONDUCTOR:
+                self.croco_sctfcoil.output()
+            else:
+                raise ValueError(
+                    "Unsupported superconducting TF turn type: "
+                    f"{data.superconducting_tfcoil.i_tf_turn_type}"
+                )
+
+        # Toroidal field coil aluminium model
+        if data.tfcoil.i_tf_sup == TFConductorModel.HELIUM_COOLED_ALUMINIUM:
+            self.aluminium_tf_coil.output()
+
+        # Tight aspect ratio machine model
+        if (
+            data.physics.itart == 1
+            and data.tfcoil.i_tf_sup != TFConductorModel.SUPERCONDUCTING
+        ):
+            self.tfcoil.output()
+
+        # Poloidal field coil model
+        self.pfcoil.output()
+
+        # Structure Model
+        self.structure.output()
+
+        # Blanket model
+        # Blanket switch values
+        # No.  |  model
+        # ---- | ------
+        # 1    |  CCFE HCPB model
+        # 2    |  KIT HCPB model
+        # 3    |  CCFE HCPB model with Tritium Breeding Ratio calculation
+        # 4    |  KIT HCLL model
+        # 5    |  DCLL model
+
+        self.shield.output()
+        self.vacuum_vessel.output()
+
+        # First wall geometry
+        self.fw.output()
+
+        if data.fwbs.i_blanket_type == BlktModelTypes.CCFE_HCPB:
+            # CCFE HCPB model
+            self.ccfe_hcpb.output()
+
+        elif data.fwbs.i_blanket_type == BlktModelTypes.DCLL:
+            # DCLL model
+            self.dcll.output()
+
+        # FISPACT and LOCA model (not used)- removed
+
+        # Power model
+        self.power.output()
+
+        # Vacuum model
+        self.vacuum.output()
+
+        # Buildings model
+        self.buildings.output()
+
+        # Water usage in secondary cooling system
+        self.water_use.output()
+
+        # stop capturing warnings so that Outfile does not end up with
+        # a lot of non-model logs
+        logging_model_handler.stop_capturing()
+
 
 # setup handlers for writing to terminal (on warnings+)
 # or writing to the log file (on info+)