diff --git a/documentation/source/development/add-vars.md b/documentation/source/development/add-vars.md
index e925ea093d..703a648e52 100644
--- a/documentation/source/development/add-vars.md
+++ b/documentation/source/development/add-vars.md
@@ -1,6 +1,6 @@
 # Guide for adding variables and constraints
 
-A guide for how to add new inputs, constraints, figures of merit, and scan variables. 
+A guide for how to add new inputs, constraints, figures of merit, and scan variables.
 
 **At all times ensure new variables and functions adhere to the [`PROCESS` style guide](../development/standards.md).**
 
@@ -23,7 +23,7 @@ It is advised that a variable is either used to define a particular PROCESS run
 ## Add a new variable
 You may need to add a variable to PROCESS when changing or creating models. In most cases, you will want to add your variable to an existing data structure. Creating an entirely new data structure is beyond the scope of this guide, so please seek support from the PROCESS maintainers.
 
-For example, if you are adding a new variable that relates to the blanket model, you would add the variable to `process/data_structure/blanket_variables.py` as part of the `BlanketData` dataclass. 
+For example, if you are adding a new variable that relates to the blanket model, you would add the variable to `process/data_structure/blanket_variables.py` as part of the `BlanketData` dataclass.
 
 ```python
 @dataclass(slots=True)
@@ -116,7 +116,7 @@ class FiguresOfMerit(IntEnum):
     ...
     BLANKET_FIGURE_OF_MERIT = (20, "my FOM description")
 ```
-Here `20` will be the identifier of the figure of merit, and **must** be unique. 
+Here `20` will be the identifier of the figure of merit, and **must** be unique.
 
 Finally, add the equation to `process/core/solver/objectives.py`:
 ```python
@@ -136,26 +136,22 @@ Remember, setting `minmax = -20` would minimise instead of maximise our new vari
 
 ## Add a scan variable
 
-After following the instructions to add an input variable, you can then make a scan variable.
+After following the instruction to add an input variable, you can make the variable a scan variable by following these steps:
 
-First, add the variable to the `ScanVariables` enum in `process/core/scan.py`. 
-```python
-class ScanVariables(Enum):
-    ...
-    blanket_scan_variable = ScanVariable(
-        "my_new_blanket_variable", "A blanket variable", 82
-    )
-```
-Here, `82` is the identifier of the scan variable and must be unique. 
+1. Update the `ScanVariables` enum in the `scan.py` file by adding a new case statement connecting the variable to the scan integer switch, the variable name and a short description.
 
-Next, increment the parameter `IPNSCNV` in `process/data_structure/scan_variables.py` and be sure to add a description of the scan variable in the docstring of the `nsweep` variable.
+2. Add a comment in the corresponding variable file in the data_structure directory, eg, `data_structure/[XX]_variables.py`, to add the variable description indicating the scan switch number.
 
-Finally, in `process/core/scan.py`, add the scan variable to the `Scan.scan_select()` method.
+`SCAN_VARIABLES` case example:
+
+Add the variable to the `ScanVariables` enum in `process/core/scan.py`.
 ```python
