diff --git a/examples/single_model_evaluation.ex.py b/examples/single_model_evaluation.ex.py
index af37205f95..7497cd7b88 100644
--- a/examples/single_model_evaluation.ex.py
+++ b/examples/single_model_evaluation.ex.py
@@ -51,7 +51,7 @@
 # Print initial values of interest
 def print_values():
     print(
-        f"W frac = {data_structure.impurity_radiation_module.f_nd_impurity_electron_array[13]:.3e}"
+        f"W frac = {single_run.data.impurity_radiation.f_nd_impurity_electron_array[13]:.3e}"
     )
     print(f"p_plasma_rad_mw = {data_structure.physics_variables.p_plasma_rad_mw:.3e}")
     print(
@@ -65,7 +65,7 @@ def print_values():
 # Now try increasing the tungsten impurity fraction to see if there's a change in the divertor power.
 
 # %%
-data_structure.impurity_radiation_module.f_nd_impurity_electron_array[13] = 5.0e-5
+single_run.data.impurity_radiation.f_nd_impurity_electron_array[13] = 5.0e-5
 single_run.models.physics.run()
 print_values()
 
@@ -90,9 +90,7 @@ def run_impurities(w_imp_fracs):
     # Loop over W impurity values, evaluate model and store responses at each point
     for i, imp_frac in enumerate(w_imp_fracs):
         # Set W impurity fraction, then run physics model
-        data_structure.impurity_radiation_module.f_nd_impurity_electron_array[13] = (
-            imp_frac
-        )
+        single_run.data.impurity_radiation.f_nd_impurity_electron_array[13] = imp_frac
         single_run.models.physics.run()
 
         # Evaluate constraint equation 15 (L-H threshold constraint)
diff --git a/process/core/init.py b/process/core/init.py
index d5ed1f892e..3ed9408a30 100644
--- a/process/core/init.py
+++ b/process/core/init.py
@@ -16,9 +16,7 @@
 from process.core.log import logging_model_handler
 from process.core.solver import iteration_variables
 from process.core.solver.constraints import ConstraintManager
-from process.data_structure.impurity_radiation_module import (
-    init_impurity_radiation_module,
-)
+from process.data_structure.impurity_radiation_variables import N_IMPURITIES
 from process.data_structure.physics_variables import (
     init_physics_module,
     init_physics_variables,
@@ -241,7 +239,6 @@ def init_all_module_vars():
     logging_model_handler.clear_logs()
     data_structure.numerics.init_numerics()
     data_structure.global_variables.init_global_variables()
-    init_impurity_radiation_module()
     init_physics_module()
     init_physics_variables()
     init_scan_variables()
@@ -355,16 +352,16 @@ def check_process(inputs, data):  # noqa: ARG001
         data.buildings.triv = 0.0
         data.heat_transport.p_tritium_plant_electric_mw = 0.0
 
-    if data_structure.impurity_radiation_module.f_nd_impurity_electrons[1] != 0.1:  # noqa: RUF069
+    if data.impurity_radiation.f_nd_impurity_electrons[1] != 0.1:  # noqa: RUF069
         raise ProcessValidationError(
             "The thermal alpha/electron density ratio should be controlled using f_nd_alpha_electron (itv 109) and not f_nd_impurity_electrons(2)."
             "f_nd_impurity_electrons(2) should be removed from the input file, or set to the default value 0.1D0."
         )
 
     # Impurity fractions
-    for imp in range(data_structure.impurity_radiation_module.N_IMPURITIES):
-        data_structure.impurity_radiation_module.f_nd_impurity_electron_array[imp] = (
-            data_structure.impurity_radiation_module.f_nd_impurity_electrons[imp]
+    for imp in range(N_IMPURITIES):
+        data.impurity_radiation.f_nd_impurity_electron_array[imp] = (
+            data.impurity_radiation.f_nd_impurity_electrons[imp]
         )
 
     # Stop the run if oacdcp is used as an optimisation variable
diff --git a/process/core/input.py b/process/core/input.py
index 565a130876..679996af66 100644
--- a/process/core/input.py
+++ b/process/core/input.py
@@ -16,6 +16,7 @@
     ProcessValueError,
 )
 from process.core.solver.constraints import ConstraintManager
+from process.data_structure.impurity_radiation_variables import N_IMPURITIES
 from process.data_structure.pfcoil_variables import N_PF_GROUPS_MAX
 
 if TYPE_CHECKING:
@@ -143,10 +144,10 @@ def __post_init__(self):
         data_structure.physics_variables, float, range=(0.0, 1.0)
     ),
     "radius_plasma_core_norm": InputVariable(
-        data_structure.impurity_radiation_module, float, range=(0.0, 1.0)
+        "impurity_radiation", float, range=(0.0, 1.0)
     ),
     "f_p_plasma_core_rad_reduction": InputVariable(
-        data_structure.impurity_radiation_module, float, range=(0.0, 1.0)
+        "impurity_radiation", float, range=(0.0, 1.0)
     ),
     "c_beta": InputVariable(data_structure.physics_variables, float, range=(0.0, 1.0)),
     "csawth": InputVariable(data_structure.physics_variables, float, range=(0.0, 10.0)),
@@ -185,9 +186,7 @@ def __post_init__(self):
         data_structure.physics_variables, float, range=(-1.0, 5.0)
     ),
     # TODO: does f_nd_impurity_electrons require an additional range?
-    "f_nd_impurity_electrons": InputVariable(
-        data_structure.impurity_radiation_module, float, array=True
-    ),
+    "f_nd_impurity_electrons": InputVariable("impurity_radiation", float, array=True),
     "fkzohm": InputVariable(data_structure.physics_variables, float, range=(0.5, 2.0)),
     "abktflnc": InputVariable("costs", float, range=(0.1, 100.0)),
     "adivflnc": InputVariable("costs", float, range=(0.1, 100.0)),
