From f9ca4cfb3d565430f1354bf46d552e33c97ba0f8 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 11:51:03 +0100
Subject: [PATCH 01/23] Add DetailedPhysics class for enhanced plasma
 processing models

---
 process/main.py    |  5 ++++-
 process/physics.py | 10 ++++++++++
 2 files changed, 14 insertions(+), 1 deletion(-)

diff --git a/process/main.py b/process/main.py
index e16b7f2109..7755416e06 100644
--- a/process/main.py
+++ b/process/main.py
@@ -94,7 +94,7 @@
 )
 from process.log import logging_model_handler, show_errors
 from process.pfcoil import PFCoil
-from process.physics import Physics
+from process.physics import DetailedPhysics, Physics
 from process.plasma_geometry import PlasmaGeom
 from process.plasma_profiles import PlasmaProfile
 from process.power import Power
@@ -678,6 +678,9 @@ def __init__(self):
         self.physics = Physics(
             plasma_profile=self.plasma_profile, current_drive=self.current_drive
         )
+        self.physics_detailed = DetailedPhysics(
+            plasma_profile=self.plasma_profile, current_drive=self.current_drive
+        )
         self.neoclassics = Neoclassics()
         self.stellarator = Stellarator(
             availability=self.availability,
diff --git a/process/physics.py b/process/physics.py
index 727859b55b..6b70a11a2e 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9046,3 +9046,13 @@ def reinke_tsep(b_plasma_toroidal_on_axis, flh, qstar, rmajor, eps, fgw, kappa,
         * (eps**0.15 * (1.0 + kappa**2.0) ** 0.34)
         * (lhat**0.29 * kappa_0 ** (-0.29) * 0.285)
     )
+
+
+class DetailedPhysics:
+    """Class to hold detailed physics models for plasma processing."""
+
+    def __init__(self, plasma_profile, current_drive):
+        self.outfile = constants.NOUT
+        self.mfile = constants.MFILE
+        self.plasma_profile = plasma_profile
+        self.current_drive = current_drive

From 55891f68f2281ca7b750354c1ce263ff65975aa8 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 11:55:32 +0100
Subject: [PATCH 02/23] :sparkle: Add Debye length calculations to
 DetailedPhysics class

---
 process/physics.py | 33 +++++++++++++++++++++++++++++++++
 1 file changed, 33 insertions(+)

diff --git a/process/physics.py b/process/physics.py
index 6b70a11a2e..423c89e01a 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9056,3 +9056,36 @@ def __init__(self, plasma_profile, current_drive):
         self.mfile = constants.MFILE
         self.plasma_profile = plasma_profile
         self.current_drive = current_drive
+
+    def calculate_debye_length(
+        self,
+        temp_plasma_species_kev: float,
+        nd_plasma_species: float,
+    ) -> float:
+        """
+        Calculate the Debye length for a plasma.
+
+        :param temp_plasma_species_kev: Species temperature in keV.
+        :type temp_plasma_species_kev: float
+        :param nd_plasma_species: Species number density (/m^3).
+        :type nd_plasma_species: float
+
+        :returns: Debye length in meters.
+        :rtype: float
+        """
+        return (
+            (constants.EPSILON0 * temp_plasma_species_kev)
+            / (nd_plasma_species * constants.ELECTRON_CHARGE**2)
+        ) ** 0.5
+
+    def calculate_debye_length_profile(
+        self,
+        temp_plasma_species_profile_kev: np.ndarray,
+        nd_plasma_species_profile: np.ndarray,
+    ) -> np.ndarray:
+        """
+        Calculate the Debye length profile for a plasma.
+        """
+        return self.calculate_debye_length(
+            temp_plasma_species_profile_kev, nd_plasma_species_profile
+        )

From 7708d9fea37721996dc543acd5ea03d2b9ffb78f Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 11:57:14 +0100
Subject: [PATCH 03/23] :sparkle: Add electron Debye length profile variable to
 physics module

---
 process/data_structure/physics_variables.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index e14f132539..bf54f1c703 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1320,6 +1320,9 @@
 n_charge_plasma_effective_mass_weighted_vol_avg: float = None
 """Plasma mass-weighted volume averaged plasma effective charge"""
 
+len_plasma_debye_electron_profile: list[float] = None
+"""Profile of electron Debye length in plasma (m)"""
+
 
 def init_physics_module():
     """Initialise the physics module"""
@@ -1640,6 +1643,7 @@ def init_physics_variables():
         n_charge_plasma_effective_vol_avg, \
         n_charge_plasma_effective_profile, \
         n_charge_plasma_effective_mass_weighted_vol_avg
+    global len_plasma_debye_electron_profile
 
     m_beam_amu = 0.0
     m_fuel_amu = 0.0
@@ -1900,3 +1904,4 @@ def init_physics_variables():
     n_charge_plasma_effective_vol_avg = 0.0
     n_charge_plasma_effective_profile = []
     n_charge_plasma_effective_mass_weighted_vol_avg = 0.0
+    len_plasma_debye_electron_profile = []

From 237b0dfec4d32fd84f8ae15c6a42e2c1eb713b76 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 13:04:25 +0100
Subject: [PATCH 04/23] Add volume averaged electron Debye length variable to
 physics module

---
 process/data_structure/physics_variables.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index bf54f1c703..322d7a1ca0 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1323,6 +1323,9 @@
 len_plasma_debye_electron_profile: list[float] = None
 """Profile of electron Debye length in plasma (m)"""
 
+len_plasma_debye_electron_vol_avg: float = None
+"""Volume averaged electron Debye length in plasma (m)"""
+
 
 def init_physics_module():
     """Initialise the physics module"""
@@ -1644,6 +1647,7 @@ def init_physics_variables():
         n_charge_plasma_effective_profile, \
         n_charge_plasma_effective_mass_weighted_vol_avg
     global len_plasma_debye_electron_profile
+    global len_plasma_debye_electron_vol_avg
 
     m_beam_amu = 0.0
     m_fuel_amu = 0.0
@@ -1905,3 +1909,4 @@ def init_physics_variables():
     n_charge_plasma_effective_profile = []
     n_charge_plasma_effective_mass_weighted_vol_avg = 0.0
     len_plasma_debye_electron_profile = []
+    len_plasma_debye_electron_vol_avg = 0.0

From 7b6af0babca984a817b8af7be7a1f27a9b155687 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 13:41:18 +0100
Subject: [PATCH 05/23] Add detailed Debye length calculations and output to
 physics module

---
 process/caller.py       |  1 +
 process/io/plot_proc.py | 24 ++++++++++++++++++++++++
 process/output.py       |  1 +
 process/physics.py      | 41 ++++++++++++++++++++++++++++++++++++++++-
 4 files changed, 66 insertions(+), 1 deletion(-)

diff --git a/process/caller.py b/process/caller.py
index 0a93243a15..b73675f71c 100644
--- a/process/caller.py
+++ b/process/caller.py
@@ -253,6 +253,7 @@ def _call_models_once(self, xc: np.ndarray) -> None:
         self.models.build.run()
 
         self.models.physics.physics()
+        self.models.physics_detailed.run()
 
         # Toroidal field coil model
 
diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index 7d175b1741..d1076f1378 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -12269,6 +12269,30 @@ def plot_ebw_ecrh_coupling_graph(axis: plt.Axes, mfile: mf.MFile, scan: int):
     axis.minorticks_on()
 
 
+def plot_debye_length_profile(axis, mfile_data, scan):
+    """Plot the Debye length profile on the given axis."""
+    len_plasma_debye_electron_profile = [
+        mfile_data.data[f"len_plasma_debye_electron_profile{i}"].get_scan(scan)
+        for i in range(500)
+    ]
+
+    axis.plot(
+        np.linspace(0, 1, len(len_plasma_debye_electron_profile)),
+        len_plasma_debye_electron_profile,
+        color="blue",
+        linestyle="-",
+        label=r"$\lambda_{Debye,e}$",
+    )
+
+    axis.set_ylabel("Debye Length [m]")
+    axis.set_xlabel("$\\rho \\ [r/a]$")
+    axis.set_yscale("log")
+    axis.grid(True, which="both", linestyle="--", alpha=0.5)
+    axis.set_xlim([0, 1.025])
+    axis.minorticks_on()
+    axis.legend()
+
+
 def main_plot(
     figs: list[Axes],
     m_file: mf.MFile,
diff --git a/process/output.py b/process/output.py
index 5cd8422f73..765871265b 100644
--- a/process/output.py
+++ b/process/output.py
@@ -43,6 +43,7 @@ def write(models, _outfile):
     # Writing the output from physics.f90 into OUT.DAT + MFILE.DAT
     models.physics.calculate_effective_charge_ionisation_profiles()
     models.physics.outplas()
+    models.physics_detailed.output_detailed_physics()
 
     # TODO what is this? Not in caller.f90?
     models.current_drive.output_current_drive()
diff --git a/process/physics.py b/process/physics.py
index 423c89e01a..98958eeee6 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9057,6 +9057,23 @@ def __init__(self, plasma_profile, current_drive):
         self.plasma_profile = plasma_profile
         self.current_drive = current_drive
 
+    def run(self):
+        # ---------------------------
+        #  Debye length calculation
+        # ---------------------------
+
+        physics_variables.len_plasma_debye_electron_vol_avg = self.calculate_debye_length_profile(
+            temp_plasma_species_profile_kev=physics_variables.temp_plasma_electron_vol_avg_kev,
+            nd_plasma_species_profile=physics_variables.nd_plasma_electrons_vol_avg,
+        )
+
+        physics_variables.len_plasma_debye_electron_profile = (
+            self.calculate_debye_length_profile(
+                temp_plasma_species_profile_kev=self.plasma_profile.teprofile.profile_y,
+                nd_plasma_species_profile=self.plasma_profile.neprofile.profile_y,
+            )
+        )
+
     def calculate_debye_length(
         self,
         temp_plasma_species_kev: float,
@@ -9074,7 +9091,7 @@ def calculate_debye_length(
         :rtype: float
         """
         return (
-            (constants.EPSILON0 * temp_plasma_species_kev)
+            (constants.EPSILON0 * temp_plasma_species_kev * constants.KILOELECTRON_VOLT)
             / (nd_plasma_species * constants.ELECTRON_CHARGE**2)
         ) ** 0.5
 
@@ -9089,3 +9106,25 @@ def calculate_debye_length_profile(
         return self.calculate_debye_length(
             temp_plasma_species_profile_kev, nd_plasma_species_profile
         )
+
+    def output_detailed_physics(self):
+        """Outputs detailed physics variables to file."""
+
+        po.oheadr(self.outfile, "Detailed Plasma")
+
+        po.osubhd(self.outfile, "Debye lengths:")
+
+        po.ovarrf(
+            self.outfile,
+            "Plasma volume averaged electron Debye length (m)",
+            "(len_plasma_debye_electron_vol_avg)",
+            physics_variables.len_plasma_debye_electron_vol_avg,
+            "OP ",
+        )
+        for i in range(len(physics_variables.len_plasma_debye_electron_profile)):
+            po.ovarre(
+                self.mfile,
+                f"Plasma electron Debye length at point {i}",
+                f"len_plasma_debye_electron_profile{i}",
+                physics_variables.len_plasma_debye_electron_profile[i],
+            )

From 58943aaf64f0ac6d851dbd4773c80d8cd92c2f02 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 13:44:41 +0100
Subject: [PATCH 06/23] :sparkle: Add Lorentz factor and relativistic particle
 speed calculations to DetailedPhysics class

---
 process/physics.py | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)

diff --git a/process/physics.py b/process/physics.py
index 98958eeee6..68d5751a2a 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9107,6 +9107,34 @@ def calculate_debye_length_profile(
             temp_plasma_species_profile_kev, nd_plasma_species_profile
         )
 
+    def calculate_lorentz_factor(self, velocity: float) -> float:
+        """
+        Calculate the Lorentz factor for a given velocity.
+        :param velocity: Velocity in m/s.
+        :type velocity: float
+        :returns: Lorentz factor (dimensionless).
+        :rtype: float
+        """
+        return 1 / (1 - (velocity / constants.SPEED_LIGHT) ** 2) ** 0.5
+
+    def calculate_relativistic_particle_speed(
+        self, e_kinetic: float, mass: float
+    ) -> float:
+        """
+        Calculate the speed of a particle given its kinetic energy and mass using relativistic mechanics.
+        :param e_kinetic: Kinetic energy in Joules.
+        :type e_kinetic: float
+        :param mass: Mass of the particle in kg.
+        :type mass: float
+        :returns: Speed of the particle in m/s.
+        :rtype: float
+        """
+        return (
+            constants.SPEED_OF_LIGHT
+            * (1 - (1 / ((e_kinetic / (mass * constants.SPEED_OF_LIGHT**2)) + 1) ** 2))
+            ** 0.5
+        )
+
     def output_detailed_physics(self):
         """Outputs detailed physics variables to file."""
 

From c97f85cf83e7ee70a761f492484053c14c3aa9d6 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 13:46:21 +0100
Subject: [PATCH 07/23] :sparkle: Add electron thermal velocity profile
 variable to physics module

---
 process/data_structure/physics_variables.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 322d7a1ca0..7001ca9258 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1326,6 +1326,9 @@
 len_plasma_debye_electron_vol_avg: float = None
 """Volume averaged electron Debye length in plasma (m)"""
 
+vel_plasma_electron_profile: list[float] = None
+"""Profile of electron thermal velocity in plasma (m/s)"""
+
 
 def init_physics_module():
     """Initialise the physics module"""
@@ -1648,6 +1651,7 @@ def init_physics_variables():
         n_charge_plasma_effective_mass_weighted_vol_avg
     global len_plasma_debye_electron_profile
     global len_plasma_debye_electron_vol_avg
+    global vel_plasma_electron_profile
 
     m_beam_amu = 0.0
     m_fuel_amu = 0.0
@@ -1910,3 +1914,4 @@ def init_physics_variables():
     n_charge_plasma_effective_mass_weighted_vol_avg = 0.0
     len_plasma_debye_electron_profile = []
     len_plasma_debye_electron_vol_avg = 0.0
+    vel_plasma_electron_profile = []

From b8f87a78070abbd47763ddb9502f6bfa2d9666a6 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 13:56:32 +0100
Subject: [PATCH 08/23] :sparkle: Add electron thermal velocity profile
 calculation and plotting functionality

---
 process/constants.py    |  2 +-
 process/io/plot_proc.py | 23 +++++++++++++++++++++++
 process/physics.py      | 26 ++++++++++++++++++++++++--
 3 files changed, 48 insertions(+), 3 deletions(-)

diff --git a/process/constants.py b/process/constants.py
index 98025ec20b..9415549e57 100644
--- a/process/constants.py
+++ b/process/constants.py
@@ -189,7 +189,7 @@
 https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=W
 """
 
-SPEED_LIGHT = 299792458
+SPEED_LIGHT = 299792458.0
 """Speed of light in vacuum (c) [m/s]
 Reference: National Institute of Standards and Technology (NIST)
 https://physics.nist.gov/cgi-bin/cuu/Value?c|search_for=light
diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index d1076f1378..4859ab83a3 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -12293,6 +12293,29 @@ def plot_debye_length_profile(axis, mfile_data, scan):
     axis.legend()
 
 
+def plot_velocity_profile(axis, mfile_data, scan):
+    """Plot the electron thermal velocity profile on the given axis."""
+    vel_plasma_electron_profile = [
+        mfile_data.data[f"vel_plasma_electron_profile{i}"].get_scan(scan)
+        for i in range(500)
+    ]
+
+    axis.plot(
+        np.linspace(0, 1, len(vel_plasma_electron_profile)),
+        vel_plasma_electron_profile,
+        color="blue",
+        linestyle="-",
+        label=r"$v_{e}$",
+    )
+
+    axis.set_ylabel("Velocity [m/s]")
+    axis.set_xlabel("$\\rho \\ [r/a]$")
+    axis.grid(True, which="both", linestyle="--", alpha=0.5)
+    axis.set_xlim([0, 1.025])
+    axis.minorticks_on()
+    axis.legend()
+
+
 def main_plot(
     figs: list[Axes],
     m_file: mf.MFile,
diff --git a/process/physics.py b/process/physics.py
index 68d5751a2a..b5dede895c 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9074,6 +9074,18 @@ def run(self):
             )
         )
 
+        # ============================
+        # Particle relativistic speeds
+        # ============================
+
+        physics_variables.vel_plasma_electron_profile = (
+            self.calculate_relativistic_particle_speed(
+                e_kinetic=self.plasma_profile.teprofile.profile_y
+                * constants.KILOELECTRON_VOLT,
+                mass=constants.ELECTRON_MASS,
+            )
+        )
+
     def calculate_debye_length(
         self,
         temp_plasma_species_kev: float,
@@ -9130,8 +9142,8 @@ def calculate_relativistic_particle_speed(
         :rtype: float
         """
         return (
-            constants.SPEED_OF_LIGHT
-            * (1 - (1 / ((e_kinetic / (mass * constants.SPEED_OF_LIGHT**2)) + 1) ** 2))
+            constants.SPEED_LIGHT
+            * (1 - (1 / ((e_kinetic / (mass * constants.SPEED_LIGHT**2)) + 1) ** 2))
             ** 0.5
         )
 
@@ -9156,3 +9168,13 @@ def output_detailed_physics(self):
                 f"len_plasma_debye_electron_profile{i}",
                 physics_variables.len_plasma_debye_electron_profile[i],
             )
+
+        po.osubhd(self.outfile, "Velocities:")
+
+        for i in range(len(physics_variables.vel_plasma_electron_profile)):
+            po.ovarre(
+                self.mfile,
+                f"Plasma electron thermal velocity at point {i}",
+                f"vel_plasma_electron_profile{i}",
+                physics_variables.vel_plasma_electron_profile[i],
+            )

From 7792d1bf311234bc2dff75ea0c9fc2162b80d836 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 16:08:49 +0100
Subject: [PATCH 09/23] Add Planck's constant and new physics calculations to
 DetailedPhysics class

---
 process/constants.py |  3 +++
 process/physics.py   | 38 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 41 insertions(+)

diff --git a/process/constants.py b/process/constants.py
index 9415549e57..100fe4d657 100644
--- a/process/constants.py
+++ b/process/constants.py
@@ -195,6 +195,9 @@
 https://physics.nist.gov/cgi-bin/cuu/Value?c|search_for=light
 """
 
+PLANCK_CONSTANT = 6.62607015e-34
+"""Planck's constant [J.s]"""
+
 D_T_ENERGY = (
     (DEUTERON_MASS + TRITON_MASS) - (ALPHA_MASS + NEUTRON_MASS)
 ) * SPEED_LIGHT**2
diff --git a/process/physics.py b/process/physics.py
index b5dede895c..484985694c 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9147,6 +9147,44 @@ def calculate_relativistic_particle_speed(
             ** 0.5
         )
 
+    def calculate_coulomb_log_from_impact(
+        self, impact_param_max: float, impact_param_min: float
+    ) -> float:
+        """
+        Calculate the Coulomb logarithm from maximum and minimum impact parameters.
+        :param impact_param_max: Maximum impact parameter in meters.
+        :type impact_param_max: float
+        :param impact_param_min: Minimum impact parameter in meters.
+        :type impact_param_min: float
+        :returns: Coulomb logarithm (dimensionless).
+        :rtype: float
+        """
+        return np.log(impact_param_max / impact_param_min)
+
+    def calculate_classical_distance_of_closest_approach(
+        self,
+        charge1: float,
+        charge2: float,
+        e_kinetic: float,
+    ) -> float:
+        """ """
+
+        return (charge1 * charge2 * constants.ELECTRON_CHARGE**2) / (
+            4 * np.pi * constants.EPSILON0 * e_kinetic
+        )
+
+    def calculate_debroglie_wavelength(self, mass: float, velocity: float) -> float:
+        """
+        Calculate the de Broglie wavelength of a particle.
+        :param mass: Mass of the particle in kg.
+        :type mass: float
+        :param velocity: Velocity of the particle in m/s.
+        :type velocity: float
+        :returns: de Broglie wavelength in meters.
+        :rtype: float
+        """
+        return constants.PLANCK_CONSTANT / (mass * velocity)
+
     def output_detailed_physics(self):
         """Outputs detailed physics variables to file."""
 

From 8666342c8897741f501a5314c8e84d9ab5693032 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 16:17:46 +0100
Subject: [PATCH 10/23] :sparkle: Add electron-electron Coulomb logarithm
 profile variable to physics module

---
 process/data_structure/physics_variables.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 7001ca9258..7e63d0d382 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1329,6 +1329,9 @@
 vel_plasma_electron_profile: list[float] = None
 """Profile of electron thermal velocity in plasma (m/s)"""
 
+plasma_coulomb_log_electron_electron_profile: list[float] = None
+"""Profile of electron-electron Coulomb logarithm in plasma"""
+
 
 def init_physics_module():
     """Initialise the physics module"""
@@ -1652,6 +1655,7 @@ def init_physics_variables():
     global len_plasma_debye_electron_profile
     global len_plasma_debye_electron_vol_avg
     global vel_plasma_electron_profile
+    global plasma_coulomb_log_electron_electron_profile
 
     m_beam_amu = 0.0
     m_fuel_amu = 0.0
@@ -1915,3 +1919,4 @@ def init_physics_variables():
     len_plasma_debye_electron_profile = []
     len_plasma_debye_electron_vol_avg = 0.0
     vel_plasma_electron_profile = []
+    plasma_coulomb_log_electron_electron_profile = []

From 10e1a7ce7992eb096fc1e973157d07c1ffaa2080 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 16:40:12 +0100
Subject: [PATCH 11/23] :sparkle: Add plasma frequency calculation to
 DetailedPhysics class

---
 process/physics.py | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/process/physics.py b/process/physics.py
index 484985694c..b66a46ba88 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9184,6 +9184,28 @@ def calculate_debroglie_wavelength(self, mass: float, velocity: float) -> float:
         :rtype: float
         """
         return constants.PLANCK_CONSTANT / (mass * velocity)
+    
+    def calculate_plasma_frequency(
+        nd_particle: float, m_particle: float, z_particle: float
+    ) -> float:
+        """
+        Calculate the plasma frequency for a particle species.
+        :param nd_particle: Number density of the particle species (/m^3).
+        :type nd_particle: float
+        :param m_particle: Mass of the particle species (kg).
+        :type m_particle: float
+        :param Z_particle: Charge state of the particle species (dimensionless).
+        :type Z_particle: float
+        :returns: Plasma frequency in Hz.
+        :rtype: float
+        """
+        return (
+            (
+                (nd_particle * z_particle**2 * constants.ELECTRON_CHARGE**2)
+                / (m_particle * constants.EPSILON0)
+            )
+            ** 0.5
+        ) / (2 * np.pi)
 
     def output_detailed_physics(self):
         """Outputs detailed physics variables to file."""

From 04a50b1d5f3746fc3455797d7708a8733a480165 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 16:43:04 +0100
Subject: [PATCH 12/23] :sparkle: Add electron plasma frequency profile
 variable to physics module

---
 process/data_structure/physics_variables.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 7e63d0d382..a558296bcd 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1332,6 +1332,9 @@
 plasma_coulomb_log_electron_electron_profile: list[float] = None
 """Profile of electron-electron Coulomb logarithm in plasma"""
 
+freq_plasma_electron_profile: list[float] = None
+"""Electron plasma frequency profile (Hz)"""
+
 
 def init_physics_module():
     """Initialise the physics module"""
@@ -1656,6 +1659,7 @@ def init_physics_variables():
     global len_plasma_debye_electron_vol_avg
     global vel_plasma_electron_profile
     global plasma_coulomb_log_electron_electron_profile
+    global freq_plasma_electron_profile
 
     m_beam_amu = 0.0
     m_fuel_amu = 0.0
@@ -1920,3 +1924,4 @@ def init_physics_variables():
     len_plasma_debye_electron_vol_avg = 0.0
     vel_plasma_electron_profile = []
     plasma_coulomb_log_electron_electron_profile = []
+    freq_plasma_electron_profile = []

From a829a27d385243d388019ecb328264a4424ab54a Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 16:49:27 +0100
Subject: [PATCH 13/23] :sparkle: Add electron thermal frequency profile
 calculation to DetailedPhysics class and plotting functionality

---
 process/io/plot_proc.py | 28 +++++++++++++++++++++++++---
 process/physics.py      | 26 ++++++++++++++++++++++++--
 2 files changed, 49 insertions(+), 5 deletions(-)

diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index 4859ab83a3..5dc68c7b1e 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -12273,7 +12273,7 @@ def plot_debye_length_profile(axis, mfile_data, scan):
     """Plot the Debye length profile on the given axis."""
     len_plasma_debye_electron_profile = [
         mfile_data.data[f"len_plasma_debye_electron_profile{i}"].get_scan(scan)
-        for i in range(500)
+        for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
     ]
 
     axis.plot(
@@ -12286,7 +12286,6 @@ def plot_debye_length_profile(axis, mfile_data, scan):
 
     axis.set_ylabel("Debye Length [m]")
     axis.set_xlabel("$\\rho \\ [r/a]$")
-    axis.set_yscale("log")
     axis.grid(True, which="both", linestyle="--", alpha=0.5)
     axis.set_xlim([0, 1.025])
     axis.minorticks_on()
@@ -12297,7 +12296,7 @@ def plot_velocity_profile(axis, mfile_data, scan):
     """Plot the electron thermal velocity profile on the given axis."""
     vel_plasma_electron_profile = [
         mfile_data.data[f"vel_plasma_electron_profile{i}"].get_scan(scan)
-        for i in range(500)
+        for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
     ]
 
     axis.plot(
@@ -12316,6 +12315,29 @@ def plot_velocity_profile(axis, mfile_data, scan):
     axis.legend()
 
 
+def plot_frequency_profile(axis, mfile_data, scan):
+    """Plot the electron thermal frequency profile on the given axis."""
+    freq_plasma_electron_profile = [
+        mfile_data.data[f"freq_plasma_electron_profile{i}"].get_scan(scan)
+        for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
+    ]
+
+    axis.plot(
+        np.linspace(0, 1, len(freq_plasma_electron_profile)),
+        np.array(freq_plasma_electron_profile) / 1e9,
+        color="blue",
+        linestyle="-",
+        label=r"$\omega_{p,e}$",
+    )
+
+    axis.set_ylabel("Frequency [GHz]")
+    axis.set_xlabel("$\\rho \\ [r/a]$")
+    axis.grid(True, which="both", linestyle="--", alpha=0.5)
+    axis.set_xlim([0, 1.025])
+    axis.minorticks_on()
+    axis.legend()
+
+
 def main_plot(
     figs: list[Axes],
     m_file: mf.MFile,
diff --git a/process/physics.py b/process/physics.py
index b66a46ba88..a594aa3674 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9086,6 +9086,18 @@ def run(self):
             )
         )
 
+        # ============================
+        # Plasma frequencies
+        # ============================
+
+        physics_variables.freq_plasma_electron_profile = (
+            self.calculate_plasma_frequency(
+                nd_particle=self.plasma_profile.neprofile.profile_y,
+                m_particle=constants.ELECTRON_MASS,
+                z_particle=1.0,
+            )
+        )
+
     def calculate_debye_length(
         self,
         temp_plasma_species_kev: float,
@@ -9184,9 +9196,9 @@ def calculate_debroglie_wavelength(self, mass: float, velocity: float) -> float:
         :rtype: float
         """
         return constants.PLANCK_CONSTANT / (mass * velocity)
-    
+
     def calculate_plasma_frequency(
-        nd_particle: float, m_particle: float, z_particle: float
+        self, nd_particle: float, m_particle: float, z_particle: float
     ) -> float:
         """
         Calculate the plasma frequency for a particle species.
@@ -9238,3 +9250,13 @@ def output_detailed_physics(self):
                 f"vel_plasma_electron_profile{i}",
                 physics_variables.vel_plasma_electron_profile[i],
             )
+
+        po.osubhd(self.outfile, "Frequencies:")
+
+        for i in range(len(physics_variables.freq_plasma_electron_profile)):
+            po.ovarre(
+                self.mfile,
+                f"Plasma electron frequency at point {i}",
+                f"freq_plasma_electron_profile{i}",
+                physics_variables.freq_plasma_electron_profile[i],
+            )

From a5a8dbcd7cd89da20bb2f3f095405241bbf6b2da Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 16:55:20 +0100
Subject: [PATCH 14/23] :sparkle: Add Larmor frequency calculation method to
 DetailedPhysics class

---
 process/physics.py | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/process/physics.py b/process/physics.py
index a594aa3674..42a28395c5 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9219,6 +9219,24 @@ def calculate_plasma_frequency(
             ** 0.5
         ) / (2 * np.pi)
 
+    def calculate_larmor_frequency(
+        self, b_field: float, m_particle: float, z_particle: float
+    ) -> float:
+        """
+        Calculate the Larmor frequency for a particle species.
+        :param b_field: Magnetic field strength (T).
+        :type b_field: float
+        :param m_particle: Mass of the particle species (kg).
+        :type m_particle: float
+        :param Z_particle: Charge state of the particle species (dimensionless).
+        :type Z_particle: float
+        :returns: Larmor frequency in Hz.
+        :rtype: float
+        """
+        return (z_particle * constants.ELECTRON_CHARGE * b_field) / (
+            2 * np.pi * m_particle
+        )
+
     def output_detailed_physics(self):
         """Outputs detailed physics variables to file."""
 

From 682eb74740ab31d95bf583939ff78981b8729b31 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 16:56:44 +0100
Subject: [PATCH 15/23] Add electron Larmor frequency profile variable for
 toroidal magnetic field

---
 process/data_structure/physics_variables.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index a558296bcd..68943cd910 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1335,6 +1335,9 @@
 freq_plasma_electron_profile: list[float] = None
 """Electron plasma frequency profile (Hz)"""
 
+freq_plasma_larmor_toroidal_electron_profile: list[float] = None
+"""Profile of electron Larmor frequency in plasma due to toroidal magnetic field (Hz)"""
+
 
 def init_physics_module():
     """Initialise the physics module"""
@@ -1660,6 +1663,7 @@ def init_physics_variables():
     global vel_plasma_electron_profile
     global plasma_coulomb_log_electron_electron_profile
     global freq_plasma_electron_profile
+    global freq_plasma_larmor_toroidal_electron_profile
 
     m_beam_amu = 0.0
     m_fuel_amu = 0.0
@@ -1925,3 +1929,4 @@ def init_physics_variables():
     vel_plasma_electron_profile = []
     plasma_coulomb_log_electron_electron_profile = []
     freq_plasma_electron_profile = []
+    freq_plasma_larmor_toroidal_electron_profile = []

From a4ad341f4ca44a6f0a8ede6fcb904a6c21366a1b Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 18:44:22 +0100
Subject: [PATCH 16/23] :sparkle: Add Larmor frequency calculation for electron
 profile in DetailedPhysics class

---
 process/io/plot_proc.py | 15 +++++++--------
 process/physics.py      | 12 ++++++++++++
 2 files changed, 19 insertions(+), 8 deletions(-)

diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index 5dc68c7b1e..fc12ad5efc 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -12322,18 +12322,17 @@ def plot_frequency_profile(axis, mfile_data, scan):
         for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
     ]
 
-    axis.plot(
-        np.linspace(0, 1, len(freq_plasma_electron_profile)),
-        np.array(freq_plasma_electron_profile) / 1e9,
-        color="blue",
-        linestyle="-",
-        label=r"$\omega_{p,e}$",
-    )
+    x = np.linspace(0, 1, len(freq_plasma_electron_profile))
+    y = np.array(freq_plasma_electron_profile) / 1e9
+    # original curve
+    axis.plot(x, y, color="blue", linestyle="-", label=r"$\omega_{p,e}$")
+    # mirrored across the y-axis (drawn at negative rho)
+    axis.plot(-x, y, color="blue", linestyle="-", label="_nolegend_")
+    axis.set_xlim(-1.025, 1.025)
 
     axis.set_ylabel("Frequency [GHz]")
     axis.set_xlabel("$\\rho \\ [r/a]$")
     axis.grid(True, which="both", linestyle="--", alpha=0.5)
-    axis.set_xlim([0, 1.025])
     axis.minorticks_on()
     axis.legend()
 
diff --git a/process/physics.py b/process/physics.py
index 42a28395c5..4453a661ea 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9097,6 +9097,18 @@ def run(self):
                 z_particle=1.0,
             )
         )
+        
+        # ============================
+        # Larmor frequencies
+        # ============================
+        
+        physics_variables.freq_plasma_electron_larmor_profile = (
+            self.calculate_larmor_frequency(
+                b_field=self.plasma_profile.bprofile.profile_y,
+                m_particle=constants.ELECTRON_MASS,
+                z_particle=1.0,
+            )
+        )
 
     def calculate_debye_length(
         self,

From 5b1a58c05b09d64d9bc8d0da8358d7d7a1001774 Mon Sep 17 00:00:00 2001
From: chris-ashe <christopher.ashe@ukaea.uk>
Date: Wed, 22 Oct 2025 19:17:18 +0100
Subject: [PATCH 17/23] :sparkle: Add Larmor frequency profile for toroidal
 magnetic field in DetailedPhysics class

---
 process/io/plot_proc.py | 15 +++++++++++++++
 process/physics.py      | 17 +++++++++++++----
 2 files changed, 28 insertions(+), 4 deletions(-)

diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index fc12ad5efc..8aa376ef31 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -12321,7 +12321,22 @@ def plot_frequency_profile(axis, mfile_data, scan):
         mfile_data.data[f"freq_plasma_electron_profile{i}"].get_scan(scan)
         for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
     ]
+    freq_plasma_larmor_toroidal_electron_profile = [
+        mfile_data.data[f"freq_plasma_larmor_toroidal_electron_profile{i}"].get_scan(
+            scan
+        )
+        for i in range(
+            2 * int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))
+        )
+    ]
 
+    axis.plot(
+        np.linspace(-1, 1, len(freq_plasma_larmor_toroidal_electron_profile)),
+        np.array(freq_plasma_larmor_toroidal_electron_profile) / 1e9,
+        color="red",
+        linestyle="-",
+        label=r"$f_{Larmor,toroidal,e}$",
+    )
     x = np.linspace(0, 1, len(freq_plasma_electron_profile))
     y = np.array(freq_plasma_electron_profile) / 1e9
     # original curve
diff --git a/process/physics.py b/process/physics.py
index 4453a661ea..da2a3db314 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9097,14 +9097,14 @@ def run(self):
                 z_particle=1.0,
             )
         )
-        
+
         # ============================
         # Larmor frequencies
         # ============================
-        
-        physics_variables.freq_plasma_electron_larmor_profile = (
+
+        physics_variables.freq_plasma_larmor_toroidal_electron_profile = (
             self.calculate_larmor_frequency(
-                b_field=self.plasma_profile.bprofile.profile_y,
+                b_field=physics_variables.b_plasma_toroidal_profile,
                 m_particle=constants.ELECTRON_MASS,
                 z_particle=1.0,
             )
@@ -9290,3 +9290,12 @@ def output_detailed_physics(self):
                 f"freq_plasma_electron_profile{i}",
                 physics_variables.freq_plasma_electron_profile[i],
             )
+        for i in range(
+            len(physics_variables.freq_plasma_larmor_toroidal_electron_profile)
+        ):
+            po.ovarre(
+                self.mfile,
+                f"Plasma electron Larmor frequency at point {i}",
+                f"freq_plasma_larmor_toroidal_electron_profile{i}",
+                physics_variables.freq_plasma_larmor_toroidal_electron_profile[i],
+            )

From 19724191a84dd06be3a7cda816e7e7d1ea7c54ca Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 23 Oct 2025 10:23:45 +0100
Subject: [PATCH 18/23] :sparkle: Add calculation and plotting for plasma
 Coulomb logarithms in DetailedPhysics class

---
 process/io/plot_proc.py | 24 ++++++++++++++++++++++++
 process/physics.py      | 35 +++++++++++++++++++++++++++++++++++
 2 files changed, 59 insertions(+)

diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index 8aa376ef31..013069eb3e 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -12352,6 +12352,30 @@ def plot_frequency_profile(axis, mfile_data, scan):
     axis.legend()
 
 
+def plot_plasma_coloumb_logarithms(axis, mfile_data, scan):
+    """Plot the plasma coloumb logarithms on the given axis."""
+    plasma_coulomb_log_electron_electron_profile = [
+        mfile_data.data[f"plasma_coulomb_log_electron_electron_profile{i}"].get_scan(
+            scan
+        )
+        for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
+    ]
+
+    axis.plot(
+        np.linspace(0, 1, len(plasma_coulomb_log_electron_electron_profile)),
+        plasma_coulomb_log_electron_electron_profile,
+        color="blue",
+        linestyle="-",
+        label=r"$ln \Lambda_{e-e}$",
+    )
+
+    axis.set_ylabel("Coulomb Logarithm")
+    axis.set_xlabel("$\\rho \\ [r/a]$")
+    axis.grid(True, which="both", linestyle="--", alpha=0.5)
+    axis.minorticks_on()
+    axis.legend()
+
+
 def main_plot(
     figs: list[Axes],
     m_file: mf.MFile,
diff --git a/process/physics.py b/process/physics.py
index da2a3db314..3323298427 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9110,6 +9110,29 @@ def run(self):
             )
         )
 
+        # ============================
+        # Coulomb logarithm
+        # ============================
+
+        physics_variables.plasma_coulomb_log_electron_electron_profile = np.array([
+            self.calculate_coulomb_log_from_impact(
+                impact_param_max=physics_variables.len_plasma_debye_electron_profile[i],
+                impact_param_min=max(
+                    self.calculate_classical_distance_of_closest_approach(
+                        charge1=1,
+                        charge2=1,
+                        e_kinetic=self.plasma_profile.teprofile.profile_y[i]
+                        * constants.KILOELECTRON_VOLT,
+                    ),
+                    self.calculate_debroglie_wavelength(
+                        mass=constants.ELECTRON_MASS,
+                        velocity=physics_variables.vel_plasma_electron_profile[i],
+                    ),
+                ),
+            )
+            for i in range(len(physics_variables.len_plasma_debye_electron_profile))
+        ])
+
     def calculate_debye_length(
         self,
         temp_plasma_species_kev: float,
@@ -9299,3 +9322,15 @@ def output_detailed_physics(self):
                 f"freq_plasma_larmor_toroidal_electron_profile{i}",
                 physics_variables.freq_plasma_larmor_toroidal_electron_profile[i],
             )
+
+        po.osubhd(self.outfile, "Coulomb Logarithms:")
+
+        for i in range(
+            len(physics_variables.plasma_coulomb_log_electron_electron_profile)
+        ):
+            po.ovarre(
+                self.mfile,
+                f"Electron-electron Coulomb log at point {i}",
+                f"plasma_coulomb_log_electron_electron_profile{i}",
+                physics_variables.plasma_coulomb_log_electron_electron_profile[i],
+            )

From c43926e9b6b11392daa64eb177c3df64d36f1abc Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 8 Dec 2025 13:28:44 +0000
Subject: [PATCH 19/23] Post rebase changes

---
 process/caller.py       |   1 -
 process/io/plot_proc.py | 304 +++++++++++++++++++++++++++++++++++-----
 process/output.py       |   4 +
 process/physics.py      |   2 +-
 4 files changed, 272 insertions(+), 39 deletions(-)

diff --git a/process/caller.py b/process/caller.py
index b73675f71c..0a93243a15 100644
--- a/process/caller.py
+++ b/process/caller.py
@@ -253,7 +253,6 @@ def _call_models_once(self, xc: np.ndarray) -> None:
         self.models.build.run()
 
         self.models.physics.physics()
-        self.models.physics_detailed.run()
 
         # Toroidal field coil model
 
diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index 013069eb3e..e094c97e7b 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -12276,15 +12276,21 @@ def plot_debye_length_profile(axis, mfile_data, scan):
         for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
     ]
 
+    # Convert to micrometres (1e-6 m)
+    len_plasma_debye_electron_profile_um = [
+        length * 1e6 for length in len_plasma_debye_electron_profile
+    ]
+
     axis.plot(
-        np.linspace(0, 1, len(len_plasma_debye_electron_profile)),
-        len_plasma_debye_electron_profile,
+        np.linspace(0, 1, len(len_plasma_debye_electron_profile_um)),
+        len_plasma_debye_electron_profile_um,
         color="blue",
         linestyle="-",
         label=r"$\lambda_{Debye,e}$",
     )
 
-    axis.set_ylabel("Debye Length [m]")
+    axis.set_ylabel(r"Debye Length [$\mu$m]")
+
     axis.set_xlabel("$\\rho \\ [r/a]$")
     axis.grid(True, which="both", linestyle="--", alpha=0.5)
     axis.set_xlim([0, 1.025])
@@ -12547,9 +12553,14 @@ def main_plot(
     plot_density_limit_comparison(figs[13].add_subplot(221), m_file, scan)
     plot_confinement_time_comparison(figs[13].add_subplot(224), m_file, scan)
 
+    plot_debye_length_profile(figs[14].add_subplot(232), m_file, scan)
+    plot_velocity_profile(figs[14].add_subplot(233), m_file, scan)
+    plot_frequency_profile(figs[14].add_subplot(212), m_file, scan)
+    plot_plasma_coloumb_logarithms(figs[14].add_subplot(231), m_file, scan)
+
     # Plot poloidal cross-section
     poloidal_cross_section(
-        figs[14].add_subplot(121, aspect="equal"),
+        figs[15].add_subplot(121, aspect="equal"),
         m_file,
         scan,
         demo_ranges,
@@ -12559,7 +12570,7 @@ def main_plot(
 
     # Plot toroidal cross-section
     toroidal_cross_section(
-        figs[14].add_subplot(122, aspect="equal"),
+        figs[15].add_subplot(122, aspect="equal"),
         m_file,
         scan,
         demo_ranges,
@@ -12567,19 +12578,19 @@ def main_plot(
     )
 
     # Plot color key
-    ax17 = figs[14].add_subplot(222)
+    ax17 = figs[15].add_subplot(222)
     ax17.set_position([0.5, 0.5, 0.5, 0.5])
     color_key(ax17, m_file, scan, colour_scheme)
 
-    ax18 = figs[15].add_subplot(211)
+    ax18 = figs[16].add_subplot(211)
     ax18.set_position([0.1, 0.33, 0.8, 0.6])
     plot_radial_build(ax18, m_file, colour_scheme)
 
     # Make each axes smaller vertically to leave room for the legend
-    ax185 = figs[16].add_subplot(211)
+    ax185 = figs[17].add_subplot(211)
     ax185.set_position([0.1, 0.61, 0.8, 0.32])
 
-    ax18b = figs[16].add_subplot(212)
+    ax18b = figs[17].add_subplot(212)
     ax18b.set_position([0.1, 0.13, 0.8, 0.32])
     plot_upper_vertical_build(ax185, m_file, colour_scheme)
     plot_lower_vertical_build(ax18b, m_file, colour_scheme)
@@ -12587,50 +12598,54 @@ def main_plot(
     # Can only plot WP and turn structure if superconducting coil at the moment
     if m_file.get("i_tf_sup", scan=scan) == 1:
         # TF coil with WP
-        ax19 = figs[17].add_subplot(221, aspect="equal")
-        # Half height, a bit wider, top left
-        ax19.set_position([0.025, 0.45, 0.5, 0.5])
-        plot_superconducting_tf_wp(ax19, m_file, scan, figs[17])
+        ax19 = figs[18].add_subplot(221, aspect="equal")
+        ax19.set_position([
+            0.025,
+            0.45,
+            0.5,
+            0.5,
+        ])  # Half height, a bit wider, top left
+        plot_superconducting_tf_wp(ax19, m_file, scan, figs[18])
 
         # TF coil turn structure
-        ax20 = figs[18].add_subplot(325, aspect="equal")
+        ax20 = figs[19].add_subplot(325, aspect="equal")
         ax20.set_position([0.025, 0.5, 0.4, 0.4])
-        plot_tf_cable_in_conduit_turn(ax20, figs[18], m_file, scan)
-        plot_205 = figs[18].add_subplot(223, aspect="equal")
+        plot_tf_cable_in_conduit_turn(ax20, figs[19], m_file, scan)
+        plot_205 = figs[19].add_subplot(223, aspect="equal")
         plot_205.set_position([0.075, 0.1, 0.3, 0.3])
-        plot_cable_in_conduit_cable(plot_205, figs[18], m_file, scan)
+        plot_cable_in_conduit_cable(plot_205, figs[19], m_file, scan)
     else:
-        ax19 = figs[17].add_subplot(211, aspect="equal")
+        ax19 = figs[18].add_subplot(211, aspect="equal")
         ax19.set_position([0.06, 0.55, 0.675, 0.4])
-        plot_resistive_tf_wp(ax19, m_file, scan, figs[17])
-
+        plot_resistive_tf_wp(ax19, m_file, scan, figs[18])
     plot_tf_coil_structure(
-        figs[19].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme
+        figs[20].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme
     )
 
-    plot_plasma_outboard_toroidal_ripple_map(figs[20], m_file, scan)
+    plot_plasma_outboard_toroidal_ripple_map(figs[21], m_file, scan)
 
-    plot_tf_stress(figs[21].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file)
+    axes = figs[22].subplots(nrows=3, ncols=1, sharex=True).flatten()
+    plot_tf_stress(axes)
 
-    plot_current_profiles_over_time(figs[22].add_subplot(111), m_file, scan)
+    plot_current_profiles_over_time(figs[23].add_subplot(111), m_file, scan)
 
     plot_cs_coil_structure(
-        figs[23].add_subplot(121, aspect="equal"), figs[23], m_file, scan
+        figs[24].add_subplot(121, aspect="equal"), figs[24], m_file, scan
     )
     plot_cs_turn_structure(
-        figs[23].add_subplot(224, aspect="equal"), figs[23], m_file, scan
+        figs[24].add_subplot(224, aspect="equal"), figs[24], m_file, scan
     )
 
     plot_first_wall_top_down_cross_section(
-        figs[24].add_subplot(221, aspect="equal"), m_file, scan
+        figs[25].add_subplot(221, aspect="equal"), m_file, scan
     )
-    plot_first_wall_poloidal_cross_section(figs[24].add_subplot(122), m_file, scan)
-    plot_fw_90_deg_pipe_bend(figs[24].add_subplot(337), m_file, scan)
+    plot_first_wall_poloidal_cross_section(figs[25].add_subplot(122), m_file, scan)
+    plot_fw_90_deg_pipe_bend(figs[25].add_subplot(337), m_file, scan)
 
-    plot_blkt_pipe_bends(figs[25], m_file, scan)
-    ax_blanket = figs[25].add_subplot(122, aspect="equal")
-    plot_blanket(ax_blanket, m_file, scan, radial_build, colour_scheme)
-    plot_firstwall(ax_blanket, m_file, scan, radial_build, colour_scheme)
+    plot_blkt_pipe_bends(figs[26], m_file, scan)
+    ax_blanket = figs[26].add_subplot(122, aspect="equal")
+    plot_blanket(ax_blanket, m_file_data, scan, radial_build, colour_scheme)
+    plot_firstwall(ax_blanket, m_file_data, scan, radial_build, colour_scheme)
     ax_blanket.set_xlabel("Radial position [m]")
     ax_blanket.set_ylabel("Vertical position [m]")
     ax_blanket.set_title("Blanket and First Wall Poloidal Cross-Section")
@@ -12671,16 +12686,231 @@ def main_plot(
     )
 
     plot_main_power_flow(
-        figs[26].add_subplot(111, aspect="equal"), m_file, scan, figs[26]
+        figs[27].add_subplot(111, aspect="equal"), m_file, scan, figs[27]
     )
 
-    ax24 = figs[27].add_subplot(111)
+    ax24 = figs[28].add_subplot(111)
     # set_position([left, bottom, width, height]) -> height ~ 0.66 => ~2/3 of page height
     ax24.set_position([0.08, 0.35, 0.84, 0.57])
-    plot_system_power_profiles_over_time(ax24, m_file, scan, figs[27])
+    plot_system_power_profiles_over_time(ax24, m_file, scan, figs[28])
+
+def main(args=None):
+    # TODO The use of globals here isn't ideal, but is required to get main()
+    # working with minimal changes. Should be converted to class structure
+    args = parse_args(args)
+    colour_scheme = int(args.colour)
+    # read MFILE
+    m_file = mf.MFile(args.f) if args.f != "" else mf.MFile("MFILE.DAT")
+
+    global m_file_name
+    m_file_name = args.f if args.f != "" else "MFILE.DAT"
+
+    scan = args.n if args.n else -1
+
+    demo_ranges = bool(args.DEMO_ranges)
+
+    # Check for Copper magnets
+    if "i_tf_sup" in m_file.data:
+        i_tf_sup = int(m_file.data["i_tf_sup"].get_scan(scan))
+    else:
+        i_tf_sup = 1
 
+    # Check WP configuration
+    if "i_tf_wp_geom" in m_file.data:
+        i_tf_wp_geom = int(m_file.data["i_tf_wp_geom"].get_scan(scan))
+    else:
+        i_tf_wp_geom = 0
+
+    global dr_bore
+    global dr_cs
+    global dr_cs_tf_gap
+    global dr_tf_inboard
+    global dr_shld_vv_gap_inboard
+    global ddwi
+    global dr_shld_inboard
+    global dr_blkt_inboard
+    global dr_fw_inboard
+    global dr_fw_plasma_gap_inboard
+    global rmajor
+    global rminor
+    global dr_fw_plasma_gap_outboard
+    global dr_fw_outboard
+    global dr_blkt_outboard
+    global dr_shld_outboard
+    global ddwi
+    global dr_shld_vv_gap_outboard
+    global dr_tf_outboard
+    global r_cryostat_inboard
+    global z_cryostat_half_inside
+    global dr_cryostat
+    global j_plasma_0
+
+    dr_bore = m_file.data["dr_bore"].get_scan(scan)
+    dr_cs = m_file.data["dr_cs"].get_scan(scan)
+    dr_cs_tf_gap = m_file.data["dr_cs_tf_gap"].get_scan(scan)
+    dr_tf_inboard = m_file.data["dr_tf_inboard"].get_scan(scan)
+    dr_shld_vv_gap_inboard = m_file.data["dr_shld_vv_gap_inboard"].get_scan(scan)
+    dr_shld_inboard = m_file.data["dr_shld_inboard"].get_scan(scan)
+    dr_blkt_inboard = m_file.data["dr_blkt_inboard"].get_scan(scan)
+    dr_fw_inboard = m_file.data["dr_fw_inboard"].get_scan(scan)
+    dr_fw_plasma_gap_inboard = m_file.data["dr_fw_plasma_gap_inboard"].get_scan(scan)
+    rmajor = m_file.data["rmajor"].get_scan(scan)
+    rminor = m_file.data["rminor"].get_scan(scan)
+    dr_fw_plasma_gap_outboard = m_file.data["dr_fw_plasma_gap_outboard"].get_scan(scan)
+    dr_fw_outboard = m_file.data["dr_fw_outboard"].get_scan(scan)
+    dr_blkt_outboard = m_file.data["dr_blkt_outboard"].get_scan(scan)
+    dr_shld_outboard = m_file.data["dr_shld_outboard"].get_scan(scan)
+    dr_shld_vv_gap_outboard = m_file.data["dr_shld_vv_gap_outboard"].get_scan(scan)
+    dr_tf_outboard = m_file.data["dr_tf_outboard"].get_scan(scan)
+    r_cryostat_inboard = m_file.data["r_cryostat_inboard"].get_scan(scan)
+    z_cryostat_half_inside = m_file.data["z_cryostat_half_inside"].get_scan(scan)
+    dr_cryostat = m_file.data["dr_cryostat"].get_scan(scan)
+    j_plasma_0 = m_file.data["j_plasma_on_axis"].get_scan(scan)
+
+    # Magnets related
+    global n_tf_coils
+    global dx_tf_wp_primary_toroidal
+    global dx_tf_wp_secondary_toroidal
+    global dr_tf_wp_with_insulation
+    global dx_tf_wp_insulation
+    global dr_tf_nose_case
+    global dr_tf_plasma_case
+
+    n_tf_coils = m_file.data["n_tf_coils"].get_scan(scan)
+    if i_tf_sup == 1:  # If superconducting magnets
+        dx_tf_wp_primary_toroidal = m_file.data["dx_tf_wp_primary_toroidal"].get_scan(
+            scan
+        )
+        if i_tf_wp_geom == 1:
+            dx_tf_wp_secondary_toroidal = m_file.data[
+                "dx_tf_wp_secondary_toroidal"
+            ].get_scan(scan)
+        dr_tf_wp_with_insulation = m_file.data["dr_tf_wp_with_insulation"].get_scan(
+            scan
+        )
+        dx_tf_wp_insulation = m_file.data["dx_tf_wp_insulation"].get_scan(scan)
+        dr_tf_nose_case = m_file.data["dr_tf_nose_case"].get_scan(scan)
+
+        # To be re-inergrated to resistives when in-plane stresses is integrated
+        dr_tf_plasma_case = m_file.data["dr_tf_plasma_case"].get_scan(scan)
+
+    global dx_beam_shield
+    global radius_beam_tangency
+    global radius_beam_tangency_max
+    global dx_beam_duct
+
+    i_hcd_primary = int(m_file.data["i_hcd_primary"].get_scan(scan))
+    i_hcd_secondary = int(m_file.data["i_hcd_secondary"].get_scan(scan))
+
+    if (i_hcd_primary in [5, 8]) or (i_hcd_secondary in [5, 8]):
+        dx_beam_shield = m_file.data["dx_beam_shield"].get_scan(scan)
+        radius_beam_tangency = m_file.data["radius_beam_tangency"].get_scan(scan)
+        radius_beam_tangency_max = m_file.data["radius_beam_tangency_max"].get_scan(
+            scan
+        )
+        dx_beam_duct = m_file.data["dx_beam_duct"].get_scan(scan)
+    else:
+        dx_beam_shield = radius_beam_tangency = radius_beam_tangency_max = (
+            dx_beam_duct
+        ) = 0.0
+
+    # Pedestal profile parameters
+    global i_plasma_pedestal
+    global nd_plasma_pedestal_electron
+    global nd_plasma_separatrix_electron
+    global radius_plasma_pedestal_density_norm
+    global radius_plasma_pedestal_temp_norm
+    global tbeta
+    global temp_plasma_pedestal_kev
+    global temp_plasma_separatrix_kev
+    global alphan
+    global alphat
+    global ne0
+    global nd_plasma_fuel_ions_vol_avg
+    global nd_plasma_electrons_vol_avg
+    global te0
+    global ti
+    global te
+    global fgwped_out
+    global fgwsep_out
+    global f_temp_plasma_ion_electron
+
+    i_plasma_pedestal = m_file.data["i_plasma_pedestal"].get_scan(scan)
+    nd_plasma_pedestal_electron = m_file.data["nd_plasma_pedestal_electron"].get_scan(
+        scan
+    )
+    nd_plasma_separatrix_electron = m_file.data[
+        "nd_plasma_separatrix_electron"
+    ].get_scan(scan)
+    radius_plasma_pedestal_density_norm = m_file.data[
+        "radius_plasma_pedestal_density_norm"
+    ].get_scan(scan)
+    radius_plasma_pedestal_temp_norm = m_file.data[
+        "radius_plasma_pedestal_temp_norm"
+    ].get_scan(scan)
+    tbeta = m_file.data["tbeta"].get_scan(scan)
+    temp_plasma_pedestal_kev = m_file.data["temp_plasma_pedestal_kev"].get_scan(scan)
+    temp_plasma_separatrix_kev = m_file.data["temp_plasma_separatrix_kev"].get_scan(
+        scan
+    )
+    alphan = m_file.data["alphan"].get_scan(scan)
+    alphat = m_file.data["alphat"].get_scan(scan)
+    ne0 = m_file.data["nd_plasma_electron_on_axis"].get_scan(scan)
+    nd_plasma_fuel_ions_vol_avg = m_file.data["nd_plasma_fuel_ions_vol_avg"].get_scan(
+        scan
+    )
+    nd_plasma_electrons_vol_avg = m_file.data["nd_plasma_electrons_vol_avg"].get_scan(
+        scan
+    )
+    te0 = m_file.data["temp_plasma_electron_on_axis_kev"].get_scan(scan)
+    ti = m_file.data["temp_plasma_ion_vol_avg_kev"].get_scan(scan)
+    te = m_file.data["temp_plasma_electron_vol_avg_kev"].get_scan(scan)
+    fgwped_out = m_file.data["fgwped_out"].get_scan(scan)
+    fgwsep_out = m_file.data["fgwsep_out"].get_scan(scan)
+    f_temp_plasma_ion_electron = m_file.data["f_temp_plasma_ion_electron"].get_scan(
+        scan
+    )
+
+    # Plasma
+    global triang
+    global alphaj
+    global q0
+    global q95
+    global plasma_current_MA
+    global a_plasma_poloidal
+
+    triang = m_file.data["triang95"].get_scan(scan)
+    alphaj = m_file.data["alphaj"].get_scan(scan)
+    q0 = m_file.data["q0"].get_scan(scan)
+    q95 = m_file.data["q95"].get_scan(scan)
+    plasma_current_MA = m_file.data["plasma_current_ma"].get_scan(scan)
+    a_plasma_poloidal = m_file.data["a_plasma_poloidal"].get_scan(scan)
+
+    # Radial position  -- 0
+    # Electron density -- 1
+    # Electron temperature -- 2
+    # Ion temperature -- 3
+    # Deuterium density -- 4
+    # Tritium density -- 5
+    # BS current density(MA/m^2) -- 6
+    # CD current dens(MA/m^2) -- 7
+    # Total current dens(MA/m^2) -- 8
+    # Poloidal current(R*Bp)(T.m) -- 9
+    # Safety factor q -- 10
+    # Volume (m^3) -- 11
+    # dVolume/dr (m^2) -- 12
+    # Plasma conductivity(MA/(V.m) -- 13
+    # Alpha press(keV*10^10 m^-3) -- 14
+    # Ion dens(10^19 m^-3) -- 15
+    # Poloidal flux (Wb) -- 16
+    # rad profile
+    global f_sync_reflect
+    global b_plasma_toroidal_on_axis
+    global vol_plasma
+    f_sync_reflect = m_file.data["f_sync_reflect"].get_scan(scan)
+    b_plasma_toroidal_on_axis = m_file.data["b_plasma_toroidal_on_axis"].get_scan(scan)
+    vol_plasma = m_file.data["vol_plasma"].get_scan(scan)
 
-def create_thickness_builds(m_file, scan: int):
     # Build the dictionaries of radial and vertical build values and cumulative values
     if int(m_file.get("i_single_null", scan=scan)) == 0:
         vertical_upper = [
diff --git a/process/output.py b/process/output.py
index 765871265b..1fe44b8281 100644
--- a/process/output.py
+++ b/process/output.py
@@ -43,6 +43,10 @@ def write(models, _outfile):
     # Writing the output from physics.f90 into OUT.DAT + MFILE.DAT
     models.physics.calculate_effective_charge_ionisation_profiles()
     models.physics.outplas()
+
+    # Detailed physics, currently only done at final point as values are not used
+    # by any other functions
+    models.physics_detailed.run()
     models.physics_detailed.output_detailed_physics()
 
     # TODO what is this? Not in caller.f90?
diff --git a/process/physics.py b/process/physics.py
index 3323298427..2702f8600c 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9230,7 +9230,7 @@ def calculate_debroglie_wavelength(self, mass: float, velocity: float) -> float:
         :returns: de Broglie wavelength in meters.
         :rtype: float
         """
-        return constants.PLANCK_CONSTANT / (mass * velocity)
+        return (constants.PLANCK_CONSTANT / (2 * np.pi)) / (mass * velocity)
 
     def calculate_plasma_frequency(
         self, nd_particle: float, m_particle: float, z_particle: float

From f49171d8e4da310aff40b287925f747c218020b5 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 8 Dec 2025 15:17:58 +0000
Subject: [PATCH 20/23] :sparkle: Add documentation for Detailed Plasma Physics
 and link in mkdocs

---
 .../physics-models/detailed_physics.md        | 85 +++++++++++++++++++
 mkdocs.yml                                    |  1 +
 2 files changed, 86 insertions(+)
 create mode 100644 documentation/physics-models/detailed_physics.md

diff --git a/documentation/physics-models/detailed_physics.md b/documentation/physics-models/detailed_physics.md
new file mode 100644
index 0000000000..dd01a72ffd
--- /dev/null
+++ b/documentation/physics-models/detailed_physics.md
@@ -0,0 +1,85 @@
+# Detailed Plasma Physics
+
+It can sometimes be useful to calculate rough values for key plasma paramters that are normally used in higher fidelity codes. The `DetailedPhysics()` class stores functions that are called and the end of the run to show rough values for key plasma behavior parameters. The calculation is done at the end as no other methods currently depend on these values.
+
+## Detailed Plasma Physics | `DetailedPhysics()`
+
+
+------------------
+
+### Debye length | `calculate_debye_length()`
+
+Calculates the Debye lenght given by:
+
+$$
+\lambda_{D} = \sqrt{\frac{\epsilon_0 k_B T_e}{n e^2}}
+$$
+
+-------------------
+
+### Relativistic particle speed | `calculate_relativistic_particle_speed()`
+
+$$
+v = c \times \sqrt{\left(1- \frac{1}{\left(1+\frac{E}{mc^2}\right)^2}\right)}
+$$
+
+------------------
+
+### Coulomb Logarithm | `calculate_coulomb_log_from_impact()`
+
+Calculates the Coulomb logarithm assuming a straight line Landau-Spitzer method
+
+$$
+\ln \Lambda = \ln{\left(\frac{b_{\text{max}}}{b_{\text{min}}}\right)}
+$$
+
+The maximum impact parameter is given by the Debye length calculated by [`calculate_debye_length()`](#debye-length--calculate_debye_length)
+$$
+b_{\text{max}} = \lambda_{\text{Debye}}
+$$
+
+The minimum impact paramter is the largest of either the classical distance of closest approach or the Debye length. 
+
+$$
+\begin{split}b_{\text{min}} ≡
+\left\{
+    \begin{array}{ll}
+               λ_{\text{de Broglie}} & \mbox{if } λ_{\text{de Broglie}} ≥ ρ_⟂ \\
+               ρ_⟂         & \mbox{if } ρ_⟂ ≥ λ_{\text{de Broglie}}
+    \end{array}
+\right.\end{split}
+$$
+
+$ρ_⟂$ is the classical distance of closest approach calculated by [`calculate_classical_distance_of_closest_approach()`](#classical-distance-of-closest-approach----calculate_classical_distance_of_closest_approach)
+
+------------------
+
+### Classical distance of closest approach |   `calculate_classical_distance_of_closest_approach()`
+
+$$
+\frac{Z_1Z_2e^2}{4\pi \epsilon_0 E_{\text{kinetic}}}
+$$
+
+---------------------
+
+### DeBroglie Wavelength | `calculate_debroglie_wavelength()`
+
+$$
+\lambda_{\text{DeBroglie}} = \frac{h}{2\pi m v}
+$$
+
+----------------------
+
+### Plasma Frequency | `calculate_plasma_frequency()`
+
+$$
+\omega_p = \sqrt{\frac{n_ie^2}{\epsilon_0 m_i}}
+$$
+
+---------------------
+
+### Larmor Frequency | `calculate_larmor_frequency()`
+
+$$
+f_{\text{Larmor}} = \frac{Z_ieB}{2\pi m_i}
+$$
\ No newline at end of file
diff --git a/mkdocs.yml b/mkdocs.yml
index 22dcfea2dc..ba7cb96aa7 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -70,6 +70,7 @@ nav:
               - Confinement time: physics-models/plasma_confinement.md
               - L-H transition: physics-models/plasma_h_mode.md
               - Plasma Core Power Balance: physics-models/plasma_power_balance.md
+              - Detailed Plasma Physics: physics-models/detailed_physics.md
           - Pulsed Plant Operation: physics-models/pulsed-plant.md
       - Engineering Models:
           - Machine Build: eng-models/machine-build.md

From 08d7b6ea6633f01fa6dc7df5dcaaccabda346f4c Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 5 Jan 2026 14:31:24 +0000
Subject: [PATCH 21/23] Post rebase changes again

---
 process/io/plot_proc.py | 226 +---------------------------------------
 process/physics.py      |  34 ++----
 2 files changed, 16 insertions(+), 244 deletions(-)

diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index e094c97e7b..b64fd5ea91 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -12624,8 +12624,7 @@ def main_plot(
 
     plot_plasma_outboard_toroidal_ripple_map(figs[21], m_file, scan)
 
-    axes = figs[22].subplots(nrows=3, ncols=1, sharex=True).flatten()
-    plot_tf_stress(axes)
+    plot_tf_stress(figs[22].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file)
 
     plot_current_profiles_over_time(figs[23].add_subplot(111), m_file, scan)
 
@@ -12644,8 +12643,8 @@ def main_plot(
 
     plot_blkt_pipe_bends(figs[26], m_file, scan)
     ax_blanket = figs[26].add_subplot(122, aspect="equal")
-    plot_blanket(ax_blanket, m_file_data, scan, radial_build, colour_scheme)
-    plot_firstwall(ax_blanket, m_file_data, scan, radial_build, colour_scheme)
+    plot_blanket(ax_blanket, m_file, scan, radial_build, colour_scheme)
+    plot_firstwall(ax_blanket, m_file, scan, radial_build, colour_scheme)
     ax_blanket.set_xlabel("Radial position [m]")
     ax_blanket.set_ylabel("Vertical position [m]")
     ax_blanket.set_title("Blanket and First Wall Poloidal Cross-Section")
@@ -12694,223 +12693,8 @@ def main_plot(
     ax24.set_position([0.08, 0.35, 0.84, 0.57])
     plot_system_power_profiles_over_time(ax24, m_file, scan, figs[28])
 
-def main(args=None):
-    # TODO The use of globals here isn't ideal, but is required to get main()
-    # working with minimal changes. Should be converted to class structure
-    args = parse_args(args)
-    colour_scheme = int(args.colour)
-    # read MFILE
-    m_file = mf.MFile(args.f) if args.f != "" else mf.MFile("MFILE.DAT")
-
-    global m_file_name
-    m_file_name = args.f if args.f != "" else "MFILE.DAT"
-
-    scan = args.n if args.n else -1
-
-    demo_ranges = bool(args.DEMO_ranges)
-
-    # Check for Copper magnets
-    if "i_tf_sup" in m_file.data:
-        i_tf_sup = int(m_file.data["i_tf_sup"].get_scan(scan))
-    else:
-        i_tf_sup = 1
-
-    # Check WP configuration
-    if "i_tf_wp_geom" in m_file.data:
-        i_tf_wp_geom = int(m_file.data["i_tf_wp_geom"].get_scan(scan))
-    else:
-        i_tf_wp_geom = 0
-
-    global dr_bore
-    global dr_cs
-    global dr_cs_tf_gap
-    global dr_tf_inboard
-    global dr_shld_vv_gap_inboard
-    global ddwi
-    global dr_shld_inboard
-    global dr_blkt_inboard
-    global dr_fw_inboard
-    global dr_fw_plasma_gap_inboard
-    global rmajor
-    global rminor
-    global dr_fw_plasma_gap_outboard
-    global dr_fw_outboard
-    global dr_blkt_outboard
-    global dr_shld_outboard
-    global ddwi
-    global dr_shld_vv_gap_outboard
-    global dr_tf_outboard
-    global r_cryostat_inboard
-    global z_cryostat_half_inside
-    global dr_cryostat
-    global j_plasma_0
-
-    dr_bore = m_file.data["dr_bore"].get_scan(scan)
-    dr_cs = m_file.data["dr_cs"].get_scan(scan)
-    dr_cs_tf_gap = m_file.data["dr_cs_tf_gap"].get_scan(scan)
-    dr_tf_inboard = m_file.data["dr_tf_inboard"].get_scan(scan)
-    dr_shld_vv_gap_inboard = m_file.data["dr_shld_vv_gap_inboard"].get_scan(scan)
-    dr_shld_inboard = m_file.data["dr_shld_inboard"].get_scan(scan)
-    dr_blkt_inboard = m_file.data["dr_blkt_inboard"].get_scan(scan)
-    dr_fw_inboard = m_file.data["dr_fw_inboard"].get_scan(scan)
-    dr_fw_plasma_gap_inboard = m_file.data["dr_fw_plasma_gap_inboard"].get_scan(scan)
-    rmajor = m_file.data["rmajor"].get_scan(scan)
-    rminor = m_file.data["rminor"].get_scan(scan)
-    dr_fw_plasma_gap_outboard = m_file.data["dr_fw_plasma_gap_outboard"].get_scan(scan)
-    dr_fw_outboard = m_file.data["dr_fw_outboard"].get_scan(scan)
-    dr_blkt_outboard = m_file.data["dr_blkt_outboard"].get_scan(scan)
-    dr_shld_outboard = m_file.data["dr_shld_outboard"].get_scan(scan)
-    dr_shld_vv_gap_outboard = m_file.data["dr_shld_vv_gap_outboard"].get_scan(scan)
-    dr_tf_outboard = m_file.data["dr_tf_outboard"].get_scan(scan)
-    r_cryostat_inboard = m_file.data["r_cryostat_inboard"].get_scan(scan)
-    z_cryostat_half_inside = m_file.data["z_cryostat_half_inside"].get_scan(scan)
-    dr_cryostat = m_file.data["dr_cryostat"].get_scan(scan)
-    j_plasma_0 = m_file.data["j_plasma_on_axis"].get_scan(scan)
-
-    # Magnets related
-    global n_tf_coils
-    global dx_tf_wp_primary_toroidal
-    global dx_tf_wp_secondary_toroidal
-    global dr_tf_wp_with_insulation
-    global dx_tf_wp_insulation
-    global dr_tf_nose_case
-    global dr_tf_plasma_case
-
-    n_tf_coils = m_file.data["n_tf_coils"].get_scan(scan)
-    if i_tf_sup == 1:  # If superconducting magnets
-        dx_tf_wp_primary_toroidal = m_file.data["dx_tf_wp_primary_toroidal"].get_scan(
-            scan
-        )
-        if i_tf_wp_geom == 1:
-            dx_tf_wp_secondary_toroidal = m_file.data[
-                "dx_tf_wp_secondary_toroidal"
-            ].get_scan(scan)
-        dr_tf_wp_with_insulation = m_file.data["dr_tf_wp_with_insulation"].get_scan(
-            scan
-        )
-        dx_tf_wp_insulation = m_file.data["dx_tf_wp_insulation"].get_scan(scan)
-        dr_tf_nose_case = m_file.data["dr_tf_nose_case"].get_scan(scan)
-
-        # To be re-inergrated to resistives when in-plane stresses is integrated
-        dr_tf_plasma_case = m_file.data["dr_tf_plasma_case"].get_scan(scan)
-
-    global dx_beam_shield
-    global radius_beam_tangency
-    global radius_beam_tangency_max
-    global dx_beam_duct
-
-    i_hcd_primary = int(m_file.data["i_hcd_primary"].get_scan(scan))
-    i_hcd_secondary = int(m_file.data["i_hcd_secondary"].get_scan(scan))
-
-    if (i_hcd_primary in [5, 8]) or (i_hcd_secondary in [5, 8]):
-        dx_beam_shield = m_file.data["dx_beam_shield"].get_scan(scan)
-        radius_beam_tangency = m_file.data["radius_beam_tangency"].get_scan(scan)
-        radius_beam_tangency_max = m_file.data["radius_beam_tangency_max"].get_scan(
-            scan
-        )
-        dx_beam_duct = m_file.data["dx_beam_duct"].get_scan(scan)
-    else:
-        dx_beam_shield = radius_beam_tangency = radius_beam_tangency_max = (
-            dx_beam_duct
-        ) = 0.0
-
-    # Pedestal profile parameters
-    global i_plasma_pedestal
-    global nd_plasma_pedestal_electron
-    global nd_plasma_separatrix_electron
-    global radius_plasma_pedestal_density_norm
-    global radius_plasma_pedestal_temp_norm
-    global tbeta
-    global temp_plasma_pedestal_kev
-    global temp_plasma_separatrix_kev
-    global alphan
-    global alphat
-    global ne0
-    global nd_plasma_fuel_ions_vol_avg
-    global nd_plasma_electrons_vol_avg
-    global te0
-    global ti
-    global te
-    global fgwped_out
-    global fgwsep_out
-    global f_temp_plasma_ion_electron
-
-    i_plasma_pedestal = m_file.data["i_plasma_pedestal"].get_scan(scan)
-    nd_plasma_pedestal_electron = m_file.data["nd_plasma_pedestal_electron"].get_scan(
-        scan
-    )
-    nd_plasma_separatrix_electron = m_file.data[
-        "nd_plasma_separatrix_electron"
-    ].get_scan(scan)
-    radius_plasma_pedestal_density_norm = m_file.data[
-        "radius_plasma_pedestal_density_norm"
-    ].get_scan(scan)
-    radius_plasma_pedestal_temp_norm = m_file.data[
-        "radius_plasma_pedestal_temp_norm"
-    ].get_scan(scan)
-    tbeta = m_file.data["tbeta"].get_scan(scan)
-    temp_plasma_pedestal_kev = m_file.data["temp_plasma_pedestal_kev"].get_scan(scan)
-    temp_plasma_separatrix_kev = m_file.data["temp_plasma_separatrix_kev"].get_scan(
-        scan
-    )
-    alphan = m_file.data["alphan"].get_scan(scan)
-    alphat = m_file.data["alphat"].get_scan(scan)
-    ne0 = m_file.data["nd_plasma_electron_on_axis"].get_scan(scan)
-    nd_plasma_fuel_ions_vol_avg = m_file.data["nd_plasma_fuel_ions_vol_avg"].get_scan(
-        scan
-    )
-    nd_plasma_electrons_vol_avg = m_file.data["nd_plasma_electrons_vol_avg"].get_scan(
-        scan
-    )
-    te0 = m_file.data["temp_plasma_electron_on_axis_kev"].get_scan(scan)
-    ti = m_file.data["temp_plasma_ion_vol_avg_kev"].get_scan(scan)
-    te = m_file.data["temp_plasma_electron_vol_avg_kev"].get_scan(scan)
-    fgwped_out = m_file.data["fgwped_out"].get_scan(scan)
-    fgwsep_out = m_file.data["fgwsep_out"].get_scan(scan)
-    f_temp_plasma_ion_electron = m_file.data["f_temp_plasma_ion_electron"].get_scan(
-        scan
-    )
-
-    # Plasma
-    global triang
-    global alphaj
-    global q0
-    global q95
-    global plasma_current_MA
-    global a_plasma_poloidal
-
-    triang = m_file.data["triang95"].get_scan(scan)
-    alphaj = m_file.data["alphaj"].get_scan(scan)
-    q0 = m_file.data["q0"].get_scan(scan)
-    q95 = m_file.data["q95"].get_scan(scan)
-    plasma_current_MA = m_file.data["plasma_current_ma"].get_scan(scan)
-    a_plasma_poloidal = m_file.data["a_plasma_poloidal"].get_scan(scan)
-
-    # Radial position  -- 0
-    # Electron density -- 1
-    # Electron temperature -- 2
-    # Ion temperature -- 3
-    # Deuterium density -- 4
-    # Tritium density -- 5
-    # BS current density(MA/m^2) -- 6
-    # CD current dens(MA/m^2) -- 7
-    # Total current dens(MA/m^2) -- 8
-    # Poloidal current(R*Bp)(T.m) -- 9
-    # Safety factor q -- 10
-    # Volume (m^3) -- 11
-    # dVolume/dr (m^2) -- 12
-    # Plasma conductivity(MA/(V.m) -- 13
-    # Alpha press(keV*10^10 m^-3) -- 14
-    # Ion dens(10^19 m^-3) -- 15
-    # Poloidal flux (Wb) -- 16
-    # rad profile
-    global f_sync_reflect
-    global b_plasma_toroidal_on_axis
-    global vol_plasma
-    f_sync_reflect = m_file.data["f_sync_reflect"].get_scan(scan)
-    b_plasma_toroidal_on_axis = m_file.data["b_plasma_toroidal_on_axis"].get_scan(scan)
-    vol_plasma = m_file.data["vol_plasma"].get_scan(scan)
 
+def create_thickness_builds(m_file, scan: int):
     # Build the dictionaries of radial and vertical build values and cumulative values
     if int(m_file.get("i_single_null", scan=scan)) == 0:
         vertical_upper = [
@@ -12984,7 +12768,7 @@ def main(args=None):
 
     # create main plot
     # Increase range when adding new page
-    pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(28)]
+    pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(29)]
 
     # run main_plot
     main_plot(
diff --git a/process/physics.py b/process/physics.py
index 2702f8600c..487048c5e8 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9062,15 +9062,15 @@ def run(self):
         #  Debye length calculation
         # ---------------------------
 
-        physics_variables.len_plasma_debye_electron_vol_avg = self.calculate_debye_length_profile(
-            temp_plasma_species_profile_kev=physics_variables.temp_plasma_electron_vol_avg_kev,
-            nd_plasma_species_profile=physics_variables.nd_plasma_electrons_vol_avg,
+        physics_variables.len_plasma_debye_electron_vol_avg = self.calculate_debye_length(
+            temp_plasma_species_kev=physics_variables.temp_plasma_electron_vol_avg_kev,
+            nd_plasma_species=physics_variables.nd_plasma_electrons_vol_avg,
         )
 
         physics_variables.len_plasma_debye_electron_profile = (
-            self.calculate_debye_length_profile(
-                temp_plasma_species_profile_kev=self.plasma_profile.teprofile.profile_y,
-                nd_plasma_species_profile=self.plasma_profile.neprofile.profile_y,
+            self.calculate_debye_length(
+                temp_plasma_species_kev=self.plasma_profile.teprofile.profile_y,
+                nd_plasma_species=self.plasma_profile.neprofile.profile_y,
             )
         )
 
@@ -9154,18 +9154,6 @@ def calculate_debye_length(
             / (nd_plasma_species * constants.ELECTRON_CHARGE**2)
         ) ** 0.5
 
-    def calculate_debye_length_profile(
-        self,
-        temp_plasma_species_profile_kev: np.ndarray,
-        nd_plasma_species_profile: np.ndarray,
-    ) -> np.ndarray:
-        """
-        Calculate the Debye length profile for a plasma.
-        """
-        return self.calculate_debye_length(
-            temp_plasma_species_profile_kev, nd_plasma_species_profile
-        )
-
     def calculate_lorentz_factor(self, velocity: float) -> float:
         """
         Calculate the Lorentz factor for a given velocity.
@@ -9290,7 +9278,7 @@ def output_detailed_physics(self):
             po.ovarre(
                 self.mfile,
                 f"Plasma electron Debye length at point {i}",
-                f"len_plasma_debye_electron_profile{i}",
+                f"(len_plasma_debye_electron_profile{i})",
                 physics_variables.len_plasma_debye_electron_profile[i],
             )
 
@@ -9300,7 +9288,7 @@ def output_detailed_physics(self):
             po.ovarre(
                 self.mfile,
                 f"Plasma electron thermal velocity at point {i}",
-                f"vel_plasma_electron_profile{i}",
+                f"(vel_plasma_electron_profile{i})",
                 physics_variables.vel_plasma_electron_profile[i],
             )
 
@@ -9310,7 +9298,7 @@ def output_detailed_physics(self):
             po.ovarre(
                 self.mfile,
                 f"Plasma electron frequency at point {i}",
-                f"freq_plasma_electron_profile{i}",
+                f"(freq_plasma_electron_profile{i})",
                 physics_variables.freq_plasma_electron_profile[i],
             )
         for i in range(
@@ -9319,7 +9307,7 @@ def output_detailed_physics(self):
             po.ovarre(
                 self.mfile,
                 f"Plasma electron Larmor frequency at point {i}",
-                f"freq_plasma_larmor_toroidal_electron_profile{i}",
+                f"(freq_plasma_larmor_toroidal_electron_profile{i})",
                 physics_variables.freq_plasma_larmor_toroidal_electron_profile[i],
             )
 
@@ -9331,6 +9319,6 @@ def output_detailed_physics(self):
             po.ovarre(
                 self.mfile,
                 f"Electron-electron Coulomb log at point {i}",
-                f"plasma_coulomb_log_electron_electron_profile{i}",
+                f"(plasma_coulomb_log_electron_electron_profile{i})",
                 physics_variables.plasma_coulomb_log_electron_electron_profile[i],
             )

From 98bb56dae1a4fd89a2e5bfe779692101889026f8 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 7 Jan 2026 11:05:24 +0000
Subject: [PATCH 22/23] Refactor DetailedPhysics methods to static and add unit
 tests for plasma calculations

---
 process/data_structure/physics_variables.py |  14 +-
 process/main.py                             |   2 +-
 process/physics.py                          |  36 +++--
 tests/unit/test_physics.py                  | 166 ++++++++++++++++++++
 4 files changed, 193 insertions(+), 25 deletions(-)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 68943cd910..e7fc6ab59d 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1657,13 +1657,13 @@ def init_physics_variables():
         a_plasma_poloidal, \
         n_charge_plasma_effective_vol_avg, \
         n_charge_plasma_effective_profile, \
-        n_charge_plasma_effective_mass_weighted_vol_avg
-    global len_plasma_debye_electron_profile
-    global len_plasma_debye_electron_vol_avg
-    global vel_plasma_electron_profile
-    global plasma_coulomb_log_electron_electron_profile
-    global freq_plasma_electron_profile
-    global freq_plasma_larmor_toroidal_electron_profile
+        n_charge_plasma_effective_mass_weighted_vol_avg, \
+        len_plasma_debye_electron_profile, \
+        len_plasma_debye_electron_vol_avg, \
+        vel_plasma_electron_profile, \
+        plasma_coulomb_log_electron_electron_profile, \
+        freq_plasma_electron_profile, \
+        freq_plasma_larmor_toroidal_electron_profile
 
     m_beam_amu = 0.0
     m_fuel_amu = 0.0
diff --git a/process/main.py b/process/main.py
index 7755416e06..97c074e6b1 100644
--- a/process/main.py
+++ b/process/main.py
@@ -679,7 +679,7 @@ def __init__(self):
             plasma_profile=self.plasma_profile, current_drive=self.current_drive
         )
         self.physics_detailed = DetailedPhysics(
-            plasma_profile=self.plasma_profile, current_drive=self.current_drive
+            plasma_profile=self.plasma_profile,
         )
         self.neoclassics = Neoclassics()
         self.stellarator = Stellarator(
diff --git a/process/physics.py b/process/physics.py
index 487048c5e8..f6723ae3d6 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9051,11 +9051,10 @@ def reinke_tsep(b_plasma_toroidal_on_axis, flh, qstar, rmajor, eps, fgw, kappa,
 class DetailedPhysics:
     """Class to hold detailed physics models for plasma processing."""
 
-    def __init__(self, plasma_profile, current_drive):
+    def __init__(self, plasma_profile):
         self.outfile = constants.NOUT
         self.mfile = constants.MFILE
         self.plasma_profile = plasma_profile
-        self.current_drive = current_drive
 
     def run(self):
         # ---------------------------
@@ -9090,12 +9089,10 @@ def run(self):
         # Plasma frequencies
         # ============================
 
-        physics_variables.freq_plasma_electron_profile = (
-            self.calculate_plasma_frequency(
-                nd_particle=self.plasma_profile.neprofile.profile_y,
-                m_particle=constants.ELECTRON_MASS,
-                z_particle=1.0,
-            )
+        physics_variables.freq_plasma_electron_profile = self.calculate_plasma_frequency(
+            nd_particle=self.plasma_profile.neprofile.profile_y,
+            m_particle=constants.ELECTRON_MASS,
+            z_particle=1.0,
         )
 
         # ============================
@@ -9133,8 +9130,8 @@ def run(self):
             for i in range(len(physics_variables.len_plasma_debye_electron_profile))
         ])
 
+    @staticmethod
     def calculate_debye_length(
-        self,
         temp_plasma_species_kev: float,
         nd_plasma_species: float,
     ) -> float:
@@ -9154,7 +9151,8 @@ def calculate_debye_length(
             / (nd_plasma_species * constants.ELECTRON_CHARGE**2)
         ) ** 0.5
 
-    def calculate_lorentz_factor(self, velocity: float) -> float:
+    @staticmethod
+    def calculate_lorentz_factor(velocity: float) -> float:
         """
         Calculate the Lorentz factor for a given velocity.
         :param velocity: Velocity in m/s.
@@ -9164,9 +9162,8 @@ def calculate_lorentz_factor(self, velocity: float) -> float:
         """
         return 1 / (1 - (velocity / constants.SPEED_LIGHT) ** 2) ** 0.5
 
-    def calculate_relativistic_particle_speed(
-        self, e_kinetic: float, mass: float
-    ) -> float:
+    @staticmethod
+    def calculate_relativistic_particle_speed(e_kinetic: float, mass: float) -> float:
         """
         Calculate the speed of a particle given its kinetic energy and mass using relativistic mechanics.
         :param e_kinetic: Kinetic energy in Joules.
@@ -9196,8 +9193,8 @@ def calculate_coulomb_log_from_impact(
         """
         return np.log(impact_param_max / impact_param_min)
 
+    @staticmethod
     def calculate_classical_distance_of_closest_approach(
-        self,
         charge1: float,
         charge2: float,
         e_kinetic: float,
@@ -9208,7 +9205,8 @@ def calculate_classical_distance_of_closest_approach(
             4 * np.pi * constants.EPSILON0 * e_kinetic
         )
 
-    def calculate_debroglie_wavelength(self, mass: float, velocity: float) -> float:
+    @staticmethod
+    def calculate_debroglie_wavelength(mass: float, velocity: float) -> float:
         """
         Calculate the de Broglie wavelength of a particle.
         :param mass: Mass of the particle in kg.
@@ -9217,11 +9215,14 @@ def calculate_debroglie_wavelength(self, mass: float, velocity: float) -> float:
         :type velocity: float
         :returns: de Broglie wavelength in meters.
         :rtype: float
+
+        :note: Reduced Planck constant (h-bar) is used in the calculation as this is for scattering.
         """
         return (constants.PLANCK_CONSTANT / (2 * np.pi)) / (mass * velocity)
 
+    @staticmethod
     def calculate_plasma_frequency(
-        self, nd_particle: float, m_particle: float, z_particle: float
+        nd_particle: float, m_particle: float, z_particle: float
     ) -> float:
         """
         Calculate the plasma frequency for a particle species.
@@ -9242,8 +9243,9 @@ def calculate_plasma_frequency(
             ** 0.5
         ) / (2 * np.pi)
 
+    @staticmethod
     def calculate_larmor_frequency(
-        self, b_field: float, m_particle: float, z_particle: float
+        b_field: float, m_particle: float, z_particle: float
     ) -> float:
         """
         Calculate the Larmor frequency for a particle species.
diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py
index a262ba33c1..3caf874982 100644
--- a/tests/unit/test_physics.py
+++ b/tests/unit/test_physics.py
@@ -21,6 +21,7 @@
 )
 from process.impurity_radiation import initialise_imprad
 from process.physics import (
+    DetailedPhysics,
     Physics,
     calculate_beta_limit,
     calculate_current_coefficient_hastie,
@@ -3466,3 +3467,168 @@ def test_calculate_internal_inductance_iter_3():
         b_plasma_poloidal_vol_avg=1.0, c_plasma=1.5e7, vol_plasma=1000.0, rmajor=6.2
     )
     assert result == pytest.approx(0.9078959099585583, abs=0.00001)
+
+
+@pytest.mark.parametrize(
+    "b_field, m_particle, z_particle, expected",
+    [
+        (
+            1.0,
+            constants.ELECTRON_MASS,
+            1.0,
+            2.799249e10,
+        ),  # typical electron in 1T field
+        (2.0, constants.ELECTRON_MASS, 1.0, 5.598498e10),  # double field
+        (1.0, constants.PROTON_MASS, 1.0, 15245186.43761083),  # proton in 1T field
+        (0.5, constants.ELECTRON_MASS, 2.0, 2.799249e10),  # half field, double charge
+        (0.0, constants.ELECTRON_MASS, 1.0, 0.0),  # zero field
+    ],
+)
+def test_calculate_larmor_frequency(b_field, m_particle, z_particle, expected):
+    """Test calculate_larmor_frequency for various particles and fields."""
+    result = DetailedPhysics.calculate_larmor_frequency(
+        b_field=b_field, m_particle=m_particle, z_particle=z_particle
+    )
+    assert result == pytest.approx(expected, rel=1e-5)
+
+
+@pytest.mark.parametrize(
+    "nd_particle,m_particle,z_particle, expected",
+    (
+        (1.0e20, constants.ELECTRON_MASS, 1.0, 89786628157.96086),
+        (1.0e19, constants.PROTON_MASS, 1.0, 662608904.2919972),
+        (5.0e19, constants.PROTON_MASS, 2.0, 2963277104.987116),
+        (0.0, constants.ELECTRON_MASS, 1.0, 0.0),
+    ),
+)
+def test_calculate_plasma_frequency(nd_particle, m_particle, z_particle, expected):
+    """Parametrised tests for DetailedPhysics.calculate_plasma_frequency()."""
+
+    result = DetailedPhysics.calculate_plasma_frequency(
+        nd_particle, m_particle, z_particle
+    )
+    assert result == pytest.approx(expected, rel=1e-12, abs=1e-12)
+
+
+@pytest.mark.parametrize(
+    "mass,velocity,expected",
+    (
+        (
+            constants.ELECTRON_MASS,
+            constants.SPEED_LIGHT * 0.1,
+            3.861592674352376e-12,
+        ),
+        (
+            constants.PROTON_MASS,
+            constants.SPEED_LIGHT * 0.5,
+            4.2061782010279145e-16,
+        ),
+        (
+            constants.PROTON_MASS,
+            1e3,
+            6.304902508360882e-11,
+        ),
+    ),
+)
+def test_calculate_debroglie_wavelength(mass, velocity, expected):
+    """Test DetailedPhysics.calculate_debroglie_wavelength with several parameters."""
+    result = DetailedPhysics.calculate_debroglie_wavelength(mass, velocity)
+    assert result == pytest.approx(expected, rel=1e-12)
+
+
+@pytest.mark.parametrize(
+    "e_kev, mass, expected",
+    (
+        (0.0, constants.ELECTRON_MASS, 0.0),
+        (1.6e-19, constants.ELECTRON_MASS, 592693.0770572403),
+        (1.0e-13, constants.ELECTRON_MASS, 267699064.11978555),
+        (1.0e-10, constants.PROTON_MASS, 239716127.82335472),
+    ),
+)
+def test_calculate_relativistic_particle_speed(e_kev, mass, expected):
+    """Parametrised tests for DetailedPhysics.calculate_relativistic_particle_speed"""
+
+    result = DetailedPhysics.calculate_relativistic_particle_speed(
+        e_kinetic=e_kev, mass=mass
+    )
+    assert result == pytest.approx(expected, rel=1e-5)
+
+
+@pytest.mark.parametrize(
+    "velocity, expected_gamma",
+    (
+        (0.0, 1.0),
+        (0.6 * constants.SPEED_LIGHT, 1.25),
+        (0.99 * constants.SPEED_LIGHT, 7.088812050083354),
+    ),
+)
+def test_calculate_lorentz_factor(velocity, expected_gamma):
+    """Test DetailedPhysics.calculate_lorentz_factor for several velocities."""
+
+    result = DetailedPhysics.calculate_lorentz_factor(velocity=velocity)
+    assert result == pytest.approx(expected_gamma, rel=1e-12)
+
+
+@pytest.mark.parametrize(
+    "temp_keV, nd, expected",
+    (
+        (1.0, 1e19, 7.433941993525029e-05),
+        (10.0, 1e20, 7.433941993525029e-05),
+        (0.1, 1e18, 7.433941993525029e-05),
+    ),
+)
+def test_calculate_debye_length_parametrized(temp_keV, nd, expected):
+    """Parametrized test for DetailedPhysics.calculate_debye_length."""
+    result = DetailedPhysics.calculate_debye_length(temp_keV, nd)
+    assert result == pytest.approx(expected, rel=1e-12)
+
+
+def test_detailed_physics_run_computes_profiles():
+    # Minimal plasma profile
+    plasma = PlasmaProfile()
+    plasma.teprofile.profile_x = np.array([0.0, 0.5, 1.0])
+    plasma.teprofile.profile_y = np.array([1.0, 2.0, 3.0])  # keV
+    plasma.neprofile.profile_x = plasma.teprofile.profile_x
+    plasma.neprofile.profile_y = np.array([1.0e19, 2.0e19, 3.0e19])  # m^-3
+
+    # Set global physics variables required by DetailedPhysics.run
+    physics_variables.temp_plasma_electron_vol_avg_kev = float(
+        np.mean(plasma.teprofile.profile_y)
+    )
+    physics_variables.nd_plasma_electrons_vol_avg = float(
+        np.mean(plasma.neprofile.profile_y)
+    )
+    # toroidal field profile for larmor frequency calc
+    physics_variables.b_plasma_toroidal_profile = (
+        np.ones_like(plasma.neprofile.profile_y) * 5.0
+    )
+
+    dp = DetailedPhysics(plasma)
+
+    # Run should complete without error and populate physics_variables
+    dp.run()
+
+    n = len(plasma.teprofile.profile_y)
+    assert physics_variables.len_plasma_debye_electron_vol_avg > 0
+    assert hasattr(physics_variables, "len_plasma_debye_electron_profile")
+    assert np.shape(physics_variables.len_plasma_debye_electron_profile)[0] == n
+
+    assert np.shape(physics_variables.vel_plasma_electron_profile)[0] == n
+    assert np.all(np.isfinite(physics_variables.vel_plasma_electron_profile))
+
+    assert np.shape(physics_variables.freq_plasma_electron_profile)[0] == n
+    assert np.all(np.isfinite(physics_variables.freq_plasma_electron_profile))
+
+    assert (
+        np.shape(physics_variables.freq_plasma_larmor_toroidal_electron_profile)[0] == n
+    )
+    assert np.all(
+        np.isfinite(physics_variables.freq_plasma_larmor_toroidal_electron_profile)
+    )
+
+    assert (
+        np.shape(physics_variables.plasma_coulomb_log_electron_electron_profile)[0] == n
+    )
+    assert np.all(
+        np.isfinite(physics_variables.plasma_coulomb_log_electron_electron_profile)
+    )

From c04238083d9a23903cbe3a44a6f86de82ed6e3d3 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 7 Jan 2026 14:01:08 +0000
Subject: [PATCH 23/23] Enhance plasma calculation methods to support array
 inputs and update type hints

---
 process/physics.py | 67 ++++++++++++++++++++++++++++------------------
 1 file changed, 41 insertions(+), 26 deletions(-)

diff --git a/process/physics.py b/process/physics.py
index f6723ae3d6..82f3963961 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -9132,19 +9132,19 @@ def run(self):
 
     @staticmethod
     def calculate_debye_length(
-        temp_plasma_species_kev: float,
-        nd_plasma_species: float,
-    ) -> float:
+        temp_plasma_species_kev: float | np.ndarray,
+        nd_plasma_species: float | np.ndarray,
+    ) -> float | np.ndarray:
         """
         Calculate the Debye length for a plasma.
 
         :param temp_plasma_species_kev: Species temperature in keV.
-        :type temp_plasma_species_kev: float
+        :type temp_plasma_species_kev: float | np.ndarray
         :param nd_plasma_species: Species number density (/m^3).
-        :type nd_plasma_species: float
+        :type nd_plasma_species: float | np.ndarray
 
         :returns: Debye length in meters.
-        :rtype: float
+        :rtype: float | np.ndarray
         """
         return (
             (constants.EPSILON0 * temp_plasma_species_kev * constants.KILOELECTRON_VOLT)
@@ -9152,26 +9152,28 @@ def calculate_debye_length(
         ) ** 0.5
 
     @staticmethod
-    def calculate_lorentz_factor(velocity: float) -> float:
+    def calculate_lorentz_factor(velocity: float | np.ndarray) -> float | np.ndarray:
         """
         Calculate the Lorentz factor for a given velocity.
         :param velocity: Velocity in m/s.
-        :type velocity: float
+        :type velocity: float | np.ndarray
         :returns: Lorentz factor (dimensionless).
-        :rtype: float
+        :rtype: float | np.ndarray
         """
         return 1 / (1 - (velocity / constants.SPEED_LIGHT) ** 2) ** 0.5
 
     @staticmethod
-    def calculate_relativistic_particle_speed(e_kinetic: float, mass: float) -> float:
+    def calculate_relativistic_particle_speed(
+        e_kinetic: float | np.ndarray, mass: float
+    ) -> float | np.ndarray:
         """
         Calculate the speed of a particle given its kinetic energy and mass using relativistic mechanics.
         :param e_kinetic: Kinetic energy in Joules.
-        :type e_kinetic: float
+        :type e_kinetic: float | np.ndarray
         :param mass: Mass of the particle in kg.
         :type mass: float
         :returns: Speed of the particle in m/s.
-        :rtype: float
+        :rtype: float | np.ndarray
         """
         return (
             constants.SPEED_LIGHT
@@ -9197,24 +9199,37 @@ def calculate_coulomb_log_from_impact(
     def calculate_classical_distance_of_closest_approach(
         charge1: float,
         charge2: float,
-        e_kinetic: float,
-    ) -> float:
-        """ """
+        e_kinetic: float | np.ndarray,
+    ) -> float | np.ndarray:
+        """
+        Calculate the classical distance of closest approach for two charged particles.
+
+        :param charge1: Charge of particle 1 in units of elementary charge.
+        :type charge1: float
+        :param charge2: Charge of particle 2 in units of elementary charge.
+        :type charge2: float
+        :param e_kinetic: Kinetic energy of the particles in Joules.
+        :type e_kinetic: float | np.ndarray
+        :returns: Distance of closest approach in meters.
+        :rtype: float | np.ndarray
+        """
 
         return (charge1 * charge2 * constants.ELECTRON_CHARGE**2) / (
             4 * np.pi * constants.EPSILON0 * e_kinetic
         )
 
     @staticmethod
-    def calculate_debroglie_wavelength(mass: float, velocity: float) -> float:
+    def calculate_debroglie_wavelength(
+        mass: float, velocity: float | np.ndarray
+    ) -> float | np.ndarray:
         """
         Calculate the de Broglie wavelength of a particle.
         :param mass: Mass of the particle in kg.
         :type mass: float
         :param velocity: Velocity of the particle in m/s.
-        :type velocity: float
+        :type velocity: float | np.ndarray
         :returns: de Broglie wavelength in meters.
-        :rtype: float
+        :rtype: float | np.ndarray
 
         :note: Reduced Planck constant (h-bar) is used in the calculation as this is for scattering.
         """
@@ -9222,18 +9237,18 @@ def calculate_debroglie_wavelength(mass: float, velocity: float) -> float:
 
     @staticmethod
     def calculate_plasma_frequency(
-        nd_particle: float, m_particle: float, z_particle: float
-    ) -> float:
+        nd_particle: float | np.ndarray, m_particle: float, z_particle: float
+    ) -> float | np.ndarray:
         """
         Calculate the plasma frequency for a particle species.
         :param nd_particle: Number density of the particle species (/m^3).
-        :type nd_particle: float
+        :type nd_particle: float | np.ndarray
         :param m_particle: Mass of the particle species (kg).
         :type m_particle: float
         :param Z_particle: Charge state of the particle species (dimensionless).
         :type Z_particle: float
         :returns: Plasma frequency in Hz.
-        :rtype: float
+        :rtype: float | np.ndarray
         """
         return (
             (
@@ -9245,18 +9260,18 @@ def calculate_plasma_frequency(
 
     @staticmethod
     def calculate_larmor_frequency(
-        b_field: float, m_particle: float, z_particle: float
-    ) -> float:
+        b_field: float | np.ndarray, m_particle: float, z_particle: float
+    ) -> float | np.ndarray:
         """
         Calculate the Larmor frequency for a particle species.
         :param b_field: Magnetic field strength (T).
-        :type b_field: float
+        :type b_field: float | np.ndarray
         :param m_particle: Mass of the particle species (kg).
         :type m_particle: float
         :param Z_particle: Charge state of the particle species (dimensionless).
         :type Z_particle: float
         :returns: Larmor frequency in Hz.
-        :rtype: float
+        :rtype: float | np.ndarray
         """
         return (z_particle * constants.ELECTRON_CHARGE * b_field) / (
             2 * np.pi * m_particle