-match nwp:
-  ...
-  case 82:
-    self.data.tfcoil.my_new_blanket_variable = swp[iscn - 1]
+  class ScanVariables(ScanVariable, Enum):
+    aspect = (1, SVE.P),
+    pflux_div_heat_load_max_mw = (2, SVE.D)
+    ...
+    b_crit_upper_nbti =  (54, SVE.T)
+    dr_shld_inboard = (55, SVE.B),
 ```
 
 Please see the [scan documentation](../io/input-guide.md#input-guide) for how to set up a scan `IN.DAT`
diff --git a/process/core/io/plot/cli.py b/process/core/io/plot/cli.py
index 804b969a7d..38d70f106d 100644
--- a/process/core/io/plot/cli.py
+++ b/process/core/io/plot/cli.py
@@ -1,3 +1,5 @@
+"""PROCESS cli"""
+
 from pathlib import Path
 
 import click
@@ -241,10 +243,10 @@ def plot_scans_cli(
         list(map(float, filter(None, x_axis_max))),
         list(map(float, filter(None, x_axis_range))),
         y_axis_percent,
-        y_axis_percent2,
         list(map(float, filter(None, y_axis_max))),
-        list(map(float, filter(None, y_axis2_max))),
         list(map(float, filter(None, y_axis_range))),
+        y_axis_percent2,
+        list(map(float, filter(None, y_axis2_max))),
         list(map(float, filter(None, y_axis_range2))),
         label_name,
         twod_contour,
diff --git a/process/core/io/plot/scans.py b/process/core/io/plot/scans.py
index 80ac996f60..e3da43d6fc 100644
--- a/process/core/io/plot/scans.py
+++ b/process/core/io/plot/scans.py
@@ -22,17 +22,150 @@
 - If the file is a folder, the contained MFILE is used as an input.
 """
 
+from __future__ import annotations
+
 import math
-import sys
 from collections.abc import Iterable, Sequence
+from dataclasses import dataclass
+from enum import Enum, auto
 from pathlib import Path
+from typing import TYPE_CHECKING
 
 import matplotlib.pyplot as plt
-import matplotlib.ticker as mtick
 import numpy as np
+from matplotlib.ticker import MultipleLocator, PercentFormatter
 
 from process.core.io.mfile import MFile
 from process.core.io.variable_metadata import var_dicts as meta
+from process.core.scan import ScanVariables
+
+if TYPE_CHECKING:
+    from matproplib.axes import Axes
+
+
+@dataclass
+class AxisData:
+    """Axis data container"""
+
+    name: str
+    percent: bool
+    max_: Sequence[float]
+    range_: Sequence[float]
+    tick_size: float
+    font_size: float
+    legend_size: float = 12
+
+
+class AxisChoice(Enum):
+    """Axis choices"""
+
+    X = auto()
+    Y = auto()
+
+    def axis(self, ax):
+        """Get correct axis"""
+        return getattr(ax, f"{self.name.lower()}axis")
+
+    def set_lim(self, ax, lower, upper):
+        """Set axis limits"""
+        getattr(ax, f"set_{self.name.lower()}lim")(lower, upper)
+
+
+def get_list_padded(inp, names):
+    """Pad list contents"""
+    target_len = len(names)
+    inp_array = np.array(inp, dtype=float)
+    if (i_len := len(inp_array)) < target_len:
+        if i_len == 0:
+            return [None] * target_len
+
+        return np.concatenate((
+            inp_array,
+            np.full(target_len - i_len, inp_array[-1], dtype=float),
+        ))
+    return inp_array[:target_len]
+
+
+def value_checks(
+    scan_var: ScanVariables,
+    scan_2_var: ScanVariables | None,
+    m_file: MFile,
+    input_files: Sequence[Path],
+):
+    """Check scan variable values
+
+    Raises
+    ------
+    ValueError
+        Scan variable not in MFILE
+    ValueError
+        Multiple input files specified for plotting
+    """
+    ve_string = (
+        "`{}` does not exist in PROCESS dicts\n"
+        " The scan variable ({}) is probably an upper/lower boundary\n"
+        " Please modify 'nsweep_dict' dict with the constrained var"
+    )
+    # Check if the scan variable is present
+    if scan_var.out_name not in m_file.data:
+        raise ValueError(ve_string.format(scan_var.out_name, scan_var.name))
+
+    # Check if the second scan variable is present
+    if scan_2_var is not None:
+        if scan_2_var.out_name not in m_file.data:
+            raise ValueError(ve_string.format(scan_2_var.out_name, scan_2_var.name))
+
+        if len(input_files) > 1:
+            raise ValueError("Only one input file can be used for 2D scans")
+
+
+def array_check(output_name: str, m_file: MFile) -> bool:
+    """Check if the output variable exists in the MFILE"""
+    if output_name not in m_file.data:
+        print(
+            f"Warning : `{output_name}` does not exist in PROCESS dicts\n"
+            f"Warning : `{output_name}` will not be output"
+        )
+        return False
+    return True
+
+
+def create_o_array(
+    n_scan: int, m_file: MFile, output_name: str, conv_i: list[int]
+) -> np.ndarray:
+    """Creating output array"""
+    return np.array([m_file.get(output_name, scan=conv_i[ii]) for ii in range(n_scan)])
+
+
+def get_label(name: str) -> str:
+    """Get latex label"""
+    return meta[name].latex if name in meta else f"{name}"
+
+
+def axis_manipulation(ax: Axes, axis: AxisData, index: int, contour: np.ndarray):
+    """Manipulate axes for unified layout"""
+    an = AxisChoice[axis.name.upper()]
+
+    if len(axis.range_) > 0:
+        divisions = (axis.range_[1] - axis.range_[0]) / 10
+    if axis.percent:
+        if axis.max_[index] is None:
+            axis.max_[index] = max(np.abs(contour))
+        ticks = PercentFormatter(axis.max_[index])
+        if len(axis.range_) > 0:
+            scale = axis.max_[index] / 100
+            divisions = 5 * math.ceil(divisions / 5) * scale
+            range_ = (axis.range_[0] * scale, axis.range_[1] * scale)
+        an.axis(ax).set_major_formatter(ticks)
+
+    if len(axis.range_) > 0:
+        if axis.percent is False:
+            range_ = axis.range_
+        an.set_lim(ax, range_[0], range_[1])
+        an.axis(ax).set_major_locator(MultipleLocator(divisions))
+
+    ax.figure.tight_layout()
+    ax.tick_params(axis=an.name.lower(), labelsize=axis.tick_size)
 
 
 def plot_scan(
@@ -48,10 +181,10 @@ def plot_scan(
     x_axis_max: Sequence[float] = (),
     x_axis_range: Sequence[float] = (),
     y_axis_percent: bool = False,
-    y_axis_percent2: bool = False,
     y_axis_max: Sequence[float] = (),
-    y_axis2_max: Sequence[float] = (),
     y_axis_range: Sequence[float] = (),
+    y_axis_percent2: bool = False,
+    y_axis2_max: Sequence[float] = (),
     y_axis_range2: Sequence[float] = (),
     label_name: Sequence[str] = (),
     twod_contour: bool = False,
@@ -60,908 +193,465 @@ def plot_scan(
     """Main plot scans script."""
     outputdir = outputdir or Path.cwd()
     input_files = mfiles if isinstance(mfiles, Iterable) else [mfiles]
-    x_max_input = x_axis_max
-
-    y_max_input = y_axis_max
-    y_max2_input = y_axis2_max
 
     # If the input file is a directory, add MFILE.DAT
     for ii, if_ in enumerate(input_files):
         if if_.is_dir():
             input_files[ii] = if_ / "MFILE.DAT"
 
-    # nsweep variable dict
-    # -------------------
-    # TODO WOULD BE GREAT TO HAVE IT AUTOMATICALLY GENERATED ON THE PROCESS CMAKE!
-    #        THE SAME WAY THE DICTS ARE
-    # This needs to be kept in sync automatically; this will break frequently
-    # otherwise
-    # Rem : Some variables are not in the MFILE, making the definition rather tricky...
-    nsweep_dict = {
-        1: "aspect",
-        2: "pflux_div_heat_load_max_mw",
-        3: "p_plant_electric_net_mw",
-        4: "hfact",
-        5: "j_tf_coil_full_area",
-        6: "pflux_fw_neutron_max_mw",
-        7: "beamfus0",
-        8: "Obsolete",  # OBSOLETE
-        9: "temp_plasma_electron_vol_avg_kev",
-        10: "boundu(15)",
-        11: "beta_norm_max",
-        12: "f_c_plasma_bootstrap_max",
-        13: "boundu(10)",
-        14: "f_j_tf_wp_critical_max",
-        16: "rmajor",
-        # b_tf_inboard_max the maximum T field upper limit is the scan variable
-        17: "b_tf_inboard_peak_symmetric",
-        18: "eta_cd_norm_hcd_primary_max",
-        19: "boundl(16)",
-        20: "cnstv.t_burn_min",
-        21: "",
-        22: "f_t_plant_available",
-        23: "boundu(72)",
-        24: "p_fusion_total_max_mw",
-        25: "kappa",
-        26: "triang",
-        27: "tbrmin",
-        28: "b_plasma_toroidal_on_axis",
-        29: "radius_plasma_core_norm",
-        30: "",  # OBSOLETE
-        31: "f_t_alpha_energy_confinement_min",
-        32: "epsvmc",
-        33: "ttarget",
-        34: "qtargettotal",
-        35: "lambda_q_omp",
-        36: "lambda_target",
-        37: "lcon_factor",
-        38: "boundu(129)",
-        39: "boundu(131)",
-        40: "boundu(135)",
-        41: "dr_blkt_outboard",
-        42: "f_nd_impurity_electrons(9)",
-        43: "Obsolete",  # OBSOLETE
-        44: "alstrtf",
-        45: "temp_tf_superconductor_margin_min",
-        46: "boundu(152)",
-        47: "impurity_enrichment(9)",
-        48: "n_tf_wp_pancakes",
-        49: "n_tf_wp_layers",
-        50: "f_nd_impurity_electrons(13)",
-        51: "f_p_div_lower",
-        52: "rad_fraction_sol",
-        53: "Obsolete",  # OBSOLETE
-        54: "b_crit_upper_nbti",
-        55: "dr_shld_inboard",
-        56: "p_cryo_plant_electric_max_mw",
-        # Genuinely b_plasma_toroidal_on_axis lower bound
-        57: "b_plasma_toroidal_on_axis",
-        58: "dr_fw_plasma_gap_inboard",
-        59: "dr_fw_plasma_gap_outboard",
-        60: "sig_tf_wp_max",
-        61: "copperaoh_m2_max",
-        62: "j_cs_flat_top_end",
-        63: "dr_cs",
-        64: "f_z_cs_tf_internal",
-        65: "n_cycle_min",
-        66: "f_a_cs_turn_steel",
-        67: "t_crack_vertical",
-        68: "inlet_temp_liq",
-        69: "outlet_temp_liq",
-        70: "blpressure_liq",
-        71: "n_liq_recirc",
-        72: "bz_channel_conduct_liq",
-        73: "pnuc_fw_ratio_dcll",
-        74: "f_nuc_pow_bz_struct",
-        75: "dx_fw_module",
-        76: "eta_turbine",
-        77: "startupratio",
-        78: "fkind",
-        79: "eta_ecrh_injector_wall_plug",
-        80: "fcoolcp",
-        81: "n_tf_coil_turns",
-    }
-    # -------------------
-
     # Getting the scanned variable name
     m_file = MFile(filename=input_files[-1])
-    nsweep_ref = int(m_file.data["nsweep"].get_scan(-1))
-    scan_var_name = nsweep_dict[nsweep_ref]
-    # Get the eventual second scan variable
-    nsweep_2_ref = 0
-    is_2D_scan = False
-    scan_2_var_name = ""
-    if "nsweep_2" in m_file.data:
-        is_2D_scan = True
-        nsweep_2_ref = int(m_file.data["nsweep_2"].get_scan(-1))
-        scan_2_var_name = nsweep_dict[nsweep_2_ref]
-
-    # Checks
-    # ------
-    # Check if the nsweep dict has been updated
-    if nsweep_ref > len(nsweep_dict) + 1:
-        print(
-            f"ERROR : nsweep = {nsweep_ref} not supported by the utility\n"
-            "ERROR : Please update the 'nsweep_dict' dict"
-        )
-        sys.exit()
+    nsweep_ref = int(m_file.get("nsweep", scan=-1))
+    scan_var = ScanVariables(nsweep_ref)
 
-    # Check if the scan variable is present in the
-    if scan_var_name not in m_file.data:
-        print(
-            f"ERROR : `{scan_var_name}` does not exist in PROCESS dicts\n"
-            "ERROR : The scan variable is probably an upper/lower boundary\n"
-            "ERROR : Please modify 'nsweep_dict' dict with the constrained var"
-        )
-        sys.exit()
-
-    # Check if the second scan variable is present in the
-    if is_2D_scan and (scan_2_var_name not in m_file.data):
-        print(
-            f"ERROR : `{scan_2_var_name}` does not exist in PROCESS dicts\n"
-            "ERROR : The scan variable is probably an upper/lower boundary\n"
-            "ERROR : Please modify 'nsweep_dict' dict with the constrained var"
-        )
-        sys.exit()
-
-    # Only one input must be used for a 2D scan
-    if is_2D_scan and len(input_files) > 1:
-        print("ERROR : Only one input file can be used for 2D scans\nERROR : Exiting")
-        sys.exit()
-    # ------
-
-    # Plot settings
-    # -------------
-    # Plot cosmetic settings
-    def _format_lists(inp, output_names):
-        x_max = []
-        if len(inp) > 0:
-            for i in range(len(output_names)):
-                j = 0
-                try:
-                    x_max += [float(inp[i])]
-                    j += 1
-                except IndexError:
-                    x_max += [float(inp[j])]
-        else:
-            x_max = [None] * len(output_names)
+    # Get the eventual second scan variable
+    scan_2_var = (
+        ScanVariables(int(m_file.get("nsweep_2", scan=-1)))
+        if "nsweep_2" in m_file.data
+        else None
+    )
 
-        return x_max
+    value_checks(scan_var, scan_2_var, m_file, input_files)
 
-    legend_size = 12
-    x_max = (
-        _format_lists(x_max_input, output_names)
-        if len(x_max_input) != len(output_names)
-        else np.float64(x_max_input)
+    x_max = get_list_padded(x_axis_max, output_names)
+    x_axis = AxisData(
+        "x", x_axis_percent, x_max, x_axis_range, axis_tick_size, axis_font_size
     )
-    y_max = (
-        _format_lists(y_max_input, output_names)
-        if len(y_max_input) != len(output_names)
-        else np.float64(y_max_input)
+
+    y_max = get_list_padded(y_axis_max, output_names)
+    y_axis = AxisData(
+        "y", y_axis_percent, y_max, y_axis_range, axis_tick_size, axis_font_size
     )
-    if len(output_names2) > 0:
-        y_max2 = (
-            _format_lists(y_max2_input, output_names)
-            if len(y_max2_input) != len(output_names)
-            else np.float64(y_max2_input)
+    if scan_2_var is None:
+        y_axis2 = AxisData(
+            "y",
+            y_axis_percent2,
+            (
+                get_list_padded(y_axis2_max, output_names)
+                if len(output_names2) > 0
+                else y_axis2_max
+            ),
+            y_axis_range2,
+            axis_tick_size,
+            axis_font_size,
+        )
+
+        scan_var_array, output_arrays, output_arrays2 = oned_scan(
+            input_files,
+            nsweep_ref,
+            scan_var,
+            output_names,
+            output_names2,
+            term_output=term_output,
+        )
+        plot_1d_scan(
+            input_files,
+            m_file,
+            output_names,
+            output_names2,
+            scan_var,
+            outputdir,
+            save_format,
+            label_name,
+            scan_var_array,
+            output_arrays,
+            output_arrays2,
+            x_axis,
+            y_axis,
+            y_axis2,
+            stack_plots=stack_plots,
         )
     else:
-        y_max2 = y_max2_input
-    # -------------
-
-    # Case of a set of 1D scans
-    if not is_2D_scan:
-        # Loop over the MFILEs
-        output_arrays = {}
-        output_arrays2 = {}
-        scan_var_array = {}
-        for input_file in input_files:
-            # Opening the MFILE.DAT
-            m_file = MFile(filename=input_file)
-
-            # Check if the the scan variable is the same for all inputs
-            # ---
-            # Same scan var
-            nsweep = int(m_file.data["nsweep"].get_scan(-1))
-            if nsweep != nsweep_ref:
-                print(
-                    "ERROR : You must use inputs files with the same scan variables\n"
-                    "ERROR : Exiting"
-                )
-                sys.exit()
-
-            # No D scans
-            if "nsweep_2" in m_file.data:
-                print("ERROR : You cannot mix 1D with 2D scans\nERROR : Exiting")
-                sys.exit()
-            # ---
-
-            # Only selecting the scans that has converged
-            # ---
-            # Number of scan points
-            n_scan = int(m_file.data["isweep"].get_scan(-1))
-
-            # Converged indexes
-            conv_i = []
-            for ii in range(n_scan):
-                ifail = m_file.data["ifail"].get_scan(ii + 1)
-                if ifail == 1:
-                    conv_i.append(ii + 1)
-                else:
-                    failed_value = m_file.data[scan_var_name].get_scan(ii + 1)
-                    print(
-                        "Warning : Non-convergent scan point : "
-                        f"{scan_var_name} = {failed_value}\n"
-                        "Warning : This point will not be shown."
-                    )
+        twod_scan(
+            input_files,
+            scan_var,
+            scan_2_var,
+            output_names,
+            outputdir,
+            save_format,
+            x_axis,
+            y_axis,
+            twod_contour=twod_contour,
+        )
 
-            # Updating the number of scans
-            n_scan = len(conv_i)
-            # ---
-            # Scanned variable
-            scan_var_array[input_file] = np.zeros(n_scan)
 
-            for ii in range(n_scan):
-                scan_var_array[input_file][ii] = m_file.data[scan_var_name].get_scan(
-                    conv_i[ii]
+def oned_scan(
+    input_files: Sequence[Path],
+    nsweep_ref: int,
+    scan_var: ScanVariables,
+    output_names: Sequence[str],
+    output_names2: Sequence[str],
+    *,
+    term_output: bool,
+) -> tuple[np.ndarray, ...]:
+    """One d scan plotting setup
+
+    Raises
+    ------
+    ValueError
+        Input files must have the same scan variables
+    ValueError
+        Mixing 1D and 2D scan files
+    """
+    # input file, output_name, scan
+    output_arrays = {}
+    # input file, output_name2, scan
+    output_arrays2 = {}
+    # input_file, scan
+    scan_var_array = {}
+    for input_file in input_files:
+        # Opening the MFILE.DAT
+        m_file = MFile(filename=input_file)
+        n_scan = int(m_file.get("isweep", scan=-1))
+
+        # Check if the the scan variable is the same for all inputs
+        # Same scan var
+        nsweep = int(m_file.get("nsweep", scan=-1))
+        if nsweep != nsweep_ref:
+            raise ValueError("You must use inputs files with the same scan variables\n")
+
+        if "nsweep_2" in m_file.data:
+            raise ValueError("You cannot mix 1D with 2D scans\nERROR : Exiting")
+
+        # Only selecting the scans that has converged
+        # Converged indexes
+        conv_i = []
+        for ii in range(n_scan):
+            ifail = m_file.get("ifail", scan=ii + 1)
+            if ifail == 1:
+                conv_i.append(ii + 1)
+            else:
+                failed_value = scan_var.get_val(m_file, scan=ii + 1)
+                print(
+                    "Warning : Non-convergent scan point : "
+                    f"{scan_var.name} = {failed_value}\n"
+                    "Warning : This point will not be shown."
                 )
-            # output list declaration
-            output_arrays[input_file] = {}
-            output_arrays2[input_file] = {}
-            # First variable scan
-            for output_name in output_names:
-                ouput_array = np.zeros(n_scan)
-
-                # Check if the output variable exists in the MFILE
-                if output_name not in m_file.data:
-                    print(
-                        f"Warning : `{output_name}` does not exist in PROCESS dicts\n"
-                        f"Warning : `{output_name}` will not be output"
-                    )
-                    continue
-
-                for ii in range(n_scan):
-                    ouput_array[ii] = m_file.data[output_name].get_scan(conv_i[ii])
-                output_arrays[input_file][output_name] = ouput_array
-            # Second variable scan
-            for output_name2 in output_names2:
-                ouput_array2 = np.zeros(n_scan)
-
-                # Check if the output variable exists in the MFILE
-                if output_name2 not in m_file.data:
-                    print(
-                        f"Warning : `{output_name2}` does not exist in PROCESS dicts\n"
-                        f"Warning : `{output_name2}` will not be output"
-                    )
-                    continue
 
-                for ii in range(n_scan):
-                    ouput_array2[ii] = m_file.data[output_name2].get_scan(conv_i[ii])
-                output_arrays2[input_file][output_name2] = ouput_array2
-            # Terminal output
-            if term_output:
+        # Updating the number of scans
+        n_scan = len(conv_i)
+        scan_var_array[input_file] = np.array([
+            scan_var.get_val(m_file, scan=conv_i[ii]) for ii in range(n_scan)
+        ])
+        output_arrays[input_file] = {
+            output_name: create_o_array(n_scan, m_file, output_name, conv_i)
+            for output_name in output_names
+            if array_check(output_name, m_file)
+        }
+        output_arrays2[input_file] = {
+            output_name2: create_o_array(n_scan, m_file, output_name2, conv_i)
+            for output_name2 in output_names2
+            if array_check(output_name2, m_file)
+        }
+        # Terminal output
+        if term_output:
+            print(
+                f"\n{input_file} scan output\n\nX-axis:\n"
+                f"scan var {scan_var.name} : {scan_var_array[input_file]}\n\nY-axis:"
+                + "\n".join(
+                    f"{output_name} : {output_arrays[input_file][output_name]}"
+                    for output_name in output_names
+                    if output_name in m_file.data
+                )
+                + "\n"
+            )
+            if len(output_names2) > 0:
+                last_name = output_names2[-1]
                 print(
-                    f"\n{input_file} scan output\n\nX-axis:\n"
-                    f"scan var {scan_var_name} : {scan_var_array[input_file]}\n\nY-axis:"
+                    f"Y2-Axis\n  {last_name} : {output_arrays2[input_file][last_name]}\n"
                 )
-                for output_name in output_names:
-                    # Check if the output variable exists in the MFILE
-                    if output_name not in m_file.data:
-                        continue
-
-                    print(f"{output_name} : {output_arrays[input_file][output_name]}")
-                print()
-                if len(output_names2) > 0:
-                    print(
-                        f"Y2-Axis\n  {output_name2} : "
-                        f"{output_arrays2[input_file][output_name2]}\n"
-                    )
-        # Plot section
-        # -----------
-        for index, output_name in enumerate(output_names):
-            if not stack_plots:
-                fig, ax = plt.subplots()
-                if len(output_names2) > 0:
-                    ax2 = ax.twinx()
-            # reset counter for label_name
-            kk = 0
-
-            # Check if the output variable exists in the MFILE
-            if output_name not in m_file.data:
-                continue
-
-            # Loop over inputs
-            for input_file in input_files:
-                # Legend label formating
-                if len(label_name) == 0:
-                    labl = input_file.name
-                else:
-                    labl = label_name[kk]
-                    kk += 1
-
-                # Plot the graph
-                if len(output_names2) > 0 and not stack_plots:
-                    ax.plot(
-                        scan_var_array[input_file],
-                        output_arrays[input_file][output_name],
-                        "--o",
-                        color="blue" if len(input_files) == 1 else None,
-                        label=labl,
-                    )
-                    if len(y_axis_range) > 0:
-                        y_divisions = (y_axis_range[1] - y_axis_range[0]) / 10
-                    if y_axis_percent:
-                        if y_max[index] is None:
-                            y_max[index] = max(
-                                np.abs(output_arrays[input_file][output_name])
-                            )
-                        yticks = mtick.PercentFormatter(y_max[index])
-                        if len(y_axis_range) > 0:
-                            y_divisions = (
-                                5 * math.ceil(y_divisions / 5) * y_max[index] / 100
-                            )
-                            y_range = (
-                                y_axis_range[0] * y_max[index] / 100,
-                                y_axis_range[1] * y_max[index] / 100,
-                            )
-                        ax.yaxis.set_major_formatter(yticks)
-                    if len(y_axis_range) > 0:
-                        if y_axis_percent is False:
-                            y_range = y_axis_range
-                        ax.set_ylim(y_range[0], y_range[1])
-                        ax.yaxis.set_major_locator(mtick.MultipleLocator(y_divisions))
-                    if len(x_axis_range) > 0:
-                        x_divisions = (x_axis_range[1] - x_axis_range[0]) / 10
-                    if x_axis_percent:
-                        if x_max[index] is None:
-                            x_max[index] = max(np.abs(scan_var_array[input_file]))
-                        xticks = mtick.PercentFormatter(x_max[index])
-                        ax.xaxis.set_major_formatter(xticks)
-                        if len(x_axis_range) > 0:
-                            x_divisions = (
-                                5 * math.ceil(x_divisions / 5) * x_max[index] / 100
-                            )
-                            x_range = (
-                                x_axis_range[0] * x_max[index] / 100,
-                                x_axis_range[1] * x_max[index] / 100,
-                            )
-                    plt.rc("xtick", labelsize=axis_tick_size)
-                    plt.rc("ytick", labelsize=axis_tick_size)
-                    if len(x_axis_range) > 0:
-                        if x_axis_percent is False:
-                            x_range = x_axis_range
-                        plt.xlim(x_range[0], x_range[1])
-                        ax.xaxis.set_major_locator(mtick.MultipleLocator(x_divisions))
-                    plt.tight_layout()
-                elif stack_plots:
-                    # check stack plots will work
-                    if len(output_names) <= 1:
-                        raise ValueError(
-                            "stack_plots requires at least two output variables"
-                        )
-                    # Create subplots only once for the first output
-                    if index == 0:
-                        fig, axs = plt.subplots(
-                            len(output_names),
-                            1,
-                            figsize=(8.0, (3.5 + (1 * len(output_names)))),
-                            sharex=True,
-                        )
-                        fig.subplots_adjust(hspace=0.0)
-
-                    axs[index].plot(
-                        scan_var_array[input_file],
-                        output_arrays[input_file][output_name],
-                        "--o",
-                        color="blue" if len(output_names2) > 0 else None,
-                        label=labl,
-                    )
-                    if len(y_axis_range) > 0:
-                        y_divisions = (y_axis_range[1] - y_axis_range[0]) / 10
-                    if y_axis_percent:
-                        if y_max[index] is None:
-                            y_max[index] = max(
-                                np.abs(output_arrays[input_file][output_name])
-                            )
-                        yticks = mtick.PercentFormatter(y_max[index])
-                        if len(y_axis_range) > 0:
-                            y_divisions = (
-                                5 * math.ceil(y_divisions / 5) * y_max[index] / 100
-                            )
-                            y_range = (
-                                y_axis_range[0] * y_max[index] / 100,
-                                y_axis_range[1] * y_max[index] / 100,
-                            )
-                        axs[output_names.index(output_name)].yaxis.set_major_formatter(
-                            yticks
-                        )
-                    if len(y_axis_range) > 0:
-                        if y_axis_percent is False:
-                            y_range = y_axis_range
-                        axs[output_names.index(output_name)].set_ylim(
-                            y_range[0], y_range[1]
-                        )
-                        axs[output_names.index(output_name)].yaxis.set_major_locator(
-                            mtick.MultipleLocator(y_divisions)
-                        )
-                    if len(x_axis_range) > 0:
-                        x_divisions = (x_axis_range[1] - x_axis_range[0]) / 10
-                    if x_axis_percent:
-                        if x_max[index] is None:
-                            x_max[index] = max(np.abs(scan_var_array[input_file]))
-                        xticks = mtick.PercentFormatter(x_max[index])
-                        if len(x_axis_range) > 0:
-                            x_divisions = (
-                                5 * math.ceil(x_divisions / 5) * x_max[index] / 100
-                            )
-                            x_range = (
-                                x_axis_range[0] * x_max[index] / 100,
-                                x_axis_range[1] * x_max[index] / 100,
-                            )
-                        axs[output_names.index(output_name)].xaxis.set_major_formatter(
-                            xticks
-                        )
-                    if len(x_axis_range) > 0:
-                        if x_axis_percent is False:
-                            x_range = x_axis_range
-                        plt.xlim(x_range[0], x_range[1])
-                        axs[output_names.index(output_name)].xaxis.set_major_locator(
-                            mtick.MultipleLocator(x_divisions)
-                        )
-                    plt.rc("xtick", labelsize=axis_tick_size)
-                    plt.rc("ytick", labelsize=axis_tick_size)
-                    plt.tight_layout()
-                else:
-                    ax.plot(
-                        scan_var_array[input_file],
-                        output_arrays[input_file][output_name],
-                        "--o",
-                        color="blue" if len(output_names2) > 0 else None,
-                        label=labl,
-                    )
-                    if len(y_axis_range) > 0:
-                        y_divisions = (y_axis_range[1] - y_axis_range[0]) / 10
-                    if y_axis_percent:
-                        if y_max[index] is None:
-                            y_max[index] = max(
-                                np.abs(output_arrays[input_file][output_name])
-                            )
-                        yticks = mtick.PercentFormatter(y_max[index])
-                        if len(y_axis_range) > 0:
-                            y_divisions = (
-                                5 * math.ceil(y_divisions / 5) * y_max[index] / 100
-                            )
-                            y_range = (
-                                y_axis_range[0] * y_max[index] / 100,
-                                y_axis_range[1] * y_max[index] / 100,
-                            )
-                        ax.yaxis.set_major_formatter(yticks)
-                    if len(y_axis_range) > 0:
-                        if y_axis_percent is False:
-                            y_range = y_axis_range
-                        ax.set_ylim(y_range[0], y_range[1])
-                        ax.yaxis.set_major_locator(mtick.MultipleLocator(y_divisions))
-                    if len(x_axis_range) > 0:
-                        x_divisions = (x_axis_range[1] - x_axis_range[0]) / 10
-                    if x_axis_percent:
-                        if x_max[index] is None:
-                            x_max[index] = max(np.abs(scan_var_array[input_file]))
-                        xticks = mtick.PercentFormatter(x_max[index])
-                        if len(x_axis_range) > 0:
-                            x_divisions = (
-                                5 * math.ceil(x_divisions / 5) * x_max[index] / 100
-                            )
-                            x_range = (
-                                x_axis_range[0] * x_max[index] / 100,
-                                x_axis_range[1] * x_max[index] / 100,
-                            )
-                        ax.xaxis.set_major_formatter(xticks)
-                    if len(x_axis_range) > 0:
-                        if x_axis_percent is False:
-                            x_range = x_axis_range
-                        plt.xlim(x_range[0], x_range[1])
-                        ax.xaxis.set_major_locator(mtick.MultipleLocator(x_divisions))
-                    plt.rc("xtick", labelsize=axis_tick_size)
-                    plt.rc("ytick", labelsize=axis_tick_size)
-                    plt.tight_layout()
-                if len(output_names2) > 0:
+    return scan_var_array, output_arrays, output_arrays2
+
+
+def plot_1d_scan(
+    input_files: Sequence[Path],
+    m_file: MFile,
+    output_names: Sequence[str],
+    output_names2: Sequence[str],
+    scan_var: ScanVariables,
+    outputdir: Path,
+    save_format: str,
+    label_name: Sequence[str],
+    scan_var_array: np.ndarray,
+    output_arrays: np.ndarray,
+    output_arrays2: np.ndarray,
+    x_axis: AxisData,
+    y_axis: AxisData,
+    y_axis2: AxisData,
+    *,
+    stack_plots: bool,
+):
+    """One d scan plotting
+
+    Raises
+    ------
+    ValueError
+        Stack plots requires >=2 variables
+    """
+    if stack_plots:
+        # check stack plots will work
+        if len(output_names) <= 1:
+            raise ValueError("stack_plots requires at least two output variables")
+        # Create subplots only once for the first output
+        fig, axs = plt.subplots(
+            len(output_names),
+            1,
+            figsize=(8.0, (3.5 + (1 * len(output_names)))),
+            sharex=True,
+        )
+        fig.subplots_adjust(hspace=0.0)
+
+    colour = (  # be careful changing this
+        ("blue" if len(output_names2) > 0 else None)
+        if len(output_names2) <= 0 or stack_plots
+        else ("blue" if len(input_files) == 1 else None)
+    )
+
+    for index, output_name in enumerate(output_names):
+        # reset counter for label_name
+        kk = 0
+
+        if output_name not in m_file.data:
+            continue
+
+        if stack_plots:
+            ax_ = axs[index]
+            ax = axs[output_names.index(output_name)]
+        else:
+            fig, ax = plt.subplots()
+            if len(output_names2) > 0:
+                ax2 = ax.twinx()
+            ax_ = ax
+
+        for input_file in input_files:
+            if len(label_name) == 0:
+                labl = input_file.name
+            else:
+                labl = label_name[kk]
+                kk += 1
+
+            ax_.plot(
+                scan_var_array[input_file],
+                output_arrays[input_file][output_name],
+                "--o",
+                color=colour,
+                label=labl,
+            )
+
+            axis_manipulation(
+                ax, axis=x_axis, index=index, contour=scan_var_array[input_file]
+            )
+            axis_manipulation(
+                ax,
+                axis=y_axis,
+                index=index,
+                contour=output_arrays[input_file][output_name],
+            )
+
+            if len(output_names2) > 0:
+                for output_name2 in output_names2:
+                    yval = output_arrays2[input_file][output_name2]
+                    colour = "red" if len(input_files) == 1 else None
                     ax2.plot(
                         scan_var_array[input_file],
-                        output_arrays2[input_file][output_name2],
+                        yval,
                         "--o",
-                        color="red" if len(input_files) == 1 else None,
+                        color=colour,
                         label=labl,
                     )
                     ax2.set_ylabel(
-                        (
-                            meta[output_name2].latex
-                            if output_name2 in meta
-                            else f"{output_name2}"
-                        ),
-                        fontsize=axis_font_size,
-                        color="red" if len(input_files) == 1 else "black",
+                        get_label(output_name2),
+                        fontsize=y_axis.font_size,
+                        color=colour or "black",
                     )
-                    if len(y_axis_range2) > 0:
-                        y_divisions2 = (y_axis_range2[1] - y_axis_range2[0]) / 10
-                    if y_axis_percent2:
-                        if y_max2[index] is None:
-                            y_max2[index] = max(
-                                np.abs(output_arrays2[input_file][output_name2])
-                            )
-                        yticks2 = mtick.PercentFormatter(y_max2[index])
-                        if len(y_axis_range2) > 0:
-                            y_divisions2 = (
-                                5 * math.ceil(y_divisions2 / 5) * y_max2[index] / 100
-                            )
-                            y_range2 = (
-                                y_axis_range2[0] * y_max2[index] / 100,
-                                y_axis_range2[1] * y_max2[index] / 100,
-                            )
-                        ax2.yaxis.set_major_formatter(yticks2)
-                    if len(y_axis_range2) > 0:
-                        if y_axis_percent2 is False:
-                            y_range2 = y_axis_range2
-                        ax2.set_ylim(y_range2[0], y_range2[1])
-                        ax2.yaxis.set_major_locator(mtick.MultipleLocator(y_divisions2))
-                    plt.rc("xtick", labelsize=axis_tick_size)
-                    plt.rc("ytick", labelsize=axis_tick_size)
-                    plt.tight_layout()
-            if len(output_names2) > 0:
-                ax2.yaxis.grid(True)
-                ax.xaxis.grid(True)
-                ax.set_ylabel(
-                    (
-                        meta[output_name].latex
-                        if output_name in meta
-                        else f"{output_name}"
-                    ),
-                    fontsize=axis_font_size,
-                    color="blue" if len(input_files) == 1 else "black",
-                )
-                ax.set_xlabel(
-                    (
-                        meta[scan_var_name].latex
-                        if scan_var_name in meta
-                        else f"{scan_var_name}"
-                    ),
-                    fontsize=axis_font_size,
-                )
-                plt.rc("xtick", labelsize=axis_tick_size)
-                plt.rc("ytick", labelsize=axis_tick_size)
-                if len(input_files) != 1:
-                    plt.legend(loc="best", fontsize=legend_size)
-                plt.tight_layout()
-            elif stack_plots:
-                axs[output_names.index(output_name)].minorticks_on()
-                axs[output_names.index(output_name)].grid(True)
-                axs[output_names.index(output_name)].set_ylabel(
-                    (
-                        meta[output_name].latex
-                        if output_name in meta
-                        else f"{output_name}"
-                    ),
+                    axis_manipulation(ax2, axis=y_axis2, index=index, contour=yval)
+        if len(output_names2) > 0:
+            ax2.yaxis.grid(True)
+            ax.xaxis.grid(True)
+            ax.set_ylabel(
+                get_label(output_name),
+                fontsize=y_axis.font_size,
+                color="blue" if len(input_files) == 1 else "black",
+            )
+            ax.set_xlabel(
+                get_label(scan_var.name),
+                fontsize=x_axis.font_size,
+            )
+            if len(input_files) != 1:
+                fig.legend(loc="best", fontsize=x_axis.legend_size)
+        elif stack_plots:
+            ax.minorticks_on()
+            ax.grid(True)
+            ax.set_ylabel(get_label(output_name))
+            ax.set_xlabel(get_label(scan_var.name), fontsize=x_axis.font_size)
+
+            ymin, ymax = ax.get_ylim()
+            if ymin < 0 and ymax > 0:
+                mod_min = ymin * 1.1
+                mod_max = ymax * 1.1
+            elif ymin >= 0:
+                mod_min = ymin * 0.9
+                mod_max = ymax * 1.1
+            else:
+                mod_min = ymin * 1.1
+                mod_max = ymax * 0.9
+            ax.set_ylim(mod_min, mod_max)
+
+            if len(input_files) > 1:
+                fig.legend(
+                    loc="lower center",
+                    fontsize=x_axis.legend_size,
+                    bbox_to_anchor=(0.5, -0.5 - (0.1 * len(output_names))),
+                    fancybox=True,
+                    shadow=False,
+                    ncol=len(input_files),
+                    columnspacing=0.8,
                 )
 
-                plt.xlabel(
-                    (
-                        meta[scan_var_name].latex
-                        if scan_var_name in meta
-                        else f"{scan_var_name}"
-                    ),
-                    fontsize=axis_font_size,
-                )
-                plt.rc("xtick", labelsize=axis_tick_size)
-                plt.rc("ytick", labelsize=axis_tick_size)
-                if len(input_files) > 1:
-                    plt.legend(
-                        loc="lower center",
-                        fontsize=legend_size,
-                        bbox_to_anchor=(0.5, -0.5 - (0.1 * len(output_names))),
-                        fancybox=True,
-                        shadow=False,
-                        ncol=len(input_files),
-                        columnspacing=0.8,
-                    )
-                plt.tight_layout()
-                ymin, ymax = axs[output_names.index(output_name)].get_ylim()
-                if ymin < 0 and ymax > 0:
-                    axs[output_names.index(output_name)].set_ylim(ymin * 1.1, ymax * 1.1)
-                elif ymin >= 0:
-                    axs[output_names.index(output_name)].set_ylim(ymin * 0.9, ymax * 1.1)
-                else:
-                    axs[output_names.index(output_name)].set_ylim(ymin * 1.1, ymax * 0.9)
+        else:
+            ax.grid(True)
+            ax.set_ylabel(
+                get_label(output_name),
+                fontsize=x_axis.font_size,
+                color="red" if len(output_names2) > 0 else "black",
+            )
+            ax.set_xlabel(get_label(scan_var.name), fontsize=x_axis.font_size)
+
+            fig.suptitle(
+                f"{get_label(output_name)} vs {get_label(scan_var.name)}",
+                fontsize=x_axis.font_size,
+            )
+            if len(input_files) != 1:
+                fig.legend(loc="best", fontsize=x_axis.legend_size)
+
+        ax.tick_params(axis=x_axis.name.lower(), labelsize=x_axis.tick_size)
+        ax.tick_params(axis=y_axis.name.lower(), labelsize=y_axis.tick_size)
+        fig.tight_layout()
+
+        # Output file naming
+        # 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]
+        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]
+        else:
+            extra_str = (
+                f"{output_name}{f'_vs_{output_name2}' if len(output_names2) > 0 else ''}"
+            )
+
+        if (not stack_plots) or (stack_plots and output_names[-1] == output_name):
+            fig.savefig(
+                f"{outputdir}/scan_{scan_var.name}_vs_{extra_str}.{save_format}",
+                dpi=300,
+            )
+            plt.show()
+
+
+def twod_scan(
+    input_files: Sequence[Path],
+    scan_var: ScanVariables,
+    scan_2_var: ScanVariables,
+    output_names: Sequence[str],
+    outputdir: Path,
+    save_format: str,
+    x_axis: AxisData,
+    y_axis: AxisData,
+    *,
+    twod_contour: bool,
+):
+    """2D scan plottings"""
+    m_file = MFile(filename=input_files[0])
+
+    # Number of scan points
+    n_scan_1 = int(m_file.get("isweep", scan=-1))
+    n_scan_2 = int(m_file.get("isweep_2", scan=-1))
+    # Selecting the converged runs only
+    contour_conv_ij = []  # List of non-converged scan point numbers
+    # 2D array of converged scan point numbers (sweep = rows, sweep_2 = columns)
+    conv_ij = []
+    ii_jj = 0
+    for ii in range(n_scan_1):
+        conv_ij.append([])
+        for _jj in range(n_scan_2):
+            ii_jj += 1  # Represents the scan point number in the MFILE
+            ifail = m_file.get("ifail", scan=ii_jj)
+            if ifail == 1:
+                conv_ij[ii].append(ii_jj)  # Only appends scan number if scan converged
+                contour_conv_ij.append(ii_jj)
             else:
-                plt.grid(True)
-                plt.ylabel(
-                    (
-                        meta[output_name].latex
-                        if output_name in meta
-                        else f"{output_name}"
-                    ),
-                    fontsize=axis_font_size,
-                    color="red" if len(output_names2) > 0 else "black",
-                )
-                plt.xlabel(
-                    (
-                        meta[scan_var_name].latex
-                        if scan_var_name in meta
-                        else f"{scan_var_name}"
-                    ),
-                    fontsize=axis_font_size,
-                )
-                plt.rc("xtick", labelsize=axis_tick_size)
-                plt.rc("ytick", labelsize=axis_tick_size)
-                plt.title(
-                    f"{meta[output_name].latex if output_name in meta else {output_name}} vs "  # noqa: E501
-                    f"{meta[scan_var_name].latex if scan_var_name in meta else {scan_var_name}}",  # noqa: E501
-                    fontsize=axis_font_size,
+                failed_value_1 = scan_var.get_val(m_file, scan=ii_jj)
+                failed_value_2 = scan_2_var.get_val(m_file, scan=ii_jj)
+                print(
+                    "Warning : Non-convergent scan point : "
+                    f"({scan_var.name},{scan_2_var.name}) "
+                    f"= ({failed_value_1},{failed_value_2})\n"
+                    "Warning : This point will not be shown."
                 )
-                plt.tight_layout()
-                if len(input_files) != 1:
-                    plt.legend(loc="best", fontsize=legend_size)
-                    plt.tight_layout()
-
-            # Output file naming
-            if output_name == "plasma_current_MA":
-                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)}"  # noqa: E501
-            else:
-                extra_str = f"{output_name}{f'_vs_{output_name2}' if len(output_names2) > 0 else ''}"  # noqa: E501
 
-            if (not stack_plots) or (stack_plots and output_names[-1] == output_name):
-                plt.savefig(
-                    f"{outputdir}/scan_{scan_var_name}_vs_{extra_str}.{save_format}",
-                    dpi=300,
+    for index, output_name in enumerate(output_names):
+        if output_name not in m_file.data:
+            print(
+                f"Warning : `{output_name}` does not exist in PROCESS dicts\n"
+                f"Warning : `{output_name}` will not be output"
+            )
+            continue
+
+        fig, ax = plt.subplots()
+        x_contour = [scan_2_var.get_val(m_file, scan=i + 1) for i in range(n_scan_2)]
+
+        if twod_contour:
+            y_contour = [
+                scan_var.get_val(m_file, scan=i + 1)
+                for i in range(1, n_scan_1 * n_scan_2, n_scan_2)
+            ]
+
+            output_contour_z = np.zeros((n_scan_1, n_scan_2))
+            for i in contour_conv_ij:
+                ind1 = (i - 1) // n_scan_2
+                ind2 = (i - 1) % n_scan_2
+                output_contour_z[ind1][(ind2 if ind1 % 2 == 0 else (-ind2 - 1))] = (
+                    m_file.get(output_name, scan=i)
                 )
-                plt.show()
-                plt.clf()
-                plt.close()
-        # ------------
+            flat_output_z = output_contour_z.flatten()
+            flat_output_z.sort()
 
-    # In case of a 2D scan
-    else:
-        # Opening the MFILE.DAT
-        m_file = MFile(filename=input_files[0])
-
-        # Number of scan points
-        n_scan_1 = int(m_file.data["isweep"].get_scan(-1))
-        n_scan_2 = int(m_file.data["isweep_2"].get_scan(-1))
-        # Selecting the converged runs only
-        contour_conv_ij = []  # List of non-converged scan point numbers
-        # 2D array of converged scan point numbers (sweep = rows, sweep_2 = columns)
-        conv_ij = []
-        ii_jj = 0
-        for ii in range(n_scan_1):
-            conv_ij.append([])
-            for _jj in range(n_scan_2):
-                ii_jj += 1  # Represents the scan point number in the MFILE
-                ifail = m_file.data["ifail"].get_scan(ii_jj)
-                if ifail == 1:
-                    conv_ij[ii].append(
-                        ii_jj
-                    )  # Only appends scan number if scan converged
-                    contour_conv_ij.append(ii_jj)
-                else:
-                    failed_value_1 = m_file.data[scan_var_name].get_scan(ii_jj)
-                    failed_value_2 = m_file.data[scan_2_var_name].get_scan(ii_jj)
-                    print(
-                        "Warning : Non-convergent scan point : "
-                        f"({scan_var_name},{scan_2_var_name}) "
-                        f"= ({failed_value_1},{failed_value_2})\n"
-                        "Warning : This point will not be shown."
-                    )
+            levels = np.linspace(
+                next(filter(lambda i: i > 0.0, flat_output_z)), flat_output_z.max(), 50
+            )
+            contour = ax.contourf(x_contour, y_contour, output_contour_z, levels=levels)
 
-        # Looping over requested outputs
-        for index, output_name in enumerate(output_names):
-            # Check if the output variable exists in the MFILE
-            if output_name not in m_file.data:
-                print(
-                    f"Warning : `{output_name}` does not exist in PROCESS dicts\n"
-                    f"Warning : `{output_name}` will not be output"
-                )
-                continue
-
-            # Declaring the outputs
-            output_arrays = []
-
-            if twod_contour:
-                output_contour_z = np.zeros((n_scan_1, n_scan_2))
-                x_contour = [
-                    m_file.data[scan_2_var_name].get_scan(i + 1) for i in range(n_scan_2)
-                ]
-                y_contour = [
-                    m_file.data[scan_var_name].get_scan(i + 1)
-                    for i in range(1, n_scan_1 * n_scan_2, n_scan_2)
-                ]
-                for i in contour_conv_ij:
-                    output_contour_z[((i - 1) // n_scan_2)][
-                        (
-                            ((i - 1) % n_scan_2)
-                            if ((i - 1) // n_scan_2) % 2 == 0
-                            else (-((i - 1) % n_scan_2) - 1)
-                        )
-                    ] = m_file.data[output_name].get_scan(i)
-
-                flat_output_z = output_contour_z.flatten()
-                flat_output_z.sort()
-                fig, ax = plt.subplots()
-                levels = np.linspace(
-                    next(filter(lambda i: i > 0.0, flat_output_z)),
-                    flat_output_z.max(),
-                    50,
-                )
-                contour = ax.contourf(
-                    x_contour,
-                    y_contour,
-                    output_contour_z,
-                    levels=levels,
-                )
+            fig.colorbar(contour).set_label(
+                label=get_label(output_name), size=y_axis.font_size
+            )
+            ax.set_ylabel(get_label(scan_var.name), fontsize=y_axis.font_size)
 
-                fig.colorbar(contour).set_label(
-                    label=(
-                        meta[output_name].latex
-                        if output_name in meta
-                        else f"{output_name}"
-                    ),
-                    size=axis_font_size,
-                )
-                plt.ylabel(
-                    (
-                        meta[scan_var_name].latex
-                        if scan_var_name in meta
-                        else f"{scan_var_name}"
-                    ),
-                    fontsize=axis_font_size,
-                )
-                plt.xlabel(
-                    (
-                        meta[scan_2_var_name].latex
-                        if scan_2_var_name in meta
-                        else f"{scan_2_var_name}"
-                    ),
-                    fontsize=axis_font_size,
-                )
-                if len(y_axis_range) > 0:
-                    y_divisions = (y_axis_range[1] - y_axis_range[0]) / 10
-                if y_axis_percent:
-                    if y_max[index] is None:
-                        y_max[index] = max(np.abs(y_contour))
-                    yticks = mtick.PercentFormatter(y_max[index])
-                    if len(y_axis_range) > 0:
-                        y_divisions = 5 * math.ceil(y_divisions / 5) * y_max[index] / 100
-                        y_range = (
-                            y_axis_range[0] * y_max[index] / 100,
-                            y_axis_range[1] * y_max[index] / 100,
-                        )
-                    ax.yaxis.set_major_formatter(yticks)
-                if len(y_axis_range) > 0:
-                    if y_axis_percent is False:
-                        y_range = y_axis_range
-                    ax.set_ylim(y_range[0], y_range[1])
-                    ax.yaxis.set_major_locator(mtick.MultipleLocator(y_divisions))
-                if len(x_axis_range) > 0:
-                    x_divisions = (x_axis_range[1] - x_axis_range[0]) / 10
-                if x_axis_percent:
-                    if x_max[index] is None:
-                        x_max[index] = max(np.abs(x_contour))
-                    xticks = mtick.PercentFormatter(x_max[index])
-                    if len(x_axis_range) > 0:
-                        x_divisions = 5 * math.ceil(x_divisions / 5) * x_max[index] / 100
-                        x_range = (
-                            x_axis_range[0] * x_max[index] / 100,
-                            x_axis_range[1] * x_max[index] / 100,
-                        )
-                    ax.xaxis.set_major_formatter(xticks)
-                if len(x_axis_range) > 0:
-                    if x_axis_percent is False:
-                        x_range = x_axis_range
-                    plt.xlim(x_range[0], x_range[1])
-                    ax.xaxis.set_major_locator(mtick.MultipleLocator(x_divisions))
-                plt.rc("xtick", labelsize=axis_tick_size)
-                plt.rc("ytick", labelsize=axis_tick_size)
-                plt.tight_layout()
-                plt.savefig(
-                    outputdir
-                    / f"scan_{output_name}_vs_{scan_var_name}_{scan_2_var_name}.{save_format}"  # noqa: E501
-                )
-                plt.grid(True)
-                plt.show()
-                plt.clf()
+        else:
+            y_contour = [m_file.get(output_name, scan=i + 1) for i in range(n_scan_2)]
 
-            else:
-                # Converged indexes, for normal 2D line plot
-                fig, ax = plt.subplots()
-                for conv_j in (
-                    conv_ij
-                ):  # conv_j is an array element containing the converged scan numbers
-                    # Scanned variables
-                    scan_1_var_array = np.zeros(len(conv_j))
-                    scan_2_var_array = np.zeros(len(conv_j))
-                    output_array = np.zeros(len(conv_j))
-                    for jj in range(len(conv_j)):
-                        scan_1_var_array[jj] = m_file.data[scan_var_name].get_scan(
-                            conv_j[jj]
-                        )
-                        scan_2_var_array[jj] = m_file.data[scan_2_var_name].get_scan(
-                            conv_j[jj]
-                        )
-                        output_array[jj] = m_file.data[output_name].get_scan(conv_j[jj])
-
-                    # Label formating
-                    labl = f"{meta[scan_var_name].latex if scan_var_name in meta else {scan_var_name}} = {scan_1_var_array[0]}"  # noqa: E501
-
-                    # Plot the graph
-                    ax.plot(scan_2_var_array, output_array, "--o", label=labl)
-
-                plt.grid(True)
-                plt.ylabel(
+            # conv_j is an array element containing the converged scan numbers
+            for conv_j in conv_ij:
+                # Scanned variables
+                scan_1_var_array, scan_2_var_array, output_array = np.array([
                     (
-                        meta[output_name].latex
-                        if output_name in meta
-                        else f"{output_name}"
-                    ),
-                    fontsize=axis_font_size,
-                )
-                plt.xlabel(
-                    (
-                        meta[scan_2_var_name].latex
-                        if scan_2_var_name in meta
-                        else f"{scan_2_var_name}"
-                    ),
-                    fontsize=axis_font_size,
-                )
-                plt.legend(loc="best", fontsize=legend_size)
-                y_data = [
-                    m_file.data[output_name].get_scan(i + 1) for i in range(n_scan_2)
-                ]
-                if len(y_axis_range) > 0:
-                    y_divisions = (y_axis_range[1] - y_axis_range[0]) / 10
-                if y_axis_percent:
-                    if y_max[index] is None:
-                        y_max[index] = max(np.abs(y_data))
-                    yticks = mtick.PercentFormatter(y_max[index])
-                    if len(y_axis_range) > 0:
-                        y_divisions = 5 * math.ceil(y_divisions / 5) * y_max[index] / 100
-                        y_range = (
-                            y_axis_range[0] * y_max[index] / 100,
-                            y_axis_range[1] * y_max[index] / 100,
-                        )
-                    ax.yaxis.set_major_formatter(yticks)
-                if len(y_axis_range) > 0:
-                    if y_axis_percent is False:
-                        y_range = y_axis_range
-                    ax.set_ylim(y_range[0], y_range[1])
-                    ax.yaxis.set_major_locator(mtick.MultipleLocator(y_divisions))
-                x_data = [
-                    m_file.data[scan_2_var_name].get_scan(i + 1) for i in range(n_scan_2)
-                ]
-                if len(x_axis_range) > 0:
-                    x_divisions = (x_axis_range[1] - x_axis_range[0]) / 10
-                if x_axis_percent:
-                    if x_max[index] is None:
-                        x_max[index] = max(np.abs(x_data))
-                    xticks = mtick.PercentFormatter(x_max[index])
-                    if len(x_axis_range) > 0:
-                        x_divisions = 5 * math.ceil(x_divisions / 5) * x_max[index] / 100
-                        x_range = (
-                            x_axis_range[0] * x_max[index] / 100,
-                            x_axis_range[1] * x_max[index] / 100,
-                        )
-                    ax.xaxis.set_major_formatter(xticks)
-                if len(x_axis_range) > 0:
-                    if x_axis_percent is False:
-                        x_range = x_axis_range
-                    plt.xlim(x_range[0], x_range[1])
-                    ax.xaxis.set_major_locator(mtick.MultipleLocator(x_divisions))
-                plt.rc("xtick", labelsize=8)
-                plt.rc("ytick", labelsize=8)
-                plt.tight_layout()
-                plt.savefig(
-                    outputdir
-                    / f"scan_{output_name}_vs_{scan_var_name}_{scan_2_var_name}.{save_format}"  # noqa: E501
-                )
+                        scan_var.get_val(m_file, scan=conv_j[jj]),
+                        scan_2_var.get_val(m_file, scan=conv_j[jj]),
+                        m_file.get(output_name, scan=conv_j[jj]),
+                    )
+                    for jj in range(len(conv_j))
+                ]).T
+
+                label = f"{get_label(scan_var.name)} = {scan_1_var_array[0]}"
+                ax.plot(scan_2_var_array, output_array, "--o", label=label)
+
+            ax.set_ylabel(get_label(output_name), fontsize=y_axis.font_size)
+            fig.legend(loc="best", fontsize=x_axis.legend_size)
+
+        ax.set_xlabel(get_label(scan_2_var.name), fontsize=x_axis.font_size)
 
-                # Display plot (used in Jupyter notebooks)
-                plt.show()
-                plt.clf()
+        axis_manipulation(ax, axis=x_axis, index=index, contour=x_contour)
+        axis_manipulation(ax, axis=y_axis, index=index, contour=y_contour)
+        ax.grid(True)
+        fname = f"scan_{output_name}_vs_{scan_var.name}_{scan_2_var.name}.{save_format}"
+        fig.savefig(outputdir / fname)
+        plt.show()
diff --git a/process/core/scan.py b/process/core/scan.py
index d77fa81df5..e74a960c0d 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -1,9 +1,12 @@
+"""Scanning mechanics"""
+
 from __future__ import annotations
 
 import logging
 import time
-from dataclasses import astuple, dataclass
+from dataclasses import dataclass, field
 from enum import Enum
+from types import DynamicClassAttribute
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -12,6 +15,7 @@
 from process.core import constants, process_output
 from process.core.caller import write_output_files
 from process.core.exceptions import ProcessValueError
+from process.core.io.data_structure_dicts import get_dicts
 from process.core.log import logging_model_handler, show_errors
 from process.core.solver import constraints
 from process.core.solver.solver_handler import SolverHandler
@@ -22,169 +26,204 @@
 if TYPE_CHECKING:
     from process.core.model import DataStructure, Model
 
+
 logger = logging.getLogger(__name__)
 
 
 @dataclass
 class ScanVariable:
-    variable_name: str
-    variable_description: str
-    variable_num: int
+    """Scan variable container"""
+
+    number: int
+    area: Area = field(repr=False)
+    _out_name_: str | None = None
+
+
+class Area(Enum):
+    """Scan variable data structure area shorthand
 
-    def __iter__(self):
-        return iter(astuple(self)[:2])
+    this mirrors data_structure.area
+    """
 
+    P = "physics"
+    D = "divertor"
+    C = "constraints"
+    T = "tfcoil"
+    TR = "rebco"
+    CD = "current_drive"
+    NUM = "numerics"
+    CST = "costs"
+    IR = "impurity_radiation"
+    B = "build"
+    HT = "heat_transport"
+    PF = "pf_coil"
+    CS = "cs_fatigue"
+    FWBS = "fwbs"
+
+
+class ScanVariables(ScanVariable, Enum):
+    """Scan variable options"""
 
-class ScanVariables(Enum):
     @classmethod
-    def _missing_(cls, var):
-        if isinstance(var, int):
+    def _missing_(cls, value):
+        if isinstance(value, int):
             for sv in cls:
-                if sv.value.variable_num == var:
+                if sv.number == value:
                     return sv
-        return super()._missing_(var)
-
-    aspect = ScanVariable("aspect", "Aspect_ratio", 1)
-    pflux_div_heat_load_max_mw = ScanVariable(
-        "pflux_div_heat_load_max_mw", "Div_heat_limit_(MW/m2)", 2
-    )
-    p_plant_electric_net_required_mw = ScanVariable(
-        "p_plant_electric_net_required_mw", "Net_electric_power_(MW)", 3
-    )
-    hfact = ScanVariable("hfact", "Confinement_H_factor", 4)
-    j_tf_coil_full_area = ScanVariable(
-        "j_tf_coil_full_area", "TF_inboard_leg_J_(MA/m2)", 5
-    )
-    pflux_fw_neutron_max_mw = ScanVariable(
-        "pflux_fw_neutron_max_mw", "Allow._wall_load_(MW/m2)", 6
-    )
-    beamfus0 = ScanVariable("beamfus0", "Beam_bkgrd_multiplier", 7)
-    temp_plasma_electron_vol_avg_kev = ScanVariable(
-        "temp_plasma_electron_vol_avg_kev", "Electron_temperature_keV", 9
-    )
-    boundu15 = ScanVariable("boundu(15)", "Volt-second_upper_bound", 10)
-    beta_norm_max = ScanVariable("beta_norm_max", "Beta_coefficient", 11)
-    f_c_plasma_bootstrap_max = ScanVariable(
-        "f_c_plasma_bootstrap_max", "Bootstrap_fraction", 12
-    )
-    boundu10 = ScanVariable("boundu(10)", "H_factor_upper_bound", 13)
-    f_j_tf_wp_critical_max = ScanVariable(
-        "f_j_tf_wp_critical_max", "TFC_Iop_/_Icrit_margin", 14
-    )
-    rmajor = ScanVariable("rmajor", "Plasma_major_radius_(m)", 16)
-    b_tf_inboard_max = ScanVariable("b_tf_inboard_max", "Max_toroidal_field_(T)", 17)
-    eta_cd_norm_hcd_primary_max = ScanVariable(
-        "eta_cd_norm_hcd_primary_max", "Maximum_CD_gamma", 18
-    )
-    boundl16 = ScanVariable("boundl(16)", "CS_thickness_lower_bound", 19)
-    t_burn_min = ScanVariable("t_burn_min", "Minimum_burn_time_(s)", 20)
-    f_t_plant_available = ScanVariable(
-        "f_t_plant_available", "Plant_availability_factor", 22
-    )
-    p_fusion_total_max_mw = ScanVariable(
-        "p_fusion_total_max_mw", "Fusion_power_limit_(MW)", 24
-    )
-    kappa = ScanVariable("kappa", "Plasma_elongation", 25)
-    triang = ScanVariable("triang", "Plasma_triangularity", 26)
-    tbrmin = ScanVariable("tbrmin", "Min_tritium_breed._ratio", 27)
-    b_plasma_toroidal_on_axis = ScanVariable(
-        "b_plasma_toroidal_on_axis", "Tor._field_on_axis_(T)", 28
-    )
-    coreradius = ScanVariable("coreradius", "Core_radius", 29)
-    f_t_alpha_energy_confinement_min = ScanVariable(
-        "f_t_alpha_energy_confinement_min", "t_alpha_confinement/taueff_lower_limit", 31
-    )
-    epsvmc = ScanVariable("epsvmc", "VMCON error tolerance", 32)
-    boundu129 = ScanVariable("boundu(129)", " Neon upper limit", 38)
-    boundu131 = ScanVariable("boundu(131)", " Argon upper limit", 39)
-    boundu135 = ScanVariable("boundu(135)", " Xenon upper limit", 40)
-    dr_blkt_outboard = ScanVariable("dr_blkt_outboard", "Outboard blanket thick.", 41)
-    f_nd_impurity_electrons9 = ScanVariable(
-        "f_nd_impurity_electrons(9)", "Argon fraction", 42
-    )
-    sig_tf_case_max = ScanVariable(
-        "sig_tf_case_max", "Allowable_stress_in_tf_coil_case_Tresca_(pa)", 44
-    )
-    temp_tf_superconductor_margin_min = ScanVariable(
-        "temp_tf_superconductor_margin_min", "Minimum_allowable_temperature_margin", 45
-    )
-    boundu152 = ScanVariable(
-        "boundu(152)", "Max allowable f_nd_plasma_separatrix_greenwald", 46
-    )
-    n_tf_wp_pancakes = ScanVariable("n_tf_wp_pancakes", "TF Coil - n_tf_wp_pancakes", 48)
-    n_tf_wp_layers = ScanVariable("n_tf_wp_layers", "TF Coil - n_tf_wp_layers", 49)
-    f_nd_impurity_electrons13 = ScanVariable(
-        "f_nd_impurity_electrons(13)", "Xenon fraction", 50
-    )
-    f_p_div_lower = ScanVariable("f_p_div_lower", "lower_divertor_power_fraction", 51)
-    rad_fraction_sol = ScanVariable("rad_fraction_sol", "SoL radiation fraction", 52)
-    boundu157 = ScanVariable("boundu(157)", "Max allowable fvssu", 53)
-    Bc2_0K = ScanVariable("Bc2(0K)", "GL_NbTi Bc2(0K)", 54)
-    dr_shld_inboard = ScanVariable("dr_shld_inboard", "Inboard neutronic shield", 55)
-    p_cryo_plant_electric_max_mw = ScanVariable(
-        "p_cryo_plant_electric_max_mw", "max allowable p_cryo_plant_electric_mw", 56
-    )
-    boundl2 = ScanVariable("boundl(2)", "b_plasma_toroidal_on_axis minimum", 57)
-    dr_fw_plasma_gap_inboard = ScanVariable(
-        "dr_fw_plasma_gap_inboard", "Inboard FW-plasma sep gap", 58
-    )
-    dr_fw_plasma_gap_outboard = ScanVariable(
-        "dr_fw_plasma_gap_outboard", "Outboard FW-plasma sep gap", 59
-    )
-    sig_tf_wp_max = ScanVariable(
-        "sig_tf_wp_max", "Allowable_stress_in_tf_coil_conduit_Tresca_(pa)", 60
-    )
-    copperaoh_m2_max = ScanVariable(
-        "copperaoh_m2_max", "Max CS coil current / copper area", 61
-    )
-    coheof = ScanVariable("coheof", "CS coil current density at EOF (A/m2)", 62)
-    dr_cs = ScanVariable("dr_cs", "CS coil thickness (m)", 63)
-    ohhghf = ScanVariable("ohhghf", "CS height (m)", 64)
-    n_cycle_min = ScanVariable("n_cycle_min", "CS stress cycles min", 65)
-    oh_steel_frac = ScanVariable("oh_steel_frac", "CS steel fraction", 66)
-    t_crack_vertical = ScanVariable(
-        "t_crack_vertical", "Initial crack vertical size (m)", 67
-    )
-    inlet_temp_liq = ScanVariable(
-        "inlet_temp_liq", "Inlet Temperature Liquid Metal Breeder/Coolant (K)", 68
-    )
-    outlet_temp_liq = ScanVariable(
-        "outlet_temp_liq", "Outlet Temperature Liquid Metal Breeder/Coolant (K)", 69
-    )
-    blpressure_liq = ScanVariable(
-        "blpressure_liq", "Blanket liquid metal breeder/coolant pressure (Pa)", 70
-    )
-    n_liq_recirc = ScanVariable(
-        "n_liq_recirc",
-        "Selected number of liquid metal breeder recirculations per day",
-        71,
-    )
-    bz_channel_conduct_liq = ScanVariable(
-        "bz_channel_conduct_liq",
-        "Conductance of liquid metal breeder duct walls (A V-1 m-1)",
-        72,
-    )
-    pnuc_fw_ratio_dcll = ScanVariable(
-        "pnuc_fw_ratio_dcll",
-        "Ratio of FW nuclear power as fraction of total (FW+BB)",
-        73,
-    )
-    f_nuc_pow_bz_struct = ScanVariable(
-        "f_nuc_pow_bz_struct",
-        "Fraction of BZ power cooled by primary coolant for dual-coolant blanket",
-        74,
-    )
-    dx_fw_module = ScanVariable(
-        "dx_fw_module", "dx_fw_module of first wall cooling channels (m)", 75
-    )
-    eta_turbine = ScanVariable("eta_turbine", "Thermal conversion eff.", 76)
-    startupratio = ScanVariable("startupratio", "Gyrotron redundancy", 77)
-    fkind = ScanVariable("fkind", "Multiplier for Nth of a kind costs", 78)
-    eta_ecrh_injector_wall_plug = ScanVariable(
-        "eta_ecrh_injector_wall_plug", "ECH wall plug to injector efficiency", 79
-    )
-    fcoolcp = ScanVariable("fcoolcp", "Coolant fraction of TF", 80)
-    n_tf_coil_turns = ScanVariable("n_tf_coil_turns", "Number of turns in TF", 81)
+            raise ProcessValueError("Illegal scan variable number", nwp=value)
+
+        if isinstance(value, str):
+            return cls[value.replace("(", "__").replace(")", "")]
+
+        return super()._missing_(value)
+
+    @DynamicClassAttribute
+    def fname(self):
+        """Full name"""
+        if "__" in self.name:
+            return self.name.replace("__", "(") + ")"
+        return self.name
+
+    @DynamicClassAttribute
+    def out_name(self):
+        """Output name"""
+        if self._out_name_ is None:
+            return self.fname
+        return self._out_name_
+
+    def set(self, data: DataStructure, sweep_val: float):
+        """Set value of scan variable
+
+        Raises
+        ------
+        ProcessValueError
+            Scanning f_t_plant_available if i_plant_availability=1"
+
+        """
+        var_area = getattr(data, self.area.value)
+
+        if (
+            self.number == 22
+            and AvailabilityModel(var_area.i_plant_availability)
+            != AvailabilityModel.USER_INPUT
+        ):
+            raise ProcessValueError(
+                "Do not scan f_t_plant_available if i_plant_availability=1"
+            )
+
+        if "__" in self.name:
+            name, index = self.name.split("__")
+            getattr(var_area, name)[int(index) - 1] = sweep_val
+            if name == "f_nd_impurity_electrons":
+                var_area.f_nd_impurity_electron_array[int(index - 1)] = sweep_val
+        else:
+            setattr(var_area, self.name, sweep_val)
+            name = self.name
+
+        self._data_ = getattr(var_area, name)
+        self._description_ = get_dicts()["DICT_DESCRIPTIONS"][name]
+
+    @DynamicClassAttribute
+    def data(self):
+        """Variable value
+
+        Raises
+        ------
+        ValueError
+            data not set
+        """
+        if hasattr(self, "_data_"):
+            return self._data_
+        raise ValueError("Data not available")
+
+    @DynamicClassAttribute
+    def description(self):
+        """Variable description
+
+        Raises
+        ------
+        ValueError
+            description not set
+        """
+        if hasattr(self, "_description_"):
+            return self._description_
+        raise ValueError("Description not available")
+
+    def get_val(self, mfile, scan):
+        """Get value from mfile"""
+        # TODO this will fail for boundu/l we should write the scan variable to
+        # the mfile and use that directly (also replacign output names)
+        return mfile.get(self.out_name, scan=scan)
+
+    aspect = (1, Area.P)
+    pflux_div_heat_load_max_mw = (2, Area.D)
+    p_plant_electric_net_required_mw = (3, Area.C)
+    hfact = (4, Area.P)
+    j_tf_coil_full_area = (5, Area.T)
+    pflux_fw_neutron_max_mw = (6, Area.C)
+    beamfus0 = (7, Area.P)
+    temp_plasma_electron_vol_avg_kev = (9, Area.P)
+    boundu__15 = (10, Area.NUM)
+    beta_norm_max = (11, Area.P)
+    f_c_plasma_bootstrap_max = (12, Area.CD)
+    boundu__10 = (13, Area.NUM)
+    f_j_tf_wp_critical_max = (14, Area.C)  # TODO is this needed
+    rmajor = (16, Area.P)
+    b_tf_inboard_max = (17, Area.C, "b_tf_inboard_peak_symmetric")
+    eta_cd_norm_hcd_primary_max = (18, Area.C)
+    boundl__16 = (19, Area.NUM)
+    t_burn_min = (20, Area.C)
+    f_t_plant_available = (22, Area.CST)
+    p_fusion_total_max_mw = (24, Area.C)
+    kappa = (25, Area.P)
+    triang = (26, Area.P)
+    tbrmin = (27, Area.C)
+    b_plasma_toroidal_on_axis = (28, Area.P)
+    coreradius = (29, Area.IR)
+    f_alpha_energy_confinement_min = (31, Area.C)
+    epsvmc = (32, Area.NUM)
+    boundu__129 = (38, Area.NUM)
+    boundu__131 = (39, Area.NUM)
+    boundu__135 = (40, Area.NUM)
+    dr_blkt_outboard = (41, Area.B)
+    f_nd_impurity_electrons__9 = (42, Area.IR)
+    sig_tf_case_max = (44, Area.T)
+    temp_tf_superconductor_margin_min = (45, Area.T)
+    boundu__152 = (46, Area.NUM)
+    n_tf_wp_pancakes = (48, Area.T)
+    n_tf_wp_layers = (49, Area.T)
+    f_nd_impurity_electrons__13 = (50, Area.IR)
+    f_p_div_lower = (51, Area.P)
+    rad_fraction_sol = (52, Area.P)
+    boundu__157 = (53, Area.NUM)
+    b_crit_upper_nbti = (54, Area.T)
+    dr_shld_inboard = (55, Area.B)
+    p_cryo_plant_electric_max_mw = (56, Area.HT)
+    boundl__2 = (57, Area.NUM)
+    dr_fw_plasma_gap_inboard = (58, Area.B)
+    dr_fw_plasma_gap_outboard = (59, Area.B)
+    sig_tf_wp_max = (60, Area.T)
+    copperaoh_m2_max = (61, Area.TR)
+    coheof = (62, Area.PF)
+    dr_cs = (63, Area.B)
+    ohhghf = (64, Area.PF)
+    n_cycle_min = (65, Area.CS)
+    oh_steel_frac = (66, Area.PF)
+    t_crack_vertical = (67, Area.CS)
+    inlet_temp_liq = (68, Area.FWBS)
+    outlet_temp_liq = (69, Area.FWBS)
+    blpressure_liq = (70, Area.FWBS)
+    n_liq_recirc = (71, Area.FWBS)
+    bz_channel_conduct_liq = (72, Area.FWBS)
+    pnuc_fw_ratio_dcll = (73, Area.FWBS)
+    f_nuc_pow_bz_struct = (74, Area.FWBS)
+    dx_fw_module = (75, Area.FWBS)
+    eta_turbine = (76, Area.HT)
+    startupratio = (77, Area.CST)
+    fkind = (78, Area.CST)
+    eta_ecrh_injector_wall_plug = (79, Area.CD)
+    fcoolcp = (80, Area.T)
+    n_tf_coil_turns = (81, Area.T)
 
 
 class Scan:
@@ -215,6 +254,11 @@ def run_scan(self):
         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.
+
+        Raises
+        ------
+        ProcessValueError
+            isweep value greater than IPNSCNS
         """
         if self.data.scan.isweep == 0:
             # Solve single problem, rather than an array of problems (scan)
@@ -877,7 +921,7 @@ def scan_1d(self):
         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_name, _ = self.scan_select(
+        nsweep_var = self.scan_select(
             self.data.scan.nsweep, self.data.scan.sweep, self.data.scan.isweep
         )
         converged_count = 0
@@ -891,13 +935,13 @@ def scan_1d(self):
             if scan_1d_ifail_dict[iscan] == 1:
                 converged_count += 1
                 print(
-                    f"Scan {iscan:02d}: {nsweep_var_name} = {sweep_values[iscan - 1]} "
+                    f"Scan {iscan:02d}: {nsweep_var.fname} = {sweep_values[iscan - 1]} "
                     + " " * offsets[iscan - 1]
                     + "\u001b[32mCONVERGED \u001b[0m"
                 )
             else:
                 print(
-                    f"Scan {iscan:02d}: {nsweep_var_name} = {sweep_values[iscan - 1]} "
+                    f"Scan {iscan:02d}: {nsweep_var.fname} = {sweep_values[iscan - 1]} "
                     + " " * offsets[iscan - 1]
                     + "\u001b[31mUNCONVERGED \u001b[0m"
                 )
@@ -936,10 +980,10 @@ def scan_2d(self):
         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_name, _ = self.scan_select(
+        nsweep_var = self.scan_select(
             self.data.scan.nsweep, self.data.scan.sweep, self.data.scan.isweep
         )
-        nsweep_2_var_name, _ = self.scan_select(
+        nsweep_2_var = self.scan_select(
             self.data.scan.nsweep_2, self.data.scan.sweep_2, self.data.scan.isweep_2
         )
         converged_count = 0
@@ -965,8 +1009,8 @@ def scan_2d(self):
                     converged_count += 1
                     print(
                         (
-                            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]
@@ -976,8 +1020,8 @@ def scan_2d(self):
                 else:
                     print(
                         (
-                            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]
@@ -991,6 +1035,7 @@ def scan_2d(self):
 
     @staticmethod
     def scan_2d_init(scan_data: ScanData):
+        """Scan 2d initialisation"""
         process_output.ovarre(
             constants.MFILE,
             "Number of first variable scan points",
@@ -1029,10 +1074,12 @@ def scan_2d_init(scan_data: ScanData):
         )
 
     def scan_1d_write_point_header(self, iscan: int):
+        """Scan 1d header"""
         self.data.globals.iscan_global = iscan
-        self.data.globals.vlabel, self.data.globals.xlabel = self.scan_select(
-            self.data.scan.nsweep, self.data.scan.sweep, 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
 
         process_output.oblnkl(constants.NOUT)
         process_output.ostars(constants.NOUT, 110)
@@ -1055,17 +1102,20 @@ def scan_1d_write_point_header(self, iscan: int):
         )
 
     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
+        sv_1 = self.scan_select(self.data.scan.nsweep, self.data.scan.sweep, iscan_1)
 
-        self.data.globals.vlabel, self.data.globals.xlabel = self.scan_select(
-            self.data.scan.nsweep, self.data.scan.sweep, iscan_1
-        )
-        self.data.globals.vlabel_2, self.data.globals.xlabel_2 = self.scan_select(
-            self.data.scan.nsweep_2, self.data.scan.sweep_2, iscan_r
-        )
+        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
 
         process_output.oblnkl(constants.NOUT)
         process_output.ostars(constants.NOUT, 110)
@@ -1093,6 +1143,7 @@ def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2):
 
     @staticmethod
     def scan_1d_write_plot(scan_data: ScanData):
+        """Scan 1d plotter"""
         if scan_data.first_call_1d:
             process_output.ovarre(
                 constants.MFILE,
@@ -1109,159 +1160,8 @@ def scan_1d_write_plot(scan_data: ScanData):
 
             scan_data.first_call_1d = False
 
-    def scan_select(self, nwp, swp, iscn):
-        match nwp:
-            case 1:
-                self.data.physics.aspect = swp[iscn - 1]
-            case 2:
-                self.data.divertor.pflux_div_heat_load_max_mw = swp[iscn - 1]
-            case 3:
-                self.data.constraints.p_plant_electric_net_required_mw = swp[iscn - 1]
-            case 4:
-                self.data.physics.hfact = swp[iscn - 1]
-            case 5:
-                self.data.tfcoil.j_tf_coil_full_area = swp[iscn - 1]
-            case 6:
-                self.data.constraints.pflux_fw_neutron_max_mw = swp[iscn - 1]
-            case 7:
-                self.data.physics.beamfus0 = swp[iscn - 1]
-            case 9:
-                self.data.physics.temp_plasma_electron_vol_avg_kev = swp[iscn - 1]
-            case 10:
-                self.data.numerics.boundu[14] = swp[iscn - 1]
-            case 11:
-                self.data.physics.beta_norm_max = swp[iscn - 1]
-            case 12:
-                self.data.current_drive.f_c_plasma_bootstrap_max = swp[iscn - 1]
-            case 13:
-                self.data.numerics.boundu[9] = swp[iscn - 1]
-            case 16:
-                self.data.physics.rmajor = swp[iscn - 1]
-            case 17:
-                self.data.constraints.b_tf_inboard_max = swp[iscn - 1]
-            case 18:
-                self.data.constraints.eta_cd_norm_hcd_primary_max = swp[iscn - 1]
-            case 19:
-                self.data.numerics.boundl[15] = swp[iscn - 1]
-            case 20:
-                self.data.constraints.t_burn_min = swp[iscn - 1]
-            case 22:
-                if (
-                    AvailabilityModel(self.data.costs.i_plant_availability)
-                    != AvailabilityModel.USER_INPUT
-                ):
-                    raise ProcessValueError(
-                        "Do not scan f_t_plant_available if i_plant_availability is not "
-                        "set to 1 (user input) in the input file."
-                    )
-                self.data.costs.f_t_plant_available = swp[iscn - 1]
-            case 24:
-                self.data.constraints.p_fusion_total_max_mw = swp[iscn - 1]
-            case 25:
-                self.data.physics.kappa = swp[iscn - 1]
-            case 26:
-                self.data.physics.triang = swp[iscn - 1]
-            case 27:
-                self.data.constraints.tbrmin = swp[iscn - 1]
-            case 28:
-                self.data.physics.b_plasma_toroidal_on_axis = swp[iscn - 1]
-            case 29:
-                self.data.impurity_radiation.radius_plasma_core_norm = swp[iscn - 1]
-            case 31:
-                self.data.constraints.f_t_alpha_energy_confinement_min = swp[iscn - 1]
-            case 32:
-                self.data.numerics.epsvmc = swp[iscn - 1]
-            case 38:
-                self.data.numerics.boundu[128] = swp[iscn - 1]
-            case 39:
-                self.data.numerics.boundu[130] = swp[iscn - 1]
-            case 40:
-                self.data.numerics.boundu[134] = swp[iscn - 1]
-            case 41:
-                self.data.build.dr_blkt_outboard = swp[iscn - 1]
-            case 42:
-                self.data.impurity_radiation.f_nd_impurity_electrons[8] = swp[iscn - 1]
-                self.data.impurity_radiation.f_nd_impurity_electron_array[8] = (
-                    self.data.impurity_radiation.f_nd_impurity_electrons[8]
-                )
-            case 44:
-                self.data.tfcoil.sig_tf_case_max = swp[iscn - 1]
-            case 45:
-                self.data.tfcoil.temp_tf_superconductor_margin_min = swp[iscn - 1]
-            case 46:
-                self.data.numerics.boundu[151] = swp[iscn - 1]
-            case 48:
-                self.data.tfcoil.n_tf_wp_pancakes = int(swp[iscn - 1])
-            case 49:
-                self.data.tfcoil.n_tf_wp_layers = int(swp[iscn - 1])
-            case 50:
-                self.data.impurity_radiation.f_nd_impurity_electrons[12] = swp[iscn - 1]
-                self.data.impurity_radiation.f_nd_impurity_electron_array[12] = (
-                    self.data.impurity_radiation.f_nd_impurity_electrons[12]
-                )
-            case 51:
-                self.data.physics.f_p_div_lower = swp[iscn - 1]
-            case 52:
-                self.data.physics.rad_fraction_sol = swp[iscn - 1]
-            case 53:
-                self.data.numerics.boundu[156] = swp[iscn - 1]
-            case 54:
-                self.data.tfcoil.b_crit_upper_nbti = swp[iscn - 1]
-            case 55:
-                self.data.build.dr_shld_inboard = swp[iscn - 1]
-            case 56:
-                self.data.heat_transport.p_cryo_plant_electric_max_mw = swp[iscn - 1]
-            case 57:
-                self.data.numerics.boundl[1] = swp[iscn - 1]
-            case 58:
-                self.data.build.dr_fw_plasma_gap_inboard = swp[iscn - 1]
-            case 59:
-                self.data.build.dr_fw_plasma_gap_outboard = swp[iscn - 1]
-            case 60:
-                self.data.tfcoil.sig_tf_wp_max = swp[iscn - 1]
-            case 61:
-                self.data.rebco.copperaoh_m2_max = swp[iscn - 1]
-            case 62:
-                self.data.pf_coil.j_cs_flat_top_end = swp[iscn - 1]
-            case 63:
-                self.data.build.dr_cs = swp[iscn - 1]
-            case 64:
-                self.data.pf_coil.f_z_cs_tf_internal = swp[iscn - 1]
-            case 65:
-                self.data.cs_fatigue.n_cycle_min = swp[iscn - 1]
-            case 66:
-                self.data.pf_coil.f_a_cs_turn_steel = swp[iscn - 1]
-            case 67:
-                self.data.cs_fatigue.t_crack_vertical = swp[iscn - 1]
-            case 68:
-                self.data.fwbs.inlet_temp_liq = swp[iscn - 1]
-            case 69:
-                self.data.fwbs.outlet_temp_liq = swp[iscn - 1]
-            case 70:
-                self.data.fwbs.blpressure_liq = swp[iscn - 1]
-            case 71:
-                self.data.fwbs.n_liq_recirc = swp[iscn - 1]
-            case 72:
-                self.data.fwbs.bz_channel_conduct_liq = swp[iscn - 1]
-            case 73:
-                self.data.fwbs.pnuc_fw_ratio_dcll = swp[iscn - 1]
-            case 74:
-                self.data.fwbs.f_nuc_pow_bz_struct = swp[iscn - 1]
-            case 75:
-                self.data.fwbs.dx_fw_module = swp[iscn - 1]
-            case 76:
-                self.data.heat_transport.eta_turbine = swp[iscn - 1]
-            case 77:
-                self.data.costs.startupratio = swp[iscn - 1]
-            case 78:
-                self.data.costs.fkind = swp[iscn - 1]
-            case 79:
-                self.data.current_drive.eta_ecrh_injector_wall_plug = swp[iscn - 1]
-            case 80:
-                self.data.tfcoil.fcoolcp = swp[iscn - 1]
-            case 81:
-                self.data.tfcoil.n_tf_coil_turns = swp[iscn - 1]
-            case _:
-                raise ProcessValueError("Illegal scan variable number", nwp=nwp)
-
-        return ScanVariables(int(nwp)).value
+    def scan_select(self, nsweep, sweep, iscan) -> ScanVariables:
+        """Select a scan"""
+        sv = ScanVariables(nsweep)
+        sv.set(self.data, sweep[iscan - 1])
+        return sv
diff --git a/process/data_structure/scan_variables.py b/process/data_structure/scan_variables.py
index 04aeb2dc9f..5e0180b1b6 100644
--- a/process/data_structure/scan_variables.py
+++ b/process/data_structure/scan_variables.py
@@ -32,85 +32,7 @@ class ScanData:
     """Number of 2D scan points to calculate"""
 
     nsweep: int = 1
-    """Switch denoting quantity to scan:<UL>
-    <LI> 1  aspect
-    <LI> 2  pflux_div_heat_load_max_mw
-    <LI> 3  p_plant_electric_net_required_mw
-    <LI> 4  hfact
-    <LI> 5  j_tf_coil_full_area
-    <LI> 6  pflux_fw_neutron_max_mw
-    <LI> 7  beamfus0
-    <LI> 8  NOT USED
-    <LI> 9  temp_plasma_electron_vol_avg_kev
-    <LI> 10 NOT USED
-    <LI> 11 beta_norm_max
-    <LI> 12 f_c_plasma_bootstrap_max
-    <LI> 13 boundu(10: hfact)
-    <LI> 14 f_j_tf_wp_critical_max
-    <LI> 16 rmajor
-    <LI> 15 NOT USED
-    <LI> 17 b_tf_inboard_max
-    <LI> 18 eta_cd_norm_hcd_primary_max
-    <LI> 19 boundl(16: dr_cs)
-    <LI> 20 t_burn_min
-    <LI> 21 NOT USED
-    <LI> 22 f_t_plant_available (N.B. requires i_plant_availability=0)
-    <LI> 23 NOT USED
-    <LI> 24 p_fusion_total_max_mw
-    <LI> 25 kappa
-    <LI> 26 triang
-    <LI> 27 tbrmin (for blktmodel > 0 only)
-    <LI> 28 b_plasma_toroidal_on_axis
-    <LI> 29 radius_plasma_core_norm
-    <LI> 30 fimpvar # OBSOLETE
-    <LI> 31 f_t_alpha_energy_confinement_min
-    <LI> 32 epsvmc
-    <LI> 33 ttarget
-    <LI> 34 qtargettotal
-    <LI> 35 lambda_q_omp
-    <LI> 36 lambda_target
-    <LI> 37 lcon_factor
-    <LI> 38 Neon upper limit
-    <LI> 39 Argon upper limit
-    <LI> 40 Xenon upper limit
-    <LI> 41 dr_blkt_outboard
-    <LI> 42 Argon fraction f_nd_impurity_electrons(9)
-    <LI> 43 normalised minor radius at which electron cyclotron current drive is maximum
-    <LI> 44 Allowable maximum shear stress (Tresca) in tf coil structural material
-    <LI> 45 Minimum allowable temperature margin ; tf coils
-    <LI> 46 boundu(150) f_nd_plasma_separatrix_greenwald
-    <LI> 47 impurity_enrichment(9) Argon impurity enrichment
-    <LI> 48 TF coil - n_tf_wp_pancakes (integer turn winding pack)
-    <LI> 49 TF coil - n_tf_wp_layers (integer turn winding pack)
-    <LI> 50 Xenon fraction f_nd_impurity_electrons(13)
-    <LI> 51 Power fraction to lower DN Divertor f_p_div_lower
-    <LI> 52 SoL radiation fraction
-    <LI> 54 GL_nbti upper critical field at 0 Kelvin
-    <LI> 55 `dr_shld_inboard` : Inboard neutron shield thickness
-    <LI> 56 p_cryo_plant_electric_max_mw: Maximum cryogenic power (ixx=164, ixc=87)
-    <LI> 57 `b_plasma_toroidal_on_axis` lower boundary
-    <LI> 58 `dr_fw_plasma_gap_inboard` : Inboard plasma-first wall gap
-    <LI> 59 `dr_fw_plasma_gap_outboard` : Outboard plasma-first wall gap
-    <LI> 60 sig_tf_wp_max: Allowable stress in TF Coil conduit (Tresca)
-    <LI> 61 copperaoh_m2_max : CS coil current / copper area
-    <LI> 62 j_cs_flat_top_end : CS coil current density at EOF
-    <LI> 63 dr_cs : CS thickness (m)
-    <LI> 64 f_z_cs_tf_internal : CS height (m)
-    <LI> 65 n_cycle_min : Minimum cycles for CS stress model constraint 90
-    <LI> 66 f_a_cs_turn_steel: Steel fraction in CS coil
-    <LI> 67 t_crack_vertical: Initial crack vertical dimension (m) </UL>
-    <LI> 68 `inlet_temp_liq' : Inlet temperature of blanket liquid metal coolant/breeder (K)
-    <LI> 69 `outlet_temp_liq' : Outlet temperature of blanket liquid metal coolant/breeder (K)
-    <LI> 70 `blpressure_liq' : Blanket liquid metal breeder/coolant pressure (Pa)
-    <LI> 71 `n_liq_recirc' : Selected number of liquid metal breeder recirculations per day
-    <LI> 72 `bz_channel_conduct_liq' : Conductance of liquid metal breeder duct walls (A V-1 m-1)
-    <LI> 73 `pnuc_fw_ratio_dcll' : Ratio of FW nuclear power as fraction of total (FW+BB)
-    <LI> 74 `f_nuc_pow_bz_struct' : Fraction of BZ power cooled by primary coolant for dual-coolant blanket
-    <LI> 75 dx_fw_module : pitch of first wall cooling channels (m)
-    <LI> 76 eta_turbine : Thermal conversion eff.
-    <LI> 77 startupratio : Gyrotron redundancy
-    <LI> 78 fkind : Multiplier for Nth of a kind costs
-    <LI> 79 eta_ecrh_injector_wall_plug : ECH wall plug to injector efficiency
+    """Switch denoting quantity to scan
     """
 
     nsweep_2: int = 3