@@ -1379,7 +1378,7 @@ def __post_init__(self):
     "impvardiv": InputVariable(
         "reinke",
         int,
-        choices=range(3, data_structure.impurity_radiation_module.N_IMPURITIES + 1),
+        choices=range(3, N_IMPURITIES + 1),
     ),
     "j_pf_coil_wp_peak": InputVariable("pf_coil", float, array=True),
     "ixc": InputVariable(
diff --git a/process/core/io/mfile/base.py b/process/core/io/mfile/base.py
index 1d27bd6eee..f43f92306d 100644
--- a/process/core/io/mfile/base.py
+++ b/process/core/io/mfile/base.py
@@ -597,7 +597,7 @@ def get_mfile_initial_ixc_values(file_path: Path, data: DataStructure):
     """
     from process.main import SingleRun  # noqa:PLC0415
 
-    SingleRun(file_path.as_posix())
+    SingleRun(file_path.as_posix(), data_structure=data)
     iteration_variables.load_iteration_variables(data)
 
     iteration_variable_names = []
diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index ac1e0951d5..a0180e00c4 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -21,7 +21,7 @@
 from process.core import constants
 from process.core.io.mfile import MFile, MFileErrorClass
 from process.core.solver.objectives import OBJECTIVE_NAMES
-from process.data_structure import impurity_radiation_module
+from process.data_structure.impurity_radiation_variables import N_IMPURITIES
 from process.data_structure.pfcoil_variables import NFIXMX
 from process.models.build import Build
 from process.models.geometry.blanket import (
@@ -12878,7 +12878,7 @@ def plot_ion_charge_profile(axis: plt.Axes, mfile: MFile, scan: int):
 
     n_charge_plasma_profile = []
     avg_ionisation_percentages = []
-    for imp in range(impurity_radiation_module.N_IMPURITIES):
+    for imp in range(N_IMPURITIES):
         if imp_frac[imp] > 1.0e-30:
             profile = [
                 mfile.get(f"n_charge_plasma_profile{imp}_{i}", scan=scan)
diff --git a/process/core/model.py b/process/core/model.py
index e2cdb48397..1277717e7d 100644
--- a/process/core/model.py
+++ b/process/core/model.py
@@ -16,6 +16,7 @@
 from process.data_structure.fwbs_variables import FWBSData
 from process.data_structure.heat_transport_variables import HeatTransportData
 from process.data_structure.ife_variables import IFEData
+from process.data_structure.impurity_radiation_variables import ImpurityRadiationData
 from process.data_structure.neoclassics_variables import NeoclassicsData
 from process.data_structure.pf_power_variables import PFPowerData
 from process.data_structure.pfcoil_variables import PFCoilData
@@ -63,6 +64,7 @@ class DataStructure:
     stellarator_config: StellaratorConfigData = initialise_later
     pf_power: PFPowerData = initialise_later
     neoclassics: NeoclassicsData = initialise_later
+    impurity_radiation: ImpurityRadiationData = initialise_later
 
     def __post_init__(self):
         for f in fields(self):
diff --git a/process/core/scan.py b/process/core/scan.py
index a6327b6277..87cf5cff6e 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -17,7 +17,6 @@
 from process.core.solver.solver_handler import SolverHandler
 from process.data_structure import (
     global_variables,
-    impurity_radiation_module,
     numerics,
     physics_variables,
     rebco_variables,
@@ -1119,7 +1118,7 @@ def scan_select(self, nwp, swp, iscn):
             case 28:
                 physics_variables.b_plasma_toroidal_on_axis = swp[iscn - 1]
             case 29:
-                impurity_radiation_module.radius_plasma_core_norm = swp[iscn - 1]
+                self.data.impurity_radiation.radius_plasma_core_norm = swp[iscn - 1]
             case 31:
                 self.data.constraints.f_alpha_energy_confinement_min = swp[iscn - 1]
             case 32:
@@ -1133,9 +1132,9 @@ def scan_select(self, nwp, swp, iscn):
             case 41:
                 self.data.build.dr_blkt_outboard = swp[iscn - 1]
             case 42:
-                impurity_radiation_module.f_nd_impurity_electrons[8] = swp[iscn - 1]
-                impurity_radiation_module.f_nd_impurity_electron_array[8] = (
-                    impurity_radiation_module.f_nd_impurity_electrons[8]
+                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:
                 tfcoil_variables.sig_tf_case_max = swp[iscn - 1]
@@ -1148,9 +1147,9 @@ def scan_select(self, nwp, swp, iscn):
             case 49:
                 tfcoil_variables.n_tf_wp_layers = int(swp[iscn - 1])
             case 50:
-                impurity_radiation_module.f_nd_impurity_electrons[12] = swp[iscn - 1]
-                impurity_radiation_module.f_nd_impurity_electron_array[12] = (
-                    impurity_radiation_module.f_nd_impurity_electrons[12]
+                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:
                 physics_variables.f_p_div_lower = swp[iscn - 1]
diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py
index 17392b4390..9a7169eb49 100644
--- a/process/core/solver/iteration_variables.py
+++ b/process/core/solver/iteration_variables.py
@@ -125,7 +125,7 @@ class IterationVariable:
     122: IterationVariable("f_a_cs_turn_steel", "pf_coil", 0.001, 0.950),
     125: IterationVariable(
         "f_nd_impurity_electrons(03)",
-        data_structure.impurity_radiation_module,
+        "impurity_radiation",
         1e-8,
         0.01,
         target_name="f_nd_impurity_electron_array",
@@ -133,7 +133,7 @@ class IterationVariable:
     ),
     126: IterationVariable(
         "f_nd_impurity_electrons(04)",
-        data_structure.impurity_radiation_module,
+        "impurity_radiation",
         1e-8,
         0.01,
         target_name="f_nd_impurity_electron_array",
@@ -141,7 +141,7 @@ class IterationVariable:
     ),
     127: IterationVariable(
         "f_nd_impurity_electrons(05)",
-        data_structure.impurity_radiation_module,
+        "impurity_radiation",
         1e-8,
         0.01,
         target_name="f_nd_impurity_electron_array",
@@ -149,7 +149,7 @@ class IterationVariable:
     ),
     128: IterationVariable(
         "f_nd_impurity_electrons(06)",
-        data_structure.impurity_radiation_module,
+        "impurity_radiation",
         1e-8,
         0.01,
         target_name="f_nd_impurity_electron_array",
@@ -157,7 +157,7 @@ class IterationVariable:
     ),
     129: IterationVariable(
         "f_nd_impurity_electrons(07)",
-        data_structure.impurity_radiation_module,
+        "impurity_radiation",
         1e-8,
         0.01,
         target_name="f_nd_impurity_electron_array",
@@ -165,7 +165,7 @@ class IterationVariable:
     ),
     130: IterationVariable(
         "f_nd_impurity_electrons(08)",
-        data_structure.impurity_radiation_module,
+        "impurity_radiation",
         1e-8,
         0.01,
         target_name="f_nd_impurity_electron_array",
@@ -173,7 +173,7 @@ class IterationVariable:
     ),
     131: IterationVariable(
         "f_nd_impurity_electrons(09)",
-        data_structure.impurity_radiation_module,
+        "impurity_radiation",
         1e-8,
         0.01,
         target_name="f_nd_impurity_electron_array",
@@ -181,7 +181,7 @@ class IterationVariable:
     ),
     132: IterationVariable(
         "f_nd_impurity_electrons(10)",
-        data_structure.impurity_radiation_module,
+        "impurity_radiation",
         1e-8,
         0.01,
         target_name="f_nd_impurity_electron_array",
@@ -189,7 +189,7 @@ class IterationVariable:
     ),
     133: IterationVariable(
         "f_nd_impurity_electrons(11)",
-        data_structure.impurity_radiation_module,
+        "impurity_radiation",
         1e-8,
         0.01,
         target_name="f_nd_impurity_electron_array",
@@ -197,7 +197,7 @@ class IterationVariable:
     ),
     134: IterationVariable(
         "f_nd_impurity_electrons(12)",
-        data_structure.impurity_radiation_module,
+        "impurity_radiation",
         1e-8,
         0.01,
         target_name="f_nd_impurity_electron_array",
@@ -205,7 +205,7 @@ class IterationVariable:
     ),
     135: IterationVariable(
         "f_nd_impurity_electrons(13)",
-        data_structure.impurity_radiation_module,
+        "impurity_radiation",
         1e-8,
         0.01,
         target_name="f_nd_impurity_electron_array",
@@ -213,7 +213,7 @@ class IterationVariable:
     ),
     136: IterationVariable(
         "f_nd_impurity_electrons(14)",
-        data_structure.impurity_radiation_module,
+        "impurity_radiation",
         1e-8,
         0.01,
         target_name="f_nd_impurity_electron_array",
diff --git a/process/data_structure/__init__.py b/process/data_structure/__init__.py
index 9dc212a3be..97722324db 100644
--- a/process/data_structure/__init__.py
+++ b/process/data_structure/__init__.py
@@ -14,7 +14,7 @@
     global_variables,
     heat_transport_variables,
     ife_variables,
-    impurity_radiation_module,
+    impurity_radiation_variables,
     neoclassics_variables,
     numerics,
     pf_power_variables,
@@ -52,7 +52,7 @@
     "global_variables",
     "heat_transport_variables",
     "ife_variables",
-    "impurity_radiation_module",
+    "impurity_radiation_variables",
     "neoclassics_variables",
     "numerics",
     "pf_power_variables",
diff --git a/process/data_structure/impurity_radiation_module.py b/process/data_structure/impurity_radiation_module.py
deleted file mode 100644
index c5e0c6787f..0000000000
--- a/process/data_structure/impurity_radiation_module.py
+++ /dev/null
@@ -1,113 +0,0 @@
-import numpy as np
-
-N_IMPURITIES = 14
-"""Number of ion species in impurity radiation model
-
-( 1)  Hydrogen  (fraction calculated by code)
-( 2)  Helium
-( 3)  Beryllium
-( 4)  Carbon
-( 5)  Nitrogen
-( 6)  Oxygen
-( 7)  Neon
-( 8)  Silicon
-( 9)  Argon
-(10)  Iron
-(11)  Nickel
-(12)  Krypton
-(13)  Xenon
-(14)  Tungsten
-"""
-
-radius_plasma_core_norm: float = None
-"""Normalised radius defining the 'core' region"""
-
-f_p_plasma_core_rad_reduction: float = None
-"""Fraction of radiation from 'core' region"""
-
-f_nd_impurity_electrons: list[float] = None
-
-n_charge_impurity_profile: list[float] = None
-"""charge profile of impurities"""
-
-imp_label: list[str] = None
-
-impurity_arr_label: list[str] = None
-
-impurity_arr_z: list[float] = None
-
-m_impurity_amu_array: list[float] = None
-"""2D array of impurity atomic masses in Atomic Mass Units (amu)"""
-
-f_nd_impurity_electron_array: list[float] = None
-"""2D array of impurity relative densities (n_imp/n_e)"""
-
-impurity_arr_len_tab: list[int] = None
-
-temp_impurity_keV_array: list[float] = None
-"""2D array of impurity temperatures in kilo-electronvolts (keV)"""
-
-pden_impurity_lz_nd_temp_array: list[float] = None
-
-impurity_arr_zav: list[float] = None
-
-
-def init_impurity_radiation_module():
-    global \
-        radius_plasma_core_norm, \
-        f_p_plasma_core_rad_reduction, \
-        f_nd_impurity_electrons, \
-        n_charge_impurity_profile, \
-        imp_label, \
-        impurity_arr_label, \
-        impurity_arr_z, \
-        m_impurity_amu_array, \
-        f_nd_impurity_electron_array, \
-        impurity_arr_len_tab, \
-        temp_impurity_keV_array, \
-        pden_impurity_lz_nd_temp_array, \
-        impurity_arr_zav
-
-    radius_plasma_core_norm = 0.6
-    f_p_plasma_core_rad_reduction = 1.0
-    f_nd_impurity_electrons = np.array([
-        1.0,
-        0.1,
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-    ])
-    imp_label = np.array([
-        "H_",
-        "He",
-        "Be",
-        "C_",
-        "N_",
-        "O_",
-        "Ne",
-        "Si",
-        "Ar",
-        "Fe",
-        "Ni",
-        "Kr",
-        "Xe",
-        "W_",
-    ])
-    impurity_arr_label = np.full(N_IMPURITIES, "  ")
-    impurity_arr_z = np.zeros(N_IMPURITIES)
-    m_impurity_amu_array = np.zeros(N_IMPURITIES)
-    f_nd_impurity_electron_array = np.zeros(N_IMPURITIES)
-    n_charge_impurity_profile = np.zeros((N_IMPURITIES, 200))
-    impurity_arr_len_tab = np.full(N_IMPURITIES, 0)
-    temp_impurity_keV_array = np.zeros((N_IMPURITIES, 200))
-    pden_impurity_lz_nd_temp_array = np.zeros((N_IMPURITIES, 200))
-    impurity_arr_zav = np.zeros((N_IMPURITIES, 200))
diff --git a/process/data_structure/impurity_radiation_variables.py b/process/data_structure/impurity_radiation_variables.py
new file mode 100644
index 0000000000..ec050687fc
--- /dev/null
+++ b/process/data_structure/impurity_radiation_variables.py
@@ -0,0 +1,110 @@
+from dataclasses import dataclass, field
+
+import numpy as np
+
+N_IMPURITIES = 14
+"""Number of ion species in impurity radiation model
+
+( 1)  Hydrogen  (fraction calculated by code)
+( 2)  Helium
+( 3)  Beryllium
+( 4)  Carbon
+( 5)  Nitrogen
+( 6)  Oxygen
+( 7)  Neon
+( 8)  Silicon
+( 9)  Argon
+(10)  Iron
+(11)  Nickel
+(12)  Krypton
+(13)  Xenon
+(14)  Tungsten
+"""
+
+
+@dataclass(slots=True)
+class ImpurityRadiationData:
+    radius_plasma_core_norm: float = 0.6
+    """Normalised radius defining the 'core' region"""
+
+    f_p_plasma_core_rad_reduction: float = 1.0
+    """Fraction of radiation from 'core' region"""
+
+    f_nd_impurity_electrons: list[float] = field(
+        default_factory=lambda: np.array([
+            1.0,
+            0.1,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+        ])
+    )
+
+    n_charge_impurity_profile: list[float] = field(
+        default_factory=lambda: np.zeros((N_IMPURITIES, 200))
+    )
+    """charge profile of impurities"""
+
+    imp_label: list[str] = field(
+        default_factory=lambda: np.array([
+            "H_",
+            "He",
+            "Be",
+            "C_",
+            "N_",
+            "O_",
+            "Ne",
+            "Si",
+            "Ar",
+            "Fe",
+            "Ni",
+            "Kr",
+            "Xe",
+            "W_",
+        ])
+    )
+
+    impurity_arr_label: list[str] = field(
+        default_factory=lambda: np.full(N_IMPURITIES, "  ")
+    )
+
+    impurity_arr_z: list[float] = field(default_factory=lambda: np.zeros(N_IMPURITIES))
+
+    m_impurity_amu_array: list[float] = field(
+        default_factory=lambda: np.zeros(N_IMPURITIES)
+    )
+    """2D array of impurity atomic masses in Atomic Mass Units (amu)"""
+
+    f_nd_impurity_electron_array: list[float] = field(
+        default_factory=lambda: np.zeros(N_IMPURITIES)
+    )
+    """2D array of impurity relative densities (n_imp/n_e)"""
+
+    impurity_arr_len_tab: list[int] = field(
+        default_factory=lambda: np.full(N_IMPURITIES, 0)
+    )
+
+    temp_impurity_keV_array: list[float] = field(
+        default_factory=lambda: np.zeros((N_IMPURITIES, 200))
+    )
+    """2D array of impurity temperatures in kilo-electronvolts (keV)"""
+
+    pden_impurity_lz_nd_temp_array: list[float] = field(
+        default_factory=lambda: np.zeros((N_IMPURITIES, 200))
+    )
+
+    impurity_arr_zav: list[float] = field(
+        default_factory=lambda: np.zeros((N_IMPURITIES, 200))
+    )
+
+
+CREATE_DICTS_FROM_DATACLASS = ImpurityRadiationData
diff --git a/process/main.py b/process/main.py
index 2257d544e0..22244f0cd3 100644
--- a/process/main.py
+++ b/process/main.py
@@ -306,6 +306,7 @@ def __init__(
         *,
         filepath_out: Path | str | None = None,
         update_obsolete: bool = False,
+        data_structure: DataStructure | None = None,
     ):
         """Read input file and initialise variables.
 
@@ -321,7 +322,7 @@ def __init__(
         self.validate_input(update_obsolete)
         self.init_module_vars()
         self.set_filenames(filepath_out)
-        self.data = DataStructure()
+        self.data = data_structure or DataStructure()
         self.initialise()
         self.models = Models(self.data)
         self.solver = solver
@@ -403,7 +404,7 @@ def initialise(self):
             Path(self.output_path.as_posix().replace("OUT.DAT", "process.log"))
         )
 
-        initialise_imprad()
+        initialise_imprad(self.data)
         # Reads in input file
         init.init_process(self.data)
 
diff --git a/process/models/physics/impurity_radiation.py b/process/models/physics/impurity_radiation.py
index 61209460b8..b5138ff369 100644
--- a/process/models/physics/impurity_radiation.py
+++ b/process/models/physics/impurity_radiation.py
@@ -12,13 +12,13 @@
 
 from process.core import constants
 from process.core.exceptions import ProcessError, ProcessValueError
-from process.data_structure import impurity_radiation_module
+from process.core.model import DataStructure
 from process.models.physics.plasma_profiles import PlasmaProfile
 
 logger = logging.getLogger(__name__)
 
 
-def initialise_imprad():
+def initialise_imprad(data: DataStructure):
     """Initialises the impurity radiation data structure
 
     This routine initialises the impurity radiation data.
@@ -33,12 +33,13 @@ def initialise_imprad():
 
     init_imp_element(
         n_species_index=1,
-        name_label=impurity_radiation_module.imp_label[0],
+        name_label=data.impurity_radiation.imp_label[0],
         z=1,
         m_species_amu=constants.M_PROTIUM_AMU,  # 1.00782503223 1H
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
     f_nd_species_electron = 0.0e0
@@ -46,148 +47,161 @@ def initialise_imprad():
     #  Helium
     init_imp_element(
         n_species_index=2,
-        name_label=impurity_radiation_module.imp_label[1],
+        name_label=data.impurity_radiation.imp_label[1],
         z=2,
         m_species_amu=constants.M_HELIUM_AMU,  # 4.002602 (3He,4He) Average mass
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
     #  Beryllium
     init_imp_element(
         n_species_index=3,
-        name_label=impurity_radiation_module.imp_label[2],
+        name_label=data.impurity_radiation.imp_label[2],
         z=4,
         m_species_amu=constants.M_BERYLLIUM_AMU,  # 9.0121831 9Be
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
     #  Carbon
     init_imp_element(
         n_species_index=4,
-        name_label=impurity_radiation_module.imp_label[3],
+        name_label=data.impurity_radiation.imp_label[3],
         z=6,
         m_species_amu=constants.M_CARBON_AMU,  # 12.0096, (12C,13C,14C) Average mass
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
     #  Nitrogen
     init_imp_element(
         n_species_index=5,
-        name_label=impurity_radiation_module.imp_label[4],
+        name_label=data.impurity_radiation.imp_label[4],
         z=7,
         m_species_amu=constants.M_NITROGEN_AMU,  # 14.00643, (14N,15N) Average mass
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
     #  Oxygen
     init_imp_element(
         n_species_index=6,
-        name_label=impurity_radiation_module.imp_label[5],
+        name_label=data.impurity_radiation.imp_label[5],
         z=8,
         m_species_amu=constants.M_OXYGEN_AMU,  # 15.99903, (16O,17O,18O) Average mass
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
     #  Neon
     init_imp_element(
         n_species_index=7,
-        name_label=impurity_radiation_module.imp_label[6],
+        name_label=data.impurity_radiation.imp_label[6],
         z=10,
         m_species_amu=constants.M_NEON_AMU,  # 20.1797 (20Ne,21Ne,22Ne) Average mass
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
     #  Silicon
     init_imp_element(
         n_species_index=8,
-        name_label=impurity_radiation_module.imp_label[7],
+        name_label=data.impurity_radiation.imp_label[7],
         z=14,
         m_species_amu=constants.M_SILICON_AMU,  # 28.084 (28Si,29Si,30Si) Average mass
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
     #  Argon
     init_imp_element(
         n_species_index=9,
-        name_label=impurity_radiation_module.imp_label[8],
+        name_label=data.impurity_radiation.imp_label[8],
         z=18,
         m_species_amu=constants.M_ARGON_AMU,  # 39.948 (40Ar,36Ar,38Ar) Average mass
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
     #  Iron
     init_imp_element(
         n_species_index=10,
-        name_label=impurity_radiation_module.imp_label[9],
+        name_label=data.impurity_radiation.imp_label[9],
         z=26,
         m_species_amu=constants.M_IRON_AMU,  # 55.845 (56Fe,54Fe,57Fe,58Fe) Average mass
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
     #  Nickel
     init_imp_element(
         n_species_index=11,
-        name_label=impurity_radiation_module.imp_label[10],
+        name_label=data.impurity_radiation.imp_label[10],
         z=28,
         # 58.6934 (58Ni,60Ni,61Ni,62Ni,64Ni) Average mass
         m_species_amu=constants.M_NICKEL_AMU,
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
     #  Krypton
     init_imp_element(
         n_species_index=12,
-        name_label=impurity_radiation_module.imp_label[11],
+        name_label=data.impurity_radiation.imp_label[11],
         z=36,
         # 83.798 (84Kr,86Kr,82Kr,80Kr,78Kr) Average mass
         m_species_amu=constants.M_KRYPTON_AMU,
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
     #  Xenon
     init_imp_element(
         n_species_index=13,
-        name_label=impurity_radiation_module.imp_label[12],
+        name_label=data.impurity_radiation.imp_label[12],
         z=54,
         # 131.293 (132Xe,129Xe,131Xe,134Xe,136Xe) Average mass
         m_species_amu=constants.M_XENON_AMU,
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
     #  Tungsten
     init_imp_element(
         n_species_index=14,
-        name_label=impurity_radiation_module.imp_label[13],
+        name_label=data.impurity_radiation.imp_label[13],
         z=74,
         # 183.84 (184W,186W,182W,183W,180W) Average mass
         m_species_amu=constants.M_TUNGSTEN_AMU,
         f_nd_species_electron=f_nd_species_electron,
         len_tab=table_length,
         error=errorflag,
+        data=data,
     )
 
 
@@ -253,6 +267,7 @@ def init_imp_element(
     f_nd_species_electron: float,
     len_tab: int,
     error: int,
+    data: DataStructure,
 ):
     """Initialise the impurity radiation data for a species.
 
@@ -289,25 +304,25 @@ def init_imp_element(
     if error == 1:
         return
 
-    if n_species_index > len(impurity_radiation_module.impurity_arr_label):
+    if n_species_index > len(data.impurity_radiation.impurity_arr_label):
         raise ProcessValueError(
             "Illegal impurity number",
             number=n_species_index,
-            max=len(impurity_radiation_module.impurity_arr_label),
+            max=len(data.impurity_radiation.impurity_arr_label),
         )
 
-    impurity_radiation_module.impurity_arr_label[n_species_index - 1] = name_label
-    impurity_radiation_module.impurity_arr_z[n_species_index - 1] = z
-    impurity_radiation_module.m_impurity_amu_array[n_species_index - 1] = m_species_amu
-    impurity_radiation_module.f_nd_impurity_electron_array[n_species_index - 1] = (
+    data.impurity_radiation.impurity_arr_label[n_species_index - 1] = name_label
+    data.impurity_radiation.impurity_arr_z[n_species_index - 1] = z
+    data.impurity_radiation.m_impurity_amu_array[n_species_index - 1] = m_species_amu
+    data.impurity_radiation.f_nd_impurity_electron_array[n_species_index - 1] = (
         f_nd_species_electron
     )
-    impurity_radiation_module.impurity_arr_len_tab[n_species_index - 1] = len_tab
+    data.impurity_radiation.impurity_arr_len_tab[n_species_index - 1] = len_tab
 
     if len_tab > 200:
         print(
             f"ERROR: len_tab is {len_tab} but has a maximum value of "
-            f"{impurity_radiation_module.all_array_hotfix_len}"
+            f"{data.impurity_radiation.all_array_hotfix_len}"
         )
 
     impurity_dir = resources.files("process") / "data/lz_non_corona_14_elements/"
@@ -350,37 +365,9 @@ def init_imp_element(
             f"Cannot locate Zav for infinite confinement data in {z_file}"
         )
 
-    impurity_radiation_module.temp_impurity_keV_array[n_species_index - 1, :] = Te * 1e-3
-    impurity_radiation_module.pden_impurity_lz_nd_temp_array[n_species_index - 1, :] = lz
-    impurity_radiation_module.impurity_arr_zav[n_species_index - 1, :] = zav
-
-
-def z2index(zimp):
-    """Finds the index of the impurity element with the given charge number.
-
-    Parameters
-    ----------
-    zimp : int
-        Charge number of the impurity element
-
-    Returns
-    -------
-    int
-        Index of the impurity element in the impurity array
-
-    Raises
-    ------
-    ProcessValueError
-        If the element with the given charge is not found in the impurity array
-    """
-    for i in range(len(impurity_radiation_module.impurity_arr_label)):
-        if zimp == impurity_radiation_module.impurity_arr_z[i]:
-            return i
-
-    # Should only get here if there is a problem
-    raise ProcessValueError(
-        "Element with the given charge is not in the impurity array", zimp=zimp
-    )
+    data.impurity_radiation.temp_impurity_keV_array[n_species_index - 1, :] = Te * 1e-3
+    data.impurity_radiation.pden_impurity_lz_nd_temp_array[n_species_index - 1, :] = lz
+    data.impurity_radiation.impurity_arr_zav[n_species_index - 1, :] = zav
 
 
 def fradcore(rho, radius_plasma_core_norm, f_p_plasma_core_rad_reduction):
@@ -408,7 +395,7 @@ def fradcore(rho, radius_plasma_core_norm, f_p_plasma_core_rad_reduction):
     return fradcore
 
 
-def zav_of_te(imp_element_index, teprofile):
+def zav_of_te(imp_element_index, teprofile, data: DataStructure):
     """Calculates electron temperature dependent average atomic number
 
     Parameters
@@ -428,9 +415,9 @@ def zav_of_te(imp_element_index, teprofile):
     return _zav_of_te_compiled(
         imp_element_index,
         teprofile,
-        impurity_radiation_module.temp_impurity_keV_array,
-        impurity_radiation_module.impurity_arr_zav,
-        impurity_radiation_module.impurity_arr_len_tab,
+        data.impurity_radiation.temp_impurity_keV_array,
+        data.impurity_radiation.impurity_arr_zav,
+        data.impurity_radiation.impurity_arr_len_tab,
     )
 
 
@@ -471,7 +458,7 @@ def _zav_of_te_compiled(
     return zav_of_te
 
 
-def pimpden(imp_element_index, neprofile, teprofile):
+def pimpden(imp_element_index, neprofile, teprofile, data: DataStructure):
     """Calculates the impurity radiation density (W/m3)
 
     Parameters
@@ -490,7 +477,7 @@ def pimpden(imp_element_index, neprofile, teprofile):
     """
     # less_than_imp_temp_mask = teprofile values less than impurity temperature.
     # greater_than_imp_temp_mask = teprofile values higher than impurity temperature.
-    bins = impurity_radiation_module.temp_impurity_keV_array[imp_element_index]
+    bins = data.impurity_radiation.temp_impurity_keV_array[imp_element_index]
     indices = np.digitize(teprofile, bins)
     indices[indices >= bins.shape[0]] = bins.shape[0] - 1
     indices[indices < 0] = 0
@@ -500,10 +487,10 @@ def pimpden(imp_element_index, neprofile, teprofile):
         np.interp(
             np.log(teprofile),
             np.log(
-                impurity_radiation_module.temp_impurity_keV_array[imp_element_index, :]
+                data.impurity_radiation.temp_impurity_keV_array[imp_element_index, :]
             ),
             np.log(
-                impurity_radiation_module.pden_impurity_lz_nd_temp_array[
+                data.impurity_radiation.pden_impurity_lz_nd_temp_array[
                     imp_element_index, :
                 ]
             ),
@@ -511,7 +498,7 @@ def pimpden(imp_element_index, neprofile, teprofile):
     )
 
     pimpden = (
-        impurity_radiation_module.f_nd_impurity_electron_array[imp_element_index]
+        data.impurity_radiation.f_nd_impurity_electron_array[imp_element_index]
         * neprofile
         * neprofile
         * pimpden
@@ -519,31 +506,31 @@ def pimpden(imp_element_index, neprofile, teprofile):
 
     less_than_imp_temp_mask = (
         teprofile
-        <= impurity_radiation_module.temp_impurity_keV_array[imp_element_index, 0]
+        <= data.impurity_radiation.temp_impurity_keV_array[imp_element_index, 0]
     )
     pimpden[less_than_imp_temp_mask] = (
-        impurity_radiation_module.pden_impurity_lz_nd_temp_array[imp_element_index, 0]
+        data.impurity_radiation.pden_impurity_lz_nd_temp_array[imp_element_index, 0]
     )
 
     greater_than_imp_temp_mask = (
         teprofile
-        >= impurity_radiation_module.temp_impurity_keV_array[
+        >= data.impurity_radiation.temp_impurity_keV_array[
             imp_element_index,
-            impurity_radiation_module.impurity_arr_len_tab[imp_element_index] - 1,
+            data.impurity_radiation.impurity_arr_len_tab[imp_element_index] - 1,
         ]
     )
     #  This is okay because Bremsstrahlung will dominate at higher temp.
     pimpden[greater_than_imp_temp_mask] = (
-        impurity_radiation_module.pden_impurity_lz_nd_temp_array[
+        data.impurity_radiation.pden_impurity_lz_nd_temp_array[
             imp_element_index,
-            impurity_radiation_module.impurity_arr_len_tab[imp_element_index] - 1,
+            data.impurity_radiation.impurity_arr_len_tab[imp_element_index] - 1,
         ]
     )
 
     return pimpden
 
 
-def element2index(element: str):
+def element2index(element: str, data: DataStructure):
     """Returns the index of the `element` in the impurity array with
     a given name
 
@@ -559,7 +546,7 @@ def element2index(element: str):
     """
     try:
         return (
-            impurity_radiation_module.impurity_arr_label
+            data.impurity_radiation.impurity_arr_label
             .astype(str)
             .tolist()
             .index(element)
@@ -578,7 +565,7 @@ class ImpurityRadiation:
     radiation loss from each impurity element to find the total impurity radiation loss.
     """
 
-    def __init__(self, plasma_profile: PlasmaProfile):
+    def __init__(self, plasma_profile: PlasmaProfile, data_structure: DataStructure):
         """Initialize the ImpurityRadiation class.
 
         Parameters
@@ -586,11 +573,12 @@ def __init__(self, plasma_profile: PlasmaProfile):
         plasma_profile :
             Parameterises the density and temperature profiles.
         """
+        self.data = data_structure
         self.plasma_profile = plasma_profile
         self.rho = plasma_profile.neprofile.profile_x
         self.rhodx = plasma_profile.neprofile.profile_dx
         self.imp = np.nonzero(
-            impurity_radiation_module.f_nd_impurity_electron_array > 1.0e-30
+            self.data.impurity_radiation.f_nd_impurity_electron_array > 1.0e-30
         )[0]
 
         self.pimp_profile = np.zeros(self.plasma_profile.profile_size)
@@ -625,6 +613,7 @@ def imprad_profile(self, imp_element_index):
             imp_element_index,
             self.plasma_profile.neprofile.profile_y,
             self.plasma_profile.teprofile.profile_y,
+            self.data,
         )
 
         self.pimp_profile = np.add(self.pimp_profile, pimp)
@@ -640,8 +629,8 @@ def calculate_radiation_loss_profiles(self):
             self.rho
             * fradcore(
                 self.rho,
-                impurity_radiation_module.radius_plasma_core_norm,
-                impurity_radiation_module.f_p_plasma_core_rad_reduction,
+                self.data.impurity_radiation.radius_plasma_core_norm,
+                self.data.impurity_radiation.f_p_plasma_core_rad_reduction,
             )
         )
 
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 34ce980ec6..b0df988ca5 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -17,11 +17,8 @@
 from process.core import process_output as po
 from process.core.exceptions import ProcessValueError
 from process.core.model import Model
-from process.data_structure import (
-    impurity_radiation_module,
-    numerics,
-    physics_variables,
-)
+from process.data_structure import numerics, physics_variables
+from process.data_structure.impurity_radiation_variables import N_IMPURITIES
 from process.models.physics import impurity_radiation
 from process.models.physics.plasma_geometry import PlasmaGeom
 from process.models.physics.profiles import PlasmaProfileShapeType
@@ -712,6 +709,7 @@ def run(self):
             physics_variables.rmajor,
             physics_variables.kappa,
             physics_variables.vol_plasma,
+            self.data,
         )
         physics_variables.pden_plasma_sync_mw = radpwrdata.pden_plasma_sync_mw
         physics_variables.pden_plasma_core_rad_mw = radpwrdata.pden_plasma_core_rad_mw
@@ -1061,7 +1059,7 @@ def run(self):
                 self.outfile,
                 (
                     f" 'fzactual, frac, impvardiv = {self.data.reinke.fzactual},"
-                    f" {impurity_radiation_module.f_nd_impurity_electron_array(self.data.reinke.impvardiv)},"
+                    f" {self.data.impurity_radiation.f_nd_impurity_electron_array(self.data.reinke.impvardiv)},"
                     f" {self.data.reinke.impvardiv}"
                 ),
             )
@@ -1166,12 +1164,14 @@ def plasma_composition(self):
 
         # Sum of Zi.ni for all impurity ions (those with charge > helium)
         znimp = 0.0
-        for imp in range(impurity_radiation_module.N_IMPURITIES):
-            if impurity_radiation_module.impurity_arr_z[imp] > 2:
+        for imp in range(N_IMPURITIES):
+            if self.data.impurity_radiation.impurity_arr_z[imp] > 2:
                 znimp += impurity_radiation.zav_of_te(
-                    imp, np.array([physics_variables.temp_plasma_electron_vol_avg_kev])
+                    imp,
+                    np.array([physics_variables.temp_plasma_electron_vol_avg_kev]),
+                    self.data,
                 ).squeeze() * (
-                    impurity_radiation_module.f_nd_impurity_electron_array[imp]
+                    self.data.impurity_radiation.f_nd_impurity_electron_array[imp]
                     * physics_variables.nd_plasma_electrons_vol_avg
                 )
 
@@ -1202,8 +1202,8 @@ def plasma_composition(self):
 
         # Set hydrogen and helium relative impurity densities for
         # radiation calculations
-        impurity_radiation_module.f_nd_impurity_electron_array[
-            impurity_radiation.element2index("H_")
+        self.data.impurity_radiation.f_nd_impurity_electron_array[
+            impurity_radiation.element2index("H_", self.data)
         ] = (
             physics_variables.nd_plasma_protons_vol_avg
             + (
@@ -1214,8 +1214,8 @@ def plasma_composition(self):
             + physics_variables.nd_beam_ions
         ) / physics_variables.nd_plasma_electrons_vol_avg
 
-        impurity_radiation_module.f_nd_impurity_electron_array[
-            impurity_radiation.element2index("He")
+        self.data.impurity_radiation.f_nd_impurity_electron_array[
+            impurity_radiation.element2index("He", self.data)
         ] = (
             physics_variables.f_plasma_fuel_helium3
             * physics_variables.nd_plasma_fuel_ions_vol_avg
@@ -1227,10 +1227,10 @@ def plasma_composition(self):
 
         # Total impurity density
         physics_variables.nd_plasma_impurities_vol_avg = 0.0
-        for imp in range(impurity_radiation_module.N_IMPURITIES):
-            if impurity_radiation_module.impurity_arr_z[imp] > 2:
+        for imp in range(N_IMPURITIES):
+            if self.data.impurity_radiation.impurity_arr_z[imp] > 2:
                 physics_variables.nd_plasma_impurities_vol_avg += (
-                    impurity_radiation_module.f_nd_impurity_electron_array[imp]
+                    self.data.impurity_radiation.f_nd_impurity_electron_array[imp]
                     * physics_variables.nd_plasma_electrons_vol_avg
                 )
 
@@ -1250,21 +1250,21 @@ def plasma_composition(self):
         # Set some relative impurity density variables
         # for the benefit of other routines
         physics_variables.f_nd_plasma_carbon_electron = (
-            impurity_radiation_module.f_nd_impurity_electron_array[
-                impurity_radiation.element2index("C_")
+            self.data.impurity_radiation.f_nd_impurity_electron_array[
+                impurity_radiation.element2index("C_", self.data)
             ]
         )
         physics_variables.f_nd_plasma_oxygen_electron = (
-            impurity_radiation_module.f_nd_impurity_electron_array[
-                impurity_radiation.element2index("O_")
+            self.data.impurity_radiation.f_nd_impurity_electron_array[
+                impurity_radiation.element2index("O_", self.data)
             ]
         )
         physics_variables.f_nd_plasma_iron_argon_electron = (
-            impurity_radiation_module.f_nd_impurity_electron_array[
-                impurity_radiation.element2index("Fe")
+            self.data.impurity_radiation.f_nd_impurity_electron_array[
+                impurity_radiation.element2index("Fe", self.data)
             ]
-            + impurity_radiation_module.f_nd_impurity_electron_array[
-                impurity_radiation.element2index("Ar")
+            + self.data.impurity_radiation.f_nd_impurity_electron_array[
+                impurity_radiation.element2index("Ar", self.data)
             ]
         )
 
@@ -1274,11 +1274,13 @@ def plasma_composition(self):
         # Calculation should be sum(ni.Zi^2) / sum(ni.Zi),
         # but ne = sum(ni.Zi) through quasineutrality
         physics_variables.n_charge_plasma_effective_vol_avg = 0.0
-        for imp in range(impurity_radiation_module.N_IMPURITIES):
+        for imp in range(N_IMPURITIES):
             physics_variables.n_charge_plasma_effective_vol_avg += (
-                impurity_radiation_module.f_nd_impurity_electron_array[imp]
+                self.data.impurity_radiation.f_nd_impurity_electron_array[imp]
                 * impurity_radiation.zav_of_te(
-                    imp, np.array([physics_variables.temp_plasma_electron_vol_avg_kev])
+                    imp,
+                    np.array([physics_variables.temp_plasma_electron_vol_avg_kev]),
+                    self.data,
                 ).squeeze()
                 ** 2
             )
@@ -1337,12 +1339,12 @@ def plasma_composition(self):
             + (physics_variables.nd_plasma_protons_vol_avg * constants.M_PROTON_AMU)
             + (physics_variables.m_beam_amu * physics_variables.nd_beam_ions)
         )
-        for imp in range(impurity_radiation_module.N_IMPURITIES):
-            if impurity_radiation_module.impurity_arr_z[imp] > 2:
+        for imp in range(N_IMPURITIES):
+            if self.data.impurity_radiation.impurity_arr_z[imp] > 2:
                 physics_variables.m_ions_total_amu += (
                     physics_variables.nd_plasma_electrons_vol_avg
-                    * impurity_radiation_module.f_nd_impurity_electron_array[imp]
-                    * impurity_radiation_module.m_impurity_amu_array[imp]
+                    * self.data.impurity_radiation.f_nd_impurity_electron_array[imp]
+                    * self.data.impurity_radiation.m_impurity_amu_array[imp]
                 )
 
         physics_variables.m_ions_total_amu /= (
@@ -1383,16 +1385,17 @@ def plasma_composition(self):
                 / constants.M_TRITON_AMU
             )
         ) / physics_variables.nd_plasma_electrons_vol_avg
-        for imp in range(impurity_radiation_module.N_IMPURITIES):
-            if impurity_radiation_module.impurity_arr_z[imp] > 2:
+        for imp in range(N_IMPURITIES):
+            if self.data.impurity_radiation.impurity_arr_z[imp] > 2:
                 physics_variables.n_charge_plasma_effective_mass_weighted_vol_avg += (
-                    impurity_radiation_module.f_nd_impurity_electron_array[imp]
+                    self.data.impurity_radiation.f_nd_impurity_electron_array[imp]
                     * impurity_radiation.zav_of_te(
                         imp,
                         np.array([physics_variables.temp_plasma_electron_vol_avg_kev]),
+                        self.data,
                     ).squeeze()
                     ** 2
-                    / impurity_radiation_module.m_impurity_amu_array[imp]
+                    / self.data.impurity_radiation.m_impurity_amu_array[imp]
                 )
 
         # ======================================================================
@@ -1649,18 +1652,20 @@ def calculate_effective_charge_ionisation_profiles(self):
         zeff_profile = np.zeros_like(self.plasma_profile.teprofile.profile_y)
         for i in range(len(zeff_profile)):
             zeff_profile[i] = 0.0
-            for imp in range(impurity_radiation_module.N_IMPURITIES):
+            for imp in range(N_IMPURITIES):
                 zeff_profile[i] += (
-                    impurity_radiation_module.f_nd_impurity_electron_array[imp]
+                    self.data.impurity_radiation.f_nd_impurity_electron_array[imp]
                     * impurity_radiation.zav_of_te(
-                        imp, np.array([self.plasma_profile.teprofile.profile_y[i]])
+                        imp,
+                        np.array([self.plasma_profile.teprofile.profile_y[i]]),
+                        self.data,
                     ).squeeze()
                     ** 2
                 )
         physics_variables.n_charge_plasma_effective_profile = zeff_profile
 
         # Assign the charge profiles of each species
-        n_impurities = impurity_radiation_module.N_IMPURITIES
+        n_impurities = N_IMPURITIES
         te_profile = self.plasma_profile.teprofile.profile_y
         n_points = len(te_profile)
         # Create a 2D array: (n_impurities, n_points)
@@ -1668,9 +1673,9 @@ def calculate_effective_charge_ionisation_profiles(self):
         for imp in range(n_impurities):
             for i in range(n_points):
                 charge_profiles[imp, i] = impurity_radiation.zav_of_te(
-                    imp, np.array([te_profile[i]])
+                    imp, np.array([te_profile[i]]), self.data
                 ).squeeze()
-        impurity_radiation_module.n_charge_impurity_profile = charge_profiles
+        self.data.impurity_radiation.n_charge_impurity_profile = charge_profiles
 
     def outplas(self):
         """Subroutine to output the plasma physics information
@@ -2037,13 +2042,13 @@ def outplas(self):
             self.outfile,
             "Plasma normalised minor radius defining 'core' region (ρᵢₙₙₑᵣ)",
             "(radius_plasma_core_norm)",
-            impurity_radiation_module.radius_plasma_core_norm,
+            self.data.impurity_radiation.radius_plasma_core_norm,
         )
         po.ovarre(
             self.outfile,
             "Fractional scaling of radiation power along core profile",
             "(f_p_plasma_core_rad_reduction)",
-            impurity_radiation_module.f_p_plasma_core_rad_reduction,
+            self.data.impurity_radiation.f_p_plasma_core_rad_reduction,
         )
         po.ovarre(
             self.outfile,
@@ -2709,14 +2714,17 @@ def output_temperature_density_profile_info(self) -> None:
         po.ocmmnt(self.outfile, "Plasma ion densities / electron density:")
         po.oblnkl(self.outfile)
 
-        for imp in range(impurity_radiation_module.N_IMPURITIES):
+        for imp in range(N_IMPURITIES):
             # MDK Update f_nd_impurity_electrons, as this will make the ITV output
             # work correctly.
-            impurity_radiation_module.f_nd_impurity_electrons[imp] = (
-                impurity_radiation_module.f_nd_impurity_electron_array[imp]
+            self.data.impurity_radiation.f_nd_impurity_electrons[imp] = (
+                self.data.impurity_radiation.f_nd_impurity_electron_array[imp]
             )
             str1 = (
-                impurity_radiation_module.impurity_arr_label[imp].item().replace("_", "")
+                self.data.impurity_radiation
+                .impurity_arr_label[imp]
+                .item()
+                .replace("_", "")
                 + " concentration"
             )
             str2 = f"(f_nd_impurity_electrons({imp + 1:02}))"
@@ -2726,7 +2734,7 @@ def output_temperature_density_profile_info(self) -> None:
                     self.outfile,
                     str1,
                     str2,
-                    impurity_radiation_module.f_nd_impurity_electrons[imp],
+                    self.data.impurity_radiation.f_nd_impurity_electrons[imp],
                     "OP ",
                 )
             else:
@@ -2734,15 +2742,15 @@ def output_temperature_density_profile_info(self) -> None:
                     self.outfile,
                     str1,
                     str2,
-                    impurity_radiation_module.f_nd_impurity_electrons[imp],
+                    self.data.impurity_radiation.f_nd_impurity_electrons[imp],
                 )
-                if impurity_radiation_module.f_nd_impurity_electrons[imp] != 0.0:  # noqa: RUF069
+                if self.data.impurity_radiation.f_nd_impurity_electrons[imp] != 0.0:  # noqa: RUF069
                     for i in range(physics_variables.n_plasma_profile_elements):
                         po.ovarre(
                             self.mfile,
                             str1 + f" at point {i}",
                             f"(f_nd_impurity_electrons{imp}_{i})",
-                            impurity_radiation_module.f_nd_impurity_electrons[imp]
+                            self.data.impurity_radiation.f_nd_impurity_electrons[imp]
                             * self.plasma_profile.neprofile.profile_y[i],
                             "OP ",
                         )
@@ -2764,13 +2772,13 @@ def output_temperature_density_profile_info(self) -> None:
                 "OP ",
             )
 
-        for imp in range(impurity_radiation_module.N_IMPURITIES):
+        for imp in range(N_IMPURITIES):
             for i in range(physics_variables.n_plasma_profile_elements):
                 po.ovarre(
                     self.mfile,
                     "Impurity charge at point",
                     f"(n_charge_plasma_profile{imp}_{i})",
-                    impurity_radiation_module.n_charge_impurity_profile[imp][i],
+                    self.data.impurity_radiation.n_charge_impurity_profile[imp][i],
                     "OP ",
                 )
 
diff --git a/process/models/physics/radiation_power.py b/process/models/physics/radiation_power.py
index e7d0f20483..0b19e1451b 100644
--- a/process/models/physics/radiation_power.py
+++ b/process/models/physics/radiation_power.py
@@ -10,6 +10,7 @@
 import numpy as np
 
 import process.models.physics.impurity_radiation as impurity
+from process.core.model import DataStructure
 from process.models.physics.plasma_profiles import PlasmaProfile
 
 logger = logging.getLogger(__name__)
@@ -39,6 +40,7 @@ def calculate_radiation_powers(
     rmajor: float,
     kappa: float,
     vol_plasma: float,
+    data_structure: DataStructure,
 ) -> RadpwrData:
     """Calculate the radiation powers in MW/m^3 by calling relevant routines.
 
@@ -98,7 +100,7 @@ def calculate_radiation_powers(
           of arbitrary geometry,” Nuclear Fusion, vol. 41, no. 12, pp. 1755-1758,
           Dec. 2001, doi: https://doi.org/10.1088/0029-5515/41/12/102.
     """
-    imp_rad = impurity.ImpurityRadiation(plasma_profile)
+    imp_rad = impurity.ImpurityRadiation(plasma_profile, data_structure)
     imp_rad.calculate_imprad()
 
     pden_plasma_outer_rad_mw = (
diff --git a/process/models/stellarator/neoclassics.py b/process/models/stellarator/neoclassics.py
index 07fc4f6428..23ecb21c3d 100644
--- a/process/models/stellarator/neoclassics.py
+++ b/process/models/stellarator/neoclassics.py
@@ -4,10 +4,7 @@
 
 from process.core import constants
 from process.core.model import Model
-from process.data_structure import (
-    impurity_radiation_module,
-    physics_variables,
-)
+from process.data_structure import physics_variables
 from process.models.stellarator.stellarator import KEV
 
 logger = logging.getLogger(__name__)
@@ -297,7 +294,7 @@ def calc_neoclassics(self):
             )
             * physics_variables.vol_plasma
             / physics_variables.a_plasma_surface
-            * impurity_radiation_module.radius_plasma_core_norm
+            * self.data.impurity_radiation.radius_plasma_core_norm
         )
         q_PROCESS_r1 = (
             (
@@ -366,7 +363,7 @@ def calc_neoclassics(self):
         dndt_neo_fuel = (
             (dndt_neo_D + dndt_neo_T)
             * physics_variables.a_plasma_surface
-            * impurity_radiation_module.radius_plasma_core_norm
+            * self.data.impurity_radiation.radius_plasma_core_norm
         )
         dmdt_neo_fuel = (
             dndt_neo_fuel * physics_variables.m_fuel_amu * constants.PROTON_MASS * 1.0e6
@@ -375,7 +372,7 @@ def calc_neoclassics(self):
             4
             * dndt_neo_e
             * physics_variables.a_plasma_surface
-            * impurity_radiation_module.radius_plasma_core_norm
+            * self.data.impurity_radiation.radius_plasma_core_norm
             * physics_variables.m_fuel_amu
             * constants.PROTON_MASS
             * 1.0e6
@@ -787,7 +784,7 @@ def st_calc_eff_chi(self):
             physics_variables.vol_plasma
             * self.data.stellarator.f_st_rmajor
             * (
-                impurity_radiation_module.radius_plasma_core_norm
+                self.data.impurity_radiation.radius_plasma_core_norm
                 * physics_variables.rminor
                 / self.data.stellarator_config.stella_config_rminor_ref
             )
@@ -797,7 +794,7 @@ def st_calc_eff_chi(self):
             physics_variables.a_plasma_surface
             * self.data.stellarator.f_st_rmajor
             * (
-                impurity_radiation_module.radius_plasma_core_norm
+                self.data.impurity_radiation.radius_plasma_core_norm
                 * physics_variables.rminor
                 / self.data.stellarator_config.stella_config_rminor_ref
             )
@@ -821,8 +818,8 @@ def st_calc_eff_chi(self):
                 * physics_variables.temp_plasma_electron_on_axis_kev
                 * 1e3
                 * physics_variables.alphat
-                * impurity_radiation_module.radius_plasma_core_norm
-                * (1 - impurity_radiation_module.radius_plasma_core_norm**2)
+                * self.data.impurity_radiation.radius_plasma_core_norm
+                * (1 - self.data.impurity_radiation.radius_plasma_core_norm**2)
                 ** (physics_variables.alphan + physics_variables.alphat - 1)
             )
             * surfacescaling
diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index d2dd447170..92eb5a6b35 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -2163,6 +2163,7 @@ def st_phys(self, output):
             physics_variables.rmajor,
             physics_variables.kappa,
             physics_variables.vol_plasma,
+            self.data,
         )
         physics_variables.pden_plasma_sync_mw = radpwr_data.pden_plasma_sync_mw
         physics_variables.pden_plasma_core_rad_mw = radpwr_data.pden_plasma_core_rad_mw
diff --git a/tests/unit/models/physics/test_impurity_radiation.py b/tests/unit/models/physics/test_impurity_radiation.py
index 18aab625e1..335736faaf 100644
--- a/tests/unit/models/physics/test_impurity_radiation.py
+++ b/tests/unit/models/physics/test_impurity_radiation.py
@@ -5,14 +5,12 @@
 import numpy as np
 import pytest
 
-from process.data_structure import impurity_radiation_module
 from process.models.physics import impurity_radiation
 
 
 @pytest.fixture(autouse=True)
-def initialise_impurity_radiation():
-    impurity_radiation_module.init_impurity_radiation_module()
-    impurity_radiation.initialise_imprad()
+def initialise_impurity_radiation(process_models):
+    impurity_radiation.initialise_imprad(process_models.physics.data)
 
 
 class PimpdenParam(NamedTuple):
@@ -22,7 +20,7 @@ class PimpdenParam(NamedTuple):
     expected_pimpden: np.array = np.array
 
 
-def test_pimpden():
+def test_pimpden(process_models):
     """Tests `pimpden` function.
 
     :param imp_element_index: impurity element
@@ -78,7 +76,10 @@ def test_pimpden():
     )
 
     pimpden = impurity_radiation.pimpden(
-        pimden_parameters.imp_element_index, pimden_parameters.ne, pimden_parameters.te
+        pimden_parameters.imp_element_index,
+        pimden_parameters.ne,
+        pimden_parameters.te,
+        process_models.physics.data,
     )
     assert pytest.approx(pimpden) == pimden_parameters.expected_pimpden
 
@@ -134,7 +135,7 @@ class ZavofteParam(NamedTuple):
     expected_zav_of_te: np.array = np.array
 
 
-def test_zav_of_te():
+def test_zav_of_te(process_models):
     """Tests `Zav_of_te` function.
 
     :param imp_element_index: impurity element
@@ -174,7 +175,7 @@ def test_zav_of_te():
         ]),
     )
     zav_of_te = impurity_radiation.zav_of_te(
-        zavofteparam.imp_element_index, zavofteparam.te
+        zavofteparam.imp_element_index, zavofteparam.te, process_models.physics.data
     )
 
     assert pytest.approx(zav_of_te) == zavofteparam.expected_zav_of_te
diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py
index fd4df1e646..27d4d5e010 100644
--- a/tests/unit/models/physics/test_physics.py
+++ b/tests/unit/models/physics/test_physics.py
@@ -6,10 +6,7 @@
 import pytest
 
 from process.core import constants
-from process.data_structure import (
-    impurity_radiation_module,
-    physics_variables,
-)
+from process.data_structure import physics_variables
 from process.models.physics.impurity_radiation import initialise_imprad
 from process.models.physics.physics import (
     DetailedPhysics,
@@ -1702,7 +1699,7 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics):
     :param monkeypatch: pytest fixture used to mock module/class variables
     :type monkeypatch: _pytest.monkeypatch.monkeypatch
     """
-    initialise_imprad()
+    initialise_imprad(physics.data)
 
     monkeypatch.setattr(
         physics.data.current_drive,
@@ -1711,19 +1708,19 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics):
     )
 
     monkeypatch.setattr(
-        impurity_radiation_module,
+        physics.data.impurity_radiation,
         "f_nd_impurity_electron_array",
         plasmacompositionparam.f_nd_impurity_electron_array,
     )
 
     monkeypatch.setattr(
-        impurity_radiation_module,
+        physics.data.impurity_radiation,
         "impurity_arr_z",
         plasmacompositionparam.impurity_arr_z,
     )
 
     monkeypatch.setattr(
-        impurity_radiation_module,
+        physics.data.impurity_radiation,
         "m_impurity_amu_array",
         plasmacompositionparam.m_impurity_amu_array,
     )
@@ -1916,7 +1913,7 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics):
 
     physics.plasma_composition()
 
-    assert impurity_radiation_module.f_nd_impurity_electron_array == pytest.approx(
+    assert physics.data.impurity_radiation.f_nd_impurity_electron_array == pytest.approx(
         plasmacompositionparam.expected_impurity_arr_frac
     )
 
diff --git a/tests/unit/models/stellarator/test_stellarator.py b/tests/unit/models/stellarator/test_stellarator.py
index 15e45bebac..bd9359394b 100644
--- a/tests/unit/models/stellarator/test_stellarator.py
+++ b/tests/unit/models/stellarator/test_stellarator.py
@@ -3,11 +3,7 @@
 import numpy as np
 import pytest
 
-from process.data_structure import (
-    impurity_radiation_module,
-    physics_variables,
-    tfcoil_variables,
-)
+from process.data_structure import physics_variables, tfcoil_variables
 from process.models.stellarator.build import st_build
 from process.models.stellarator.coils.coils import bmax_from_awp, intersect
 from process.models.stellarator.coils.quench import (
@@ -2110,7 +2106,7 @@ def test_st_calc_eff_chi(stcalceffchiparam, monkeypatch, stellarator):
     monkeypatch.setattr(physics_variables, "rminor", stcalceffchiparam.rminor)
 
     monkeypatch.setattr(
-        impurity_radiation_module,
+        stellarator.data.impurity_radiation,
         "radius_plasma_core_norm",
         stcalceffchiparam.radius_plasma_core_norm,
     )