From 0a8d1c6b5fecbd9af7167e765c3177c04a142030 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 1 Apr 2026 15:13:52 +0100
Subject: [PATCH 01/26] Rename hoop_stress method to calculate_cs_hoop_stress
 and update references

---
 documentation/source/eng-models/central-solenoid.md | 2 +-
 process/models/pfcoil.py                            | 4 ++--
 tests/unit/models/test_pfcoil.py                    | 4 ++--
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/documentation/source/eng-models/central-solenoid.md b/documentation/source/eng-models/central-solenoid.md
index 354aa3ae10..c1cb06dc50 100644
--- a/documentation/source/eng-models/central-solenoid.md
+++ b/documentation/source/eng-models/central-solenoid.md
@@ -157,7 +157,7 @@ to the maximum at the midplane, $F_{z}(0)$.
 --------------------------
 
 
-### Hoop stress | `hoop_stress()`
+### Hoop stress | `calculate_cs_hoop_stress()`
 
 
 
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 71f73afe9c..c7bd931b14 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -3340,7 +3340,7 @@ def ohcalc(self):
             # Superconducting coil
 
             # New calculation from M. N. Wilson for hoop stress
-            self.data.pf_coil.sig_hoop = self.hoop_stress(
+            self.data.pf_coil.sig_hoop = self.calculate_cs_hoop_stress(
                 self.data.pf_coil.r_pf_coil_inner[self.data.pf_coil.n_cs_pf_coils - 1]
             )
 
@@ -3813,7 +3813,7 @@ def calculate_cs_self_midplane_axial_stress_time_profile(
             )
             self.data.pf_coil.stress_z_cs_self_midplane_profile[time] = stress_value
 
-    def hoop_stress(self, r):
+    def calculate_cs_hoop_stress(self, r):
         """Calculation of hoop stress of central solenoid.
 
         This routine calculates the hoop stress of the central solenoid
diff --git a/tests/unit/models/test_pfcoil.py b/tests/unit/models/test_pfcoil.py
index 5ca0ba00ed..47dd35bbbd 100644
--- a/tests/unit/models/test_pfcoil.py
+++ b/tests/unit/models/test_pfcoil.py
@@ -1970,7 +1970,7 @@ def test_vsec(pfcoil, monkeypatch):
 def test_hoop_stress(cs_coil, monkeypatch):
     """Test hoop_stress subroutine.
 
-    hoop_stress() requires specific mocked variables in order to work; these were
+    calculate_cs_hoop_stress() requires specific mocked variables in order to work; these were
     discovered using gdb to break on the first subroutine call when running the
     baseline 2018 IN.DAT.
 
@@ -2045,7 +2045,7 @@ def test_hoop_stress(cs_coil, monkeypatch):
 
     r = 2.3
     s_hoop_exp = 6.737108e8
-    s_hoop = cs_coil.hoop_stress(r)
+    s_hoop = cs_coil.calculate_cs_hoop_stress(r)
 
     assert pytest.approx(s_hoop) == s_hoop_exp
 

From 36ef7baa7192c80ee068a02e35bf2241b91c94e8 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 08:57:43 +0100
Subject: [PATCH 02/26] Enhance CS hoop stress calculation with additional
 parameters and update documentation

---
 .../source/eng-models/central-solenoid.md     | 12 +++++++
 process/models/pfcoil.py                      | 33 ++++++++++++-------
 2 files changed, 33 insertions(+), 12 deletions(-)

diff --git a/documentation/source/eng-models/central-solenoid.md b/documentation/source/eng-models/central-solenoid.md
index c1cb06dc50..4a62e1e6c4 100644
--- a/documentation/source/eng-models/central-solenoid.md
+++ b/documentation/source/eng-models/central-solenoid.md
@@ -166,6 +166,18 @@ The hoop stress is calculated using equations 4.10 and 4.11 from "Superconductin
 Wilson (1983).  This is divided by the fraction of the area occupied by steel to obtain the hoop 
 stress in the steel, $\sigma_{hoop}$.
 
+$$
+\sigma_{\theta} = \frac{K(2+v)}{3(\alpha+1)}\times \left(\alpha^2+\alpha+1+\frac{\alpha^2}{\eps^2}-\eps \frac{(1+2v)(\alpha+1)}{(2+v)}\right) \\
+- \frac{M(3+v)}{8}\left(\alpha^2+1+\frac{\alpha^2}{\eps^2}-\frac{(1+3v)}{(3+v)}\eps^2\right)
+$$
+
+Where:
+
+- $\epsilon = \frac{r}{r_{CS,inner}}$
+- $\alpha = \frac{r_{CS,outer}}{r_{CS,inner}}$
+
+For a special infinite solenoid case $B$
+
 The axial stress can be calculated using "Case studies in superconducting magnets", Y. Iwasa, p. 
 86, 3.5.2, Special Case 4: Midplane force.  This applies exactly only to a thin-walled solenoid. 
 The axial stress in the steel is given by:
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index c7bd931b14..9b784a4992 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -3813,7 +3813,15 @@ def calculate_cs_self_midplane_axial_stress_time_profile(
             )
             self.data.pf_coil.stress_z_cs_self_midplane_profile[time] = stress_value
 
-    def calculate_cs_hoop_stress(self, r):
+    def calculate_cs_hoop_stress(
+        self,
+        r_stress_point,
+        r_cs_inner=None,
+        r_cs_outer=None,
+        dz_cs_half=None,
+        j_cs=None,
+        b_cs_inner=None,
+    ) -> float:
         """Calculation of hoop stress of central solenoid.
 
         This routine calculates the hoop stress of the central solenoid
@@ -3822,39 +3830,40 @@ def calculate_cs_hoop_stress(self, r):
         Parameters
         ----------
         r : float
-            radial position a < r < b
+            radial position r_cs_inner < r < b
 
+        
+        
+        
+        
         Returns
         -------
         float
             hoop stress (MPa)
         """
-        a = self.data.pf_coil.r_pf_coil_inner[self.data.pf_coil.n_cs_pf_coils - 1]
 
-        # Outer radius of central Solenoid [m]
-        b = self.data.pf_coil.r_pf_coil_outer[self.data.pf_coil.n_cs_pf_coils - 1]
+        beta = dz_cs_half / r_cs_inner
+        alpha = r_cs_outer / r_cs_inner
 
         # alpha
-        alpha = b / a
+        alpha = r_cs_outer / r_cs_inner
 
         # epsilon
-        epsilon = r / a
+        epsilon = r_stress_point / r_cs_inner
 
-        # Field at inner radius of coil [T]
-        b_a = self.data.pf_coil.b_cs_peak_pulse_start
 
         # Field at outer radius of coil [T]
         # Assume to be 0 for now
         b_b = 0.0e0
 
         # current density [A/m^2]
-        j = self.data.pf_coil.j_cs_pulse_start
+        j_cs = pfcoil_variables.j_cs_pulse_start
 
         # K term
-        k = ((alpha * b_a - b_b) * j * a) / (alpha - 1.0e0)
+        k = ((alpha * b_cs_inner - b_b) * j_cs * r_cs_inner) / (alpha - 1.0e0)
 
         # M term
-        m = ((b_a - b_b) * j * a) / (alpha - 1.0e0)
+        m = ((b_cs_inner - b_b) * j_cs * r_cs_inner) / (alpha - 1.0e0)
 
         # calculate hoop stress terms
         hp_term_1 = k * ((2.0e0 + tfv.poisson_steel) / (3.0e0 * (alpha + 1.0e0)))

From 1fb3849df1c6f697324d1a09fd384002c074c601 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 10:56:35 +0100
Subject: [PATCH 03/26] Refactor CS hoop stress calculation to include
 additional parameters for improved accuracy

---
 .../source/eng-models/central-solenoid.md     | 14 +++++++++--
 process/models/pfcoil.py                      | 25 ++++++++++++-------
 2 files changed, 28 insertions(+), 11 deletions(-)

diff --git a/documentation/source/eng-models/central-solenoid.md b/documentation/source/eng-models/central-solenoid.md
index 4a62e1e6c4..6ab0b083e6 100644
--- a/documentation/source/eng-models/central-solenoid.md
+++ b/documentation/source/eng-models/central-solenoid.md
@@ -167,8 +167,8 @@ Wilson (1983).  This is divided by the fraction of the area occupied by steel to
 stress in the steel, $\sigma_{hoop}$.
 
 $$
-\sigma_{\theta} = \frac{K(2+v)}{3(\alpha+1)}\times \left(\alpha^2+\alpha+1+\frac{\alpha^2}{\eps^2}-\eps \frac{(1+2v)(\alpha+1)}{(2+v)}\right) \\
-- \frac{M(3+v)}{8}\left(\alpha^2+1+\frac{\alpha^2}{\eps^2}-\frac{(1+3v)}{(3+v)}\eps^2\right)
+\sigma_{\theta} = \frac{K(2+v)}{3(\alpha+1)}\times \left(\alpha^2+\alpha+1+\frac{\alpha^2}{\epsilon^2}-\epsilon \frac{(1+2v)(\alpha+1)}{(2+v)}\right) \\
+- \frac{M(3+v)}{8}\left(\alpha^2+1+\frac{\alpha^2}{\epsilon^2}-\frac{(1+3v)}{(3+v)}\epsilon^2\right)
 $$
 
 Where:
@@ -176,6 +176,16 @@ Where:
 - $\epsilon = \frac{r}{r_{CS,inner}}$
 - $\alpha = \frac{r_{CS,outer}}{r_{CS,inner}}$
 
+The terms $K$ and $M$ are given by:
+
+$$
+K = \frac{Ja\left(\alpha B_{\text{a}} - B_{\text{b}}\right)}{(\alpha-1)}
+$$
+
+$$
+M = \frac{Ja\left(B_{\text{a}} - B_{\text{b}}\right)}{(\alpha-1)}
+$$
+
 For a special infinite solenoid case $B$
 
 The axial stress can be calculated using "Case studies in superconducting magnets", Y. Iwasa, p. 
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 9b784a4992..756261256f 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -3340,8 +3340,18 @@ def ohcalc(self):
             # Superconducting coil
 
             # New calculation from M. N. Wilson for hoop stress
-            self.data.pf_coil.sig_hoop = self.calculate_cs_hoop_stress(
-                self.data.pf_coil.r_pf_coil_inner[self.data.pf_coil.n_cs_pf_coils - 1]
+            pfcoil_variables.sig_hoop = self.calculate_cs_hoop_stress(
+                r_stress_point=pfcoil_variables.r_pf_coil_inner[
+                    pfcoil_variables.n_cs_pf_coils - 1
+                ],
+                r_cs_inner=pfcoil_variables.r_pf_coil_inner[
+                    pfcoil_variables.n_cs_pf_coils - 1
+                ],
+                r_cs_outer=pfcoil_variables.r_pf_coil_outer[
+                    pfcoil_variables.n_cs_pf_coils - 1
+                ],
+                j_cs=pfcoil_variables.j_cs_pulse_start,
+                b_cs_inner=pfcoil_variables.b_cs_peak_pulse_start,
             )
 
             # New calculation from Y. Iwasa for axial stress
@@ -3818,7 +3828,6 @@ def calculate_cs_hoop_stress(
         r_stress_point,
         r_cs_inner=None,
         r_cs_outer=None,
-        dz_cs_half=None,
         j_cs=None,
         b_cs_inner=None,
     ) -> float:
@@ -3832,17 +3841,16 @@ def calculate_cs_hoop_stress(
         r : float
             radial position r_cs_inner < r < b
 
-        
-        
-        
-        
+
+
+
+
         Returns
         -------
         float
             hoop stress (MPa)
         """
 
-        beta = dz_cs_half / r_cs_inner
         alpha = r_cs_outer / r_cs_inner
 
         # alpha
@@ -3851,7 +3859,6 @@ def calculate_cs_hoop_stress(
         # epsilon
         epsilon = r_stress_point / r_cs_inner
 
-
         # Field at outer radius of coil [T]
         # Assume to be 0 for now
         b_b = 0.0e0

From ba3a4352a5abca55a4bafce4fb778254f2678e91 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 11:01:18 +0100
Subject: [PATCH 04/26] Rename sig_hoop variable to stress_hoop_cs for clarity
 and update references

---
 process/data_structure/pfcoil_variables.py |  2 +-
 process/models/pfcoil.py                   | 22 +++++++++++-----------
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
index 720e813ae5..d82a813e4f 100644
--- a/process/data_structure/pfcoil_variables.py
+++ b/process/data_structure/pfcoil_variables.py
@@ -42,7 +42,7 @@ class PFCoilData:
     )
     """Axial stress (z) in central solenoid at midplane due to its own field at each time point (Pa)"""
 
-    sig_hoop: float = 0.0
+    stress_hoop_cs: float = 0.0
 
     forc_z_cs_self_peak_midplane: float = 0.0
     """Axial force (z) on central solenoid at midplane due to its own field (when at peak current) (N)"""
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 756261256f..5102d5df83 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -2199,8 +2199,8 @@ def outpf(self):
             op.ovarre(
                 self.outfile,
                 "Hoop stress in CS steel (Pa)",
-                "(sig_hoop)",
-                self.data.pf_coil.sig_hoop,
+                "(stress_hoop_cs)",
+                pfcoil_variables.stress_hoop_cs,
                 "OP ",
             )
             op.ovarre(
@@ -3340,7 +3340,7 @@ def ohcalc(self):
             # Superconducting coil
 
             # New calculation from M. N. Wilson for hoop stress
-            pfcoil_variables.sig_hoop = self.calculate_cs_hoop_stress(
+            pfcoil_variables.stress_hoop_cs = self.calculate_cs_hoop_stress(
                 r_stress_point=pfcoil_variables.r_pf_coil_inner[
                     pfcoil_variables.n_cs_pf_coils - 1
                 ],
@@ -3383,7 +3383,7 @@ def ohcalc(self):
                     self.data.cs_fatigue.n_cycle,
                     self.data.cs_fatigue.t_crack_radial,
                 ) = self.cs_fatigue.ncycle(
-                    self.data.pf_coil.sig_hoop,
+                    pfcoil_variables.stress_hoop_cs,
                     self.data.cs_fatigue.residual_sig_hoop,
                     self.data.cs_fatigue.t_crack_vertical,
                     self.data.cs_fatigue.dz_cs_turn_conduit,
@@ -3401,17 +3401,17 @@ def ohcalc(self):
             if self.data.pf_coil.i_cs_stress == 1:
                 self.data.pf_coil.s_shear_cs_peak = max(
                     abs(
-                        self.data.pf_coil.sig_hoop
-                        - self.data.pf_coil.stress_z_cs_self_peak_midplane
+                        pfcoil_variables.stress_hoop_cs
+                        - pfcoil_variables.stress_z_cs_self_peak_midplane
                     ),
-                    abs(self.data.pf_coil.stress_z_cs_self_peak_midplane - 0.0e0),
-                    abs(0.0e0 - self.data.pf_coil.sig_hoop),
+                    abs(pfcoil_variables.stress_z_cs_self_peak_midplane - 0.0e0),
+                    abs(0.0e0 - pfcoil_variables.stress_hoop_cs),
                 )
             else:
-                self.data.pf_coil.s_shear_cs_peak = max(
-                    abs(self.data.pf_coil.sig_hoop - 0.0e0),
+                pfcoil_variables.s_shear_cs_peak = max(
+                    abs(pfcoil_variables.stress_hoop_cs - 0.0e0),
                     abs(0.0e0 - 0.0e0),
-                    abs(0.0e0 - self.data.pf_coil.sig_hoop),
+                    abs(0.0e0 - pfcoil_variables.stress_hoop_cs),
                 )
 
             # Thickness of hypothetical steel cylinders assumed to encase the CS along

From 5cc5362755b60818f5f822cc9bb61e65b13ae18f Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 11:06:22 +0100
Subject: [PATCH 05/26] Update CSCoil docstring to clarify parameters and add
 reference for hoop stress calculation

---
 process/models/pfcoil.py | 23 ++++++++++++++++-------
 1 file changed, 16 insertions(+), 7 deletions(-)

diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 5102d5df83..7111a4e413 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -3838,17 +3838,26 @@ def calculate_cs_hoop_stress(
 
         Parameters
         ----------
-        r : float
-            radial position r_cs_inner < r < b
-
-
-
-
+        r_stress_point : float
+            Radial location at which to calculate the hoop stress (m)
+        r_cs_inner : float
+            Inner radius of the central solenoid (m)
+        r_cs_outer : float
+            Outer radius of the central solenoid (m)
+        j_cs : float
+            Current density in the central solenoid (A/m^2)
+        b_cs_inner : float
+            Magnetic field at the inner radius of the central solenoid (T)
 
         Returns
         -------
         float
-            hoop stress (MPa)
+            hoop stress at the specified radial location (MPa)
+            
+        References
+        ----------
+        - M. N. Wilson, Superconducting Magnets. Oxford University Press, USA, 1983.
+‌
         """
 
         alpha = r_cs_outer / r_cs_inner

From d53cf4714bc896b0f6313fbee77081cffd0193c6 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 11:09:31 +0100
Subject: [PATCH 06/26] Rename stress_hoop_cs variable to stress_hoop_cs_inner
 for clarity and update references

---
 process/data_structure/pfcoil_variables.py |  4 +-
 process/models/pfcoil.py                   | 66 +++++++++++-----------
 2 files changed, 36 insertions(+), 34 deletions(-)

diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
index d82a813e4f..b08b31938a 100644
--- a/process/data_structure/pfcoil_variables.py
+++ b/process/data_structure/pfcoil_variables.py
@@ -42,7 +42,9 @@ class PFCoilData:
     )
     """Axial stress (z) in central solenoid at midplane due to its own field at each time point (Pa)"""
 
-    stress_hoop_cs: float = 0.0
+    stress_hoop_cs_inner: float = 0.0
+"""Hoop stress in central solenoid at inboard edge due to its own field (when at peak current) (Pa)"""
+
 
     forc_z_cs_self_peak_midplane: float = 0.0
     """Axial force (z) on central solenoid at midplane due to its own field (when at peak current) (N)"""
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 7111a4e413..1c9fb921b8 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -2199,8 +2199,8 @@ def outpf(self):
             op.ovarre(
                 self.outfile,
                 "Hoop stress in CS steel (Pa)",
-                "(stress_hoop_cs)",
-                pfcoil_variables.stress_hoop_cs,
+                "(stress_hoop_cs_inner)",
+                pfcoil_variables.stress_hoop_cs_inner,
                 "OP ",
             )
             op.ovarre(
@@ -3340,7 +3340,7 @@ def ohcalc(self):
             # Superconducting coil
 
             # New calculation from M. N. Wilson for hoop stress
-            pfcoil_variables.stress_hoop_cs = self.calculate_cs_hoop_stress(
+            pfcoil_variables.stress_hoop_cs_inner = self.calculate_cs_hoop_stress(
                 r_stress_point=pfcoil_variables.r_pf_coil_inner[
                     pfcoil_variables.n_cs_pf_coils - 1
                 ],
@@ -3383,7 +3383,7 @@ def ohcalc(self):
                     self.data.cs_fatigue.n_cycle,
                     self.data.cs_fatigue.t_crack_radial,
                 ) = self.cs_fatigue.ncycle(
-                    pfcoil_variables.stress_hoop_cs,
+                    pfcoil_variables.stress_hoop_cs_inner,
                     self.data.cs_fatigue.residual_sig_hoop,
                     self.data.cs_fatigue.t_crack_vertical,
                     self.data.cs_fatigue.dz_cs_turn_conduit,
@@ -3401,17 +3401,17 @@ def ohcalc(self):
             if self.data.pf_coil.i_cs_stress == 1:
                 self.data.pf_coil.s_shear_cs_peak = max(
                     abs(
-                        pfcoil_variables.stress_hoop_cs
+                        pfcoil_variables.stress_hoop_cs_inner
                         - pfcoil_variables.stress_z_cs_self_peak_midplane
                     ),
                     abs(pfcoil_variables.stress_z_cs_self_peak_midplane - 0.0e0),
-                    abs(0.0e0 - pfcoil_variables.stress_hoop_cs),
+                    abs(0.0e0 - pfcoil_variables.stress_hoop_cs_inner),
                 )
             else:
                 pfcoil_variables.s_shear_cs_peak = max(
-                    abs(pfcoil_variables.stress_hoop_cs - 0.0e0),
+                    abs(pfcoil_variables.stress_hoop_cs_inner - 0.0e0),
                     abs(0.0e0 - 0.0e0),
-                    abs(0.0e0 - pfcoil_variables.stress_hoop_cs),
+                    abs(0.0e0 - pfcoil_variables.stress_hoop_cs_inner),
                 )
 
             # Thickness of hypothetical steel cylinders assumed to encase the CS along
@@ -3833,31 +3833,31 @@ def calculate_cs_hoop_stress(
     ) -> float:
         """Calculation of hoop stress of central solenoid.
 
-        This routine calculates the hoop stress of the central solenoid
-        from "Superconducting magnets", M. N. Wilson OUP
-
-        Parameters
-        ----------
-        r_stress_point : float
-            Radial location at which to calculate the hoop stress (m)
-        r_cs_inner : float
-            Inner radius of the central solenoid (m)
-        r_cs_outer : float
-            Outer radius of the central solenoid (m)
-        j_cs : float
-            Current density in the central solenoid (A/m^2)
-        b_cs_inner : float
-            Magnetic field at the inner radius of the central solenoid (T)
-
-        Returns
-        -------
-        float
-            hoop stress at the specified radial location (MPa)
-            
-        References
-        ----------
-        - M. N. Wilson, Superconducting Magnets. Oxford University Press, USA, 1983.
-‌
+                This routine calculates the hoop stress of the central solenoid
+                from "Superconducting magnets", M. N. Wilson OUP
+
+                Parameters
+                ----------
+                r_stress_point : float
+                    Radial location at which to calculate the hoop stress (m)
+                r_cs_inner : float
+                    Inner radius of the central solenoid (m)
+                r_cs_outer : float
+                    Outer radius of the central solenoid (m)
+                j_cs : float
+                    Current density in the central solenoid (A/m^2)
+                b_cs_inner : float
+                    Magnetic field at the inner radius of the central solenoid (T)
+
+                Returns
+                -------
+                float
+                    hoop stress at the specified radial location (MPa)
+
+                References
+                ----------
+                - M. N. Wilson, Superconducting Magnets. Oxford University Press, USA, 1983.
+        ‌
         """
 
         alpha = r_cs_outer / r_cs_inner

From d425f4a591070f54264aed8f609350701a52024c Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 11:10:07 +0100
Subject: [PATCH 07/26] Add stress_hoop_cs_inner_profile variable and update
 initialization in PF coil module

---
 process/data_structure/pfcoil_variables.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
index b08b31938a..45ac94ac0f 100644
--- a/process/data_structure/pfcoil_variables.py
+++ b/process/data_structure/pfcoil_variables.py
@@ -45,6 +45,8 @@ class PFCoilData:
     stress_hoop_cs_inner: float = 0.0
 """Hoop stress in central solenoid at inboard edge due to its own field (when at peak current) (Pa)"""
 
+stress_hoop_cs_inner_profile: list[float] = None
+"""Hoop stress in central solenoid at inboard edge due to its own field at each time point (Pa)"""
 
     forc_z_cs_self_peak_midplane: float = 0.0
     """Axial force (z) on central solenoid at midplane due to its own field (when at peak current) (N)"""

From 6b7356068ff1fc7ea8ee0464c51e6442784b20b4 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 11:47:24 +0100
Subject: [PATCH 08/26] Add a_cs_steel_poloidal variable and update
 initialization in PF coil module

---
 process/data_structure/pfcoil_variables.py |  3 +++
 process/models/pfcoil.py                   | 20 ++++++++++----------
 2 files changed, 13 insertions(+), 10 deletions(-)

diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
index 45ac94ac0f..c6f6c557ad 100644
--- a/process/data_structure/pfcoil_variables.py
+++ b/process/data_structure/pfcoil_variables.py
@@ -107,6 +107,9 @@ class PFCoilData:
     a_cs_poloidal: float = 0.0
     """Central solenoid vertical cross-sectional area (m2)"""
 
+a_cs_steel_poloidal: float = None
+"""Central solenoid vertical cross-sectional area of steel (m2)"""
+
     a_cs_turn: float = 0.0
     """Central solenoid (OH) trun cross-sectional area (m2)"""
 
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 1c9fb921b8..f1f42fb0a9 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -3267,13 +3267,13 @@ def ohcalc(self):
         )
 
         # Peak field at the End-Of-Flattop (EOF)
-        # Occurs at inner edge of coil; bmaxoh2 and bzi are of opposite sign at EOF
+        # Occurs at inner edge of coil; b_cs_self_peak_flat_top_end and bzi are of opposite sign at EOF
 
         # Peak field due to central Solenoid itself
-        bmaxoh2 = self.calculate_cs_self_peak_magnetic_field(
-            j_cs=self.data.pf_coil.j_cs_flat_top_end,
-            r_cs_inner=self.data.pf_coil.r_pf_coil_inner[
-                self.data.pf_coil.n_cs_pf_coils - 1
+        b_cs_self_peak_flat_top_end = self.calculate_cs_self_peak_magnetic_field(
+            j_cs=pfcoil_variables.j_cs_flat_top_end,
+            r_cs_inner=pfcoil_variables.r_pf_coil_inner[
+                pfcoil_variables.n_cs_pf_coils - 1
             ],
             r_cs_outer=self.data.pf_coil.r_pf_coil_outer[
                 self.data.pf_coil.n_cs_pf_coils - 1
@@ -3285,14 +3285,14 @@ def ohcalc(self):
 
         # Peak field due to other PF coils plus plasma
         timepoint = 5
-        _bri, _bro, bzi, bzo = peak_b_field_at_pf_coil(
-            n_coil=self.data.pf_coil.n_cs_pf_coils,
+        _, _, bzi, bzo = peak_b_field_at_pf_coil(
+            n_coil=pfcoil_variables.n_cs_pf_coils,
             n_coil_group=99,
             t_b_field_peak=timepoint,
             data=self.data,
         )
 
-        self.data.pf_coil.b_cs_peak_flat_top_end = abs(bzi - bmaxoh2)
+        pfcoil_variables.b_cs_peak_flat_top_end = abs(bzi - b_cs_self_peak_flat_top_end)
 
         # Peak field on outboard side of central Solenoid
         # (self-field is assumed to be zero - long solenoid approximation)
@@ -3315,8 +3315,8 @@ def ohcalc(self):
             )
         )
         timepoint = 2
-        _bri, _bro, bzi, bzo = peak_b_field_at_pf_coil(
-            n_coil=self.data.pf_coil.n_cs_pf_coils,
+        _, _, bzi, bzo = peak_b_field_at_pf_coil(
+            n_coil=pfcoil_variables.n_cs_pf_coils,
             n_coil_group=99,
             t_b_field_peak=timepoint,
             data=self.data,

From de061aa5941377225afcb544766112f60c2c94dc Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 11:52:13 +0100
Subject: [PATCH 09/26] Refactor CSCoil to calculate a_cs_steel_poloidal
 directly and update related variables

---
 process/models/pfcoil.py | 22 +++++++++++++---------
 1 file changed, 13 insertions(+), 9 deletions(-)

diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index f1f42fb0a9..33a2765dde 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -3394,8 +3394,8 @@ def ohcalc(self):
             # equation is used for Central Solenoid stress
 
             # Area of steel in Central Solenoid
-            areaspf = (
-                self.data.pf_coil.f_a_cs_turn_steel * self.data.pf_coil.a_cs_poloidal
+            pfcoil_variables.a_cs_steel_poloidal = (
+                pfcoil_variables.f_a_cs_turn_steel * pfcoil_variables.a_cs_poloidal
             )
 
             if self.data.pf_coil.i_cs_stress == 1:
@@ -3419,17 +3419,19 @@ def ohcalc(self):
             # throughout the conductor
             self.data.pf_coil.pfcaseth[self.data.pf_coil.n_cs_pf_coils - 1] = (
                 0.25e0
-                * areaspf
-                / self.data.pf_coil.z_pf_coil_upper[self.data.pf_coil.n_cs_pf_coils - 1]
+                * pfcoil_variables.a_cs_steel_poloidal
+                / pfcoil_variables.z_pf_coil_upper[pfcoil_variables.n_cs_pf_coils - 1]
             )
 
         else:
-            areaspf = 0.0e0  # Resistive Central Solenoid - no steel needed
-            self.data.pf_coil.pfcaseth[self.data.pf_coil.n_cs_pf_coils - 1] = 0.0e0
+            pfcoil_variables.a_cs_steel_poloidal = (
+                0.0e0  # Resistive Central Solenoid - no steel needed
+            )
+            pfcoil_variables.pfcaseth[pfcoil_variables.n_cs_pf_coils - 1] = 0.0e0
 
         # Weight of steel
-        self.data.pf_coil.m_pf_coil_structure[self.data.pf_coil.n_cs_pf_coils - 1] = (
-            areaspf
+        pfcoil_variables.m_pf_coil_structure[pfcoil_variables.n_cs_pf_coils - 1] = (
+            pfcoil_variables.a_cs_steel_poloidal
             * 2.0e0
             * np.pi
             * self.data.pf_coil.r_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1]
@@ -3437,7 +3439,9 @@ def ohcalc(self):
         )
 
         # Non-steel cross-sectional area
-        self.data.pf_coil.awpoh = self.data.pf_coil.a_cs_poloidal - areaspf
+        pfcoil_variables.awpoh = (
+            pfcoil_variables.a_cs_poloidal - pfcoil_variables.a_cs_steel_poloidal
+        )
 
         # Issue #97. Fudge to ensure awpoh is positive; result is continuous, smooth and
         # monotonically decreases

From 0e7d5ae081c2aaf87da575c7972800011e85eb93 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 11:58:44 +0100
Subject: [PATCH 10/26] Update output to reference a_cs_steel_poloidal in PF
 coil module

---
 process/models/pfcoil.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 33a2765dde..efda5e8c72 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -2169,8 +2169,8 @@ def outpf(self):
             op.ovarre(
                 self.outfile,
                 "CS steel cross-sectional area (m2)",
-                "(a_cs_poloidal-awpoh)",
-                self.data.pf_coil.a_cs_poloidal - self.data.pf_coil.awpoh,
+                "(a_cs_steel_poloidal)",
+                pfcoil_variables.a_cs_steel_poloidal,
                 "OP ",
             )
             op.ovarre(

From d718cabdf58da9f9c4c3bd62f6fdcd6f9d9f920f Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 13:03:38 +0100
Subject: [PATCH 11/26] Refactor variable names from awpoh to a_cs_cable_space
 across multiple modules and update related references

---
 process/data_structure/pfcoil_variables.py |  2 +-
 process/models/costs/costs.py              | 34 ++++++-------
 process/models/pfcoil.py                   | 57 ++++++++++++----------
 tests/unit/models/test_costs_1990.py       | 12 +++--
 4 files changed, 55 insertions(+), 50 deletions(-)

diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
index c6f6c557ad..ddc628334a 100644
--- a/process/data_structure/pfcoil_variables.py
+++ b/process/data_structure/pfcoil_variables.py
@@ -113,7 +113,7 @@ class PFCoilData:
     a_cs_turn: float = 0.0
     """Central solenoid (OH) trun cross-sectional area (m2)"""
 
-    awpoh: float = 0.0
+    a_cs_cable_space: float = 0.0
     """central solenoid conductor+void area with area of steel subtracted (m2)"""
 
     b_cs_peak_flat_top_end: float = 0.0
diff --git a/process/models/costs/costs.py b/process/models/costs/costs.py
index bc68a495e3..f399887d3f 100644
--- a/process/models/costs/costs.py
+++ b/process/models/costs/costs.py
@@ -1715,12 +1715,12 @@ def acc2222(self):
                 #  Issue #328  Use CS conductor cross-sectional area (m2)
                 if self.data.pf_coil.i_pf_conductor == 0:
                     costpfsc = (
-                        self.data.costs.ucsc[self.data.pf_coil.i_cs_superconductor - 1]
-                        * self.data.pf_coil.awpoh
-                        * (1 - self.data.pf_coil.f_a_cs_void)
-                        * (1 - self.data.pf_coil.fcuohsu)
-                        / self.data.pf_coil.n_pf_coil_turns[
-                            self.data.pf_coil.n_cs_pf_coils - 1
+                        cost_variables.ucsc[pfcoil_variables.i_cs_superconductor - 1]
+                        * pfcoil_variables.a_cs_cable_space
+                        * (1 - pfcoil_variables.f_a_cs_void)
+                        * (1 - pfcoil_variables.fcuohsu)
+                        / pfcoil_variables.n_pf_coil_turns[
+                            pfcoil_variables.n_cs_pf_coils - 1
                         ]
                         * tfcoil_variables.dcond[
                             self.data.pf_coil.i_cs_superconductor - 1
@@ -1745,23 +1745,23 @@ def acc2222(self):
 
             if self.data.pf_coil.i_pf_conductor == 0:
                 costpfcu = (
-                    self.data.costs.uccu
-                    * self.data.pf_coil.awpoh
-                    * (1 - self.data.pf_coil.f_a_cs_void)
-                    * self.data.pf_coil.fcuohsu
-                    / self.data.pf_coil.n_pf_coil_turns[
-                        self.data.pf_coil.n_cs_pf_coils - 1
+                    cost_variables.uccu
+                    * pfcoil_variables.a_cs_cable_space
+                    * (1 - pfcoil_variables.f_a_cs_void)
+                    * pfcoil_variables.fcuohsu
+                    / pfcoil_variables.n_pf_coil_turns[
+                        pfcoil_variables.n_cs_pf_coils - 1
                     ]
                     * constants.den_copper
                 )
             else:
                 # MDK I don't know if this is ccorrect as we never use the resistive model
                 costpfcu = (
-                    self.data.costs.uccu
-                    * self.data.pf_coil.awpoh
-                    * (1 - self.data.pf_coil.f_a_cs_void)
-                    / self.data.pf_coil.n_pf_coil_turns[
-                        self.data.pf_coil.n_cs_pf_coils - 1
+                    cost_variables.uccu
+                    * pfcoil_variables.a_cs_cable_space
+                    * (1 - pfcoil_variables.f_a_cs_void)
+                    / pfcoil_variables.n_pf_coil_turns[
+                        pfcoil_variables.n_cs_pf_coils - 1
                     ]
                     * constants.den_copper
                 )
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index efda5e8c72..868db55aa8 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -2148,22 +2148,23 @@ def outpf(self):
             op.ovarre(
                 self.outfile,
                 "CS conductor+void cross-sectional area (m2)",
-                "(awpoh)",
-                self.data.pf_coil.awpoh,
+                "(a_cs_cable_space)",
+                pfcoil_variables.a_cs_cable_space,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
                 "   CS conductor cross-sectional area (m2)",
-                "(awpoh*(1-f_a_cs_void))",
-                self.data.pf_coil.awpoh * (1.0e0 - self.data.pf_coil.f_a_cs_void),
+                "(a_cs_cable_space*(1-f_a_cs_void))",
+                pfcoil_variables.a_cs_cable_space
+                    * (1.0e0 - pfcoil_variables.f_a_cs_void),
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
                 "   CS void cross-sectional area (m2)",
-                "(awpoh*f_a_cs_void)",
-                self.data.pf_coil.awpoh * self.data.pf_coil.f_a_cs_void,
+                "(a_cs_cable_space*f_a_cs_void)",
+                pfcoil_variables.a_cs_cable_space * pfcoil_variables.f_a_cs_void,
                 "OP ",
             )
             op.ovarre(
@@ -3439,35 +3440,33 @@ def ohcalc(self):
         )
 
         # Non-steel cross-sectional area
-        pfcoil_variables.awpoh = (
+        pfcoil_variables.a_cs_cable_space = (
             pfcoil_variables.a_cs_poloidal - pfcoil_variables.a_cs_steel_poloidal
         )
 
-        # Issue #97. Fudge to ensure awpoh is positive; result is continuous, smooth and
+        # Issue #97. Fudge to ensure a_cs_cable_space is positive; result is continuous, smooth and
         # monotonically decreases
 
         da = 0.0001e0  # 1 cm^2
-        if self.data.pf_coil.awpoh < da:
-            self.data.pf_coil.awpoh = da * da / (2.0e0 * da - self.data.pf_coil.awpoh)
+        if pfcoil_variables.a_cs_cable_space < da:
+            pfcoil_variables.a_cs_cable_space = (
+                da * da / (2.0e0 * da - pfcoil_variables.a_cs_cable_space)
+            )
 
         # Weight of conductor in central Solenoid
-        if self.data.pf_coil.i_pf_conductor == 0:
-            self.data.pf_coil.m_pf_coil_conductor[
-                self.data.pf_coil.n_cs_pf_coils - 1
-            ] = (
-                self.data.pf_coil.awpoh
-                * (1.0e0 - self.data.pf_coil.f_a_cs_void)
+        if pfcoil_variables.i_pf_conductor == 0:
+            pfcoil_variables.m_pf_coil_conductor[pfcoil_variables.n_cs_pf_coils - 1] = (
+                pfcoil_variables.a_cs_cable_space
+                * (1.0e0 - pfcoil_variables.f_a_cs_void)
                 * 2.0e0
                 * np.pi
                 * self.data.pf_coil.r_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1]
                 * tfv.dcond[self.data.pf_coil.i_cs_superconductor - 1]
             )
         else:
-            self.data.pf_coil.m_pf_coil_conductor[
-                self.data.pf_coil.n_cs_pf_coils - 1
-            ] = (
-                self.data.pf_coil.awpoh
-                * (1.0e0 - self.data.pf_coil.f_a_cs_void)
+            pfcoil_variables.m_pf_coil_conductor[pfcoil_variables.n_cs_pf_coils - 1] = (
+                pfcoil_variables.a_cs_cable_space
+                * (1.0e0 - pfcoil_variables.f_a_cs_void)
                 * 2.0e0
                 * np.pi
                 * self.data.pf_coil.r_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1]
@@ -3493,7 +3492,7 @@ def ohcalc(self):
                             self.data.pf_coil.n_cs_pf_coils - 1
                         ]
                     )
-                    / self.data.pf_coil.awpoh
+                    / pfcoil_variables.a_cs_cable_space
                 )
                 * 1.0e6,
                 self.data.pf_coil.i_cs_superconductor,
@@ -3515,8 +3514,10 @@ def ohcalc(self):
                     * (1 - self.data.pf_coil.fcuohsu)
                 )
 
-            self.data.pf_coil.j_cs_critical_flat_top_end = (
-                jcritwp * self.data.pf_coil.awpoh / self.data.pf_coil.a_cs_poloidal
+            pfcoil_variables.j_cs_critical_flat_top_end = (
+                jcritwp
+                * pfcoil_variables.a_cs_cable_space
+                / pfcoil_variables.a_cs_poloidal
             )
 
             # Allowable coil overall current density at BOP
@@ -3536,7 +3537,7 @@ def ohcalc(self):
                             self.data.pf_coil.n_cs_pf_coils - 1
                         ]
                     )
-                    / self.data.pf_coil.awpoh
+                    / pfcoil_variables.a_cs_cable_space
                 )
                 * 1.0e6,
                 self.data.pf_coil.i_cs_superconductor,
@@ -3547,8 +3548,10 @@ def ohcalc(self):
                 tfv.tcritsc,
             )
 
-            self.data.pf_coil.j_pf_wp_critical[self.data.pf_coil.n_cs_pf_coils - 1] = (
-                jcritwp * self.data.pf_coil.awpoh / self.data.pf_coil.a_cs_poloidal
+            pfcoil_variables.j_pf_wp_critical[pfcoil_variables.n_cs_pf_coils - 1] = (
+                jcritwp
+                * pfcoil_variables.a_cs_cable_space
+                / pfcoil_variables.a_cs_poloidal
             )
             self.data.pf_coil.j_cs_critical_pulse_start = (
                 self.data.pf_coil.j_pf_wp_critical[self.data.pf_coil.n_cs_pf_coils - 1]
diff --git a/tests/unit/models/test_costs_1990.py b/tests/unit/models/test_costs_1990.py
index 164c7c4145..e11834ce01 100644
--- a/tests/unit/models/test_costs_1990.py
+++ b/tests/unit/models/test_costs_1990.py
@@ -2058,7 +2058,7 @@ class Acc2222Param(NamedTuple):
 
     f_a_pf_coil_void: Any = None
 
-    awpoh: Any = None
+    a_cs_cable_space: Any = None
 
     fncmass: Any = None
 
@@ -2265,7 +2265,7 @@ class Acc2222Param(NamedTuple):
                 ),
                 order="F",
             ).transpose(),
-            awpoh=3.8004675824985918,
+            a_cs_cable_space=3.8004675824985918,
             fncmass=310716.52923547616,
             dcond=np.array(
                 np.array(
@@ -2461,7 +2461,7 @@ class Acc2222Param(NamedTuple):
                 ),
                 order="F",
             ).transpose(),
-            awpoh=3.8004675824985918,
+            a_cs_cable_space=3.8004675824985918,
             fncmass=310716.52923547616,
             dcond=np.array(
                 np.array(
@@ -2657,7 +2657,7 @@ class Acc2222Param(NamedTuple):
                 ),
                 order="F",
             ).transpose(),
-            awpoh=3.8004675824985918,
+            a_cs_cable_space=3.8004675824985918,
             fncmass=310716.52923547616,
             dcond=np.array(
                 np.array(
@@ -2767,7 +2767,9 @@ def test_acc2222(acc2222param, monkeypatch, costs):
         costs.data.pf_coil, "f_a_pf_coil_void", acc2222param.f_a_pf_coil_void
     )
 
-    monkeypatch.setattr(costs.data.pf_coil, "awpoh", acc2222param.awpoh)
+    monkeypatch.setattr(
+        pfcoil_variables, "a_cs_cable_space", acc2222param.a_cs_cable_space
+    )
 
     monkeypatch.setattr(costs.data.structure, "fncmass", acc2222param.fncmass)
 

From 6e4a21a430a0a302c4f8a584c90fd27fd45c4748 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 13:15:06 +0100
Subject: [PATCH 12/26] Add inner and outer radius variables for central
 solenoid in PF coil module

---
 process/data_structure/pfcoil_variables.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
index ddc628334a..6713426f59 100644
--- a/process/data_structure/pfcoil_variables.py
+++ b/process/data_structure/pfcoil_variables.py
@@ -354,6 +354,12 @@ class PFCoilData:
     r_cs_middle: float = 0.0
     """radius to the centre of the central solenoid (m)"""
 
+r_cs_inner: float = None
+"""inner radius of the central solenoid (m)"""
+
+r_cs_outer: float = None
+"""outer radius of the central solenoid (m)"""
+
     dz_cs_full: float = 0.0
     """Full height of the central solenoid (m)"""
 

From e58868004156b631edbc4dc333c5f0224636d38a Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 13:17:44 +0100
Subject: [PATCH 13/26] Add radial inner and outer parameters for central
 solenoid in PF coil module

---
 process/models/pfcoil.py | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 868db55aa8..b92f9d0d19 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -2145,6 +2145,20 @@ def outpf(self):
                 self.data.pf_coil.r_cs_middle,
                 "OP ",
             )
+                op.ovarre(
+                    self.outfile,
+                    "CS radial inner (m)",
+                    "(r_cs_inner)",
+                    pfcoil_variables.r_cs_inner,
+                    "OP ",
+                )
+                op.ovarre(
+                    self.outfile,
+                    "CS radial outer (m)",
+                    "(r_cs_outer)",
+                    pfcoil_variables.r_cs_outer,
+                    "OP ",
+                )
             op.ovarre(
                 self.outfile,
                 "CS conductor+void cross-sectional area (m2)",
@@ -3208,6 +3222,13 @@ def ohcalc(self):
             dr_bore=self.data.build.dr_bore,
         )
 
+        pfcoil_variables.r_cs_inner = pfcoil_variables.r_pf_coil_inner[
+            pfcoil_variables.n_cs_pf_coils - 1
+        ]
+        pfcoil_variables.r_cs_outer = pfcoil_variables.r_pf_coil_outer[
+            pfcoil_variables.n_cs_pf_coils - 1
+        ]
+
         # Maximum current (MA-turns) in central Solenoid, at either BOP or EOF
         if self.data.pf_coil.j_cs_pulse_start > self.data.pf_coil.j_cs_flat_top_end:
             sgn = 1.0e0

From 234e5104f53d74da29181f6a3b132470d0344fb8 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 13:48:02 +0100
Subject: [PATCH 14/26] Add radial hoop stress profile plotting to CSCoil class
 and update parameter types

---
 process/core/io/plot/summary.py |  16 ++++-
 process/models/pfcoil.py        | 102 +++++++++++++++++++++++++++++---
 2 files changed, 109 insertions(+), 9 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index ac1e0951d5..d00fd4f63f 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -24,6 +24,7 @@
 from process.data_structure import impurity_radiation_module
 from process.data_structure.pfcoil_variables import NFIXMX
 from process.models.build import Build
+from process.models.cs_fatigue import CsFatigue
 from process.models.geometry.blanket import (
     blanket_geometry_double_null,
     blanket_geometry_single_null,
@@ -14496,10 +14497,21 @@ def main_plot(
 
     plot_current_profiles_over_time(figs[29].add_subplot(111), m_file, scan)
 
-    plot_cs_stress_time_profile(axis=figs[30].add_subplot(311), mfile=m_file, scan=scan)
+    CSCoil.plot_stress_time_profile(
+        axis=figs[29].add_subplot(222), mfile=m_file, scan=scan
+    )
+
+    cs_coil = CSCoil(cs_fatigue=CsFatigue())
+    cs_coil.plot_cs_radial_hoop_stress_profile(
+        axis=figs[29].add_subplot(224),
+        mfile=m_file,
+        scan=scan,
+        j_cs=m_file.get("j_cs_pulse_start", scan=scan),
+        b_cs_inner=m_file.get("b_cs_peak_pulse_start", scan=scan),
+    )
 
     plot_cs_coil_structure(
-        figs[30].add_subplot(223, aspect="equal"), figs[30], m_file, scan
+        figs[30].add_subplot(121, aspect="equal"), figs[30], m_file, scan
     )
     plot_cs_turn_structure(
         figs[30].add_subplot(326, aspect="equal"), figs[30], m_file, scan
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index b92f9d0d19..7fb8088cbe 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -3853,11 +3853,11 @@ def calculate_cs_self_midplane_axial_stress_time_profile(
 
     def calculate_cs_hoop_stress(
         self,
-        r_stress_point,
-        r_cs_inner=None,
-        r_cs_outer=None,
-        j_cs=None,
-        b_cs_inner=None,
+        r_stress_point: float | np.ndarray,
+        r_cs_inner: float,
+        r_cs_outer: float,
+        j_cs: float,
+        b_cs_inner: float,
     ) -> float:
         """Calculation of hoop stress of central solenoid.
 
@@ -3866,7 +3866,7 @@ def calculate_cs_hoop_stress(
 
                 Parameters
                 ----------
-                r_stress_point : float
+                r_stress_point : float/np.ndarray
                     Radial location at which to calculate the hoop stress (m)
                 r_cs_inner : float
                     Inner radius of the central solenoid (m)
@@ -3936,7 +3936,95 @@ def calculate_cs_hoop_stress(
 
         s_hoop_nom = hp_term_1 * hp_term_2 - hp_term_3 * hp_term_4
 
-        return s_hoop_nom / self.data.pf_coil.f_a_cs_turn_steel
+        return s_hoop_nom / pfcoil_variables.f_a_cs_turn_steel
+
+    @staticmethod
+    def plot_stress_time_profile(axis: plt.Axes, mfile: mf.MFile, scan: int):
+        t_plant_pulse_coil_precharge = mfile.get(
+            "t_plant_pulse_coil_precharge", scan=scan
+        )
+        t_plant_pulse_plasma_current_ramp_up = mfile.get(
+            "t_plant_pulse_plasma_current_ramp_up", scan=scan
+        )
+        t_plant_pulse_fusion_ramp = mfile.get("t_plant_pulse_fusion_ramp", scan=scan)
+        t_plant_pulse_burn = mfile.get("t_plant_pulse_burn", scan=scan)
+        t_plant_pulse_plasma_current_ramp_down = mfile.get(
+            "t_plant_pulse_plasma_current_ramp_down", scan=scan
+        )
+
+        # Define a cumulative sum list for each point in the pulse
+        t_steps = np.cumsum([
+            0,
+            t_plant_pulse_coil_precharge,
+            t_plant_pulse_plasma_current_ramp_up,
+            t_plant_pulse_fusion_ramp,
+            t_plant_pulse_burn,
+            t_plant_pulse_plasma_current_ramp_down,
+        ])
+
+        stress_times = t_steps[
+            :6
+        ]  # Get the first 6 time points corresponding to the stress profile
+
+        stress_z_cs_self_midplane_profile = np.zeros(6)
+        for i in range(6):
+            stress_z_cs_self_midplane_profile[i] = mfile.get(
+                f"stress_z_cs_self_midplane_profile[{i}]", scan=scan
+            )
+
+        # Plot stress vs time
+        axis.plot(
+            stress_times,
+            stress_z_cs_self_midplane_profile / 1e6,
+            "o-",
+            linewidth=2,
+            markersize=8,
+            label=r"$\sigma_{{z}}$,Midplane Axial Stress",
+        )
+        axis.set_xlabel("Pulse Time (s)")
+        axis.set_ylabel("Stress (MPa)")
+        axis.minorticks_on()
+        axis.legend(loc="best")
+        axis.set_title("Central Solenoid Stress")
+        axis.grid(True, alpha=0.3)
+
+    def plot_cs_radial_hoop_stress_profile(
+        self,
+        axis: plt.Axes,
+        mfile: mf.MFile,
+        scan: int,
+        j_cs: float,
+        b_cs_inner: float,
+    ):
+        r_cs_inner = mfile.get("r_cs_inner", scan=scan)
+        r_cs_outer = mfile.get("r_cs_outer", scan=scan)
+
+        radii = np.linspace(r_cs_inner, r_cs_outer, num=10)
+        stress_values = np.array([
+            self.calculate_cs_hoop_stress(
+                r_stress_point=radius,
+                r_cs_inner=r_cs_inner,
+                r_cs_outer=r_cs_outer,
+                j_cs=j_cs,
+                b_cs_inner=b_cs_inner,
+            )
+            for radius in radii
+        ])
+
+        axis.plot(
+            radii,
+            stress_values / 1e6,
+            "o-",
+            linewidth=2,
+            markersize=8,
+            label="$\\sigma_{\\theta}$,Hoop Stress",
+        )
+        axis.set_xlabel("Radial Position (m)")
+        axis.set_ylabel("Hoop Stress (MPa)")
+        axis.minorticks_on()
+        axis.legend(loc="best")
+        axis.set_title("Central Solenoid Radial Hoop Stress Profile")
+        axis.grid(True, alpha=0.3)
 
 
 def peak_b_field_at_pf_coil(

From 130bea54eef04b778b687125850c94f32270f626 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 15:01:33 +0100
Subject: [PATCH 15/26] Refactor CSCoil class to add radial stress profile
 calculations and update plotting functions

---
 process/core/io/plot/summary.py |  12 ++-
 process/models/pfcoil.py        | 182 ++++++++++++++++++++++++++------
 2 files changed, 162 insertions(+), 32 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index d00fd4f63f..922b0c9686 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -14498,12 +14498,20 @@ def main_plot(
     plot_current_profiles_over_time(figs[29].add_subplot(111), m_file, scan)
 
     CSCoil.plot_stress_time_profile(
-        axis=figs[29].add_subplot(222), mfile=m_file, scan=scan
+        axis=figs[29].add_subplot(322), mfile=m_file, scan=scan
     )
 
     cs_coil = CSCoil(cs_fatigue=CsFatigue())
     cs_coil.plot_cs_radial_hoop_stress_profile(
-        axis=figs[29].add_subplot(224),
+        axis=figs[29].add_subplot(324),
+        mfile=m_file,
+        scan=scan,
+        j_cs=m_file.get("j_cs_pulse_start", scan=scan),
+        b_cs_inner=m_file.get("b_cs_peak_pulse_start", scan=scan),
+    )
+
+    cs_coil.plot_cs_radial_stress_profile(
+        axis=figs[29].add_subplot(326),
         mfile=m_file,
         scan=scan,
         j_cs=m_file.get("j_cs_pulse_start", scan=scan),
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 7fb8088cbe..35c7a44752 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -3374,6 +3374,7 @@ def ohcalc(self):
                 ],
                 j_cs=pfcoil_variables.j_cs_pulse_start,
                 b_cs_inner=pfcoil_variables.b_cs_peak_pulse_start,
+                f_poisson_cs_structure=tfv.poisson_steel,
             )
 
             # New calculation from Y. Iwasa for axial stress
@@ -3851,40 +3852,44 @@ def calculate_cs_self_midplane_axial_stress_time_profile(
             )
             self.data.pf_coil.stress_z_cs_self_midplane_profile[time] = stress_value
 
+    @staticmethod
+    @numba.njit(cache=True)
     def calculate_cs_hoop_stress(
-        self,
         r_stress_point: float | np.ndarray,
         r_cs_inner: float,
         r_cs_outer: float,
         j_cs: float,
         b_cs_inner: float,
+        f_poisson_cs_structure: float,
     ) -> float:
         """Calculation of hoop stress of central solenoid.
 
-                This routine calculates the hoop stress of the central solenoid
-                from "Superconducting magnets", M. N. Wilson OUP
-
-                Parameters
-                ----------
-                r_stress_point : float/np.ndarray
-                    Radial location at which to calculate the hoop stress (m)
-                r_cs_inner : float
-                    Inner radius of the central solenoid (m)
-                r_cs_outer : float
-                    Outer radius of the central solenoid (m)
-                j_cs : float
-                    Current density in the central solenoid (A/m^2)
-                b_cs_inner : float
-                    Magnetic field at the inner radius of the central solenoid (T)
-
-                Returns
-                -------
-                float
-                    hoop stress at the specified radial location (MPa)
-
-                References
-                ----------
-                - M. N. Wilson, Superconducting Magnets. Oxford University Press, USA, 1983.
+        This routine calculates the hoop stress of the central solenoid
+        from "Superconducting magnets", M. N. Wilson OUP
+
+        Parameters
+        ----------
+        r_stress_point : float/np.ndarray
+            Radial location at which to calculate the hoop stress (m)
+        r_cs_inner : float
+            Inner radius of the central solenoid (m)
+        r_cs_outer : float
+            Outer radius of the central solenoid (m)
+        j_cs : float
+            Current density in the central solenoid (A/m^2)
+        b_cs_inner : float
+            Magnetic field at the inner radius of the central solenoid (T)
+        f_poisson_cs_structure : float
+            Poisson's ratio of the central solenoid structure (dimensionless)
+
+        Returns
+        -------
+        float
+            hoop stress at the specified radial location (MPa)
+
+        References
+        ----------
+        - M. N. Wilson, Superconducting Magnets. Oxford University Press, USA, 1983.
         ‌
         """
 
@@ -3910,7 +3915,7 @@ def calculate_cs_hoop_stress(
         m = ((b_cs_inner - b_b) * j_cs * r_cs_inner) / (alpha - 1.0e0)
 
         # calculate hoop stress terms
-        hp_term_1 = k * ((2.0e0 + tfv.poisson_steel) / (3.0e0 * (alpha + 1.0e0)))
+        hp_term_1 = k * ((2.0e0 + f_poisson_cs_structure) / (3.0e0 * (alpha + 1.0e0)))
 
         hp_term_2 = (
             alpha**2
@@ -3919,25 +3924,103 @@ def calculate_cs_hoop_stress(
             + alpha**2 / epsilon**2
             - epsilon
             * (
-                ((1.0e0 + 2.0e0 * tfv.poisson_steel) * (alpha + 1.0e0))
-                / (2.0e0 + tfv.poisson_steel)
+                ((1.0e0 + 2.0e0 * f_poisson_cs_structure) * (alpha + 1.0e0))
+                / (2.0e0 + f_poisson_cs_structure)
             )
         )
 
-        hp_term_3 = m * ((3.0e0 + tfv.poisson_steel) / (8.0e0))
+        hp_term_3 = m * ((3.0e0 + f_poisson_cs_structure) / (8.0e0))
 
         hp_term_4 = (
             alpha**2
             + 1.0e0
             + alpha**2 / epsilon**2
             - epsilon**2
-            * ((1.0e0 + 3.0e0 * tfv.poisson_steel) / (3.0e0 + tfv.poisson_steel))
+            * (
+                (1.0e0 + 3.0e0 * f_poisson_cs_structure)
+                / (3.0e0 + f_poisson_cs_structure)
+            )
         )
 
         s_hoop_nom = hp_term_1 * hp_term_2 - hp_term_3 * hp_term_4
 
         return s_hoop_nom / pfcoil_variables.f_a_cs_turn_steel
 
+    @staticmethod
+    @numba.njit(cache=True)
+    def calculate_cs_radial_stress(
+        r_stress_point: float | np.ndarray,
+        r_cs_inner: float,
+        r_cs_outer: float,
+        j_cs: float,
+        b_cs_inner: float,
+        f_poisson_cs_structure: float,
+    ) -> float:
+        """Calculation of radial stress of central solenoid.
+
+        This routine calculates the radial stress of the central solenoid
+        from "Superconducting magnets", M. N. Wilson OUP
+
+        Parameters
+        ----------
+        r_stress_point : float/np.ndarray
+            Radial location at which to calculate the hoop stress (m)
+        r_cs_inner : float
+            Inner radius of the central solenoid (m)
+        r_cs_outer : float
+            Outer radius of the central solenoid (m)
+        j_cs : float
+            Current density in the central solenoid (A/m^2)
+        b_cs_inner : float
+            Magnetic field at the inner radius of the central solenoid (T)
+        f_poisson_cs_structure : float
+            Poisson's ratio of the central solenoid structure (dimensionless)
+
+        Returns
+        -------
+        float
+            radial stress at the specified radial location (MPa)
+
+        References
+        ----------
+        - M. N. Wilson, Superconducting Magnets. Oxford University Press, USA, 1983.
+        ‌
+        """
+
+        alpha = r_cs_outer / r_cs_inner
+
+        # alpha
+        alpha = r_cs_outer / r_cs_inner
+
+        # epsilon
+        epsilon = r_stress_point / r_cs_inner
+
+        # Field at outer radius of coil [T]
+        # Assume to be 0 for now
+        b_b = 0.0e0
+
+        # current density [A/m^2]
+        j_cs = pfcoil_variables.j_cs_pulse_start
+
+        # K term
+        k = ((alpha * b_cs_inner - b_b) * j_cs * r_cs_inner) / (alpha - 1.0e0)
+
+        # M term
+        m = ((b_cs_inner - b_b) * j_cs * r_cs_inner) / (alpha - 1.0e0)
+
+        # calculate hoop stress terms
+        hp_term_1 = k * ((2.0e0 + f_poisson_cs_structure) / (3.0e0 * (alpha + 1.0e0)))
+
+        hp_term_2 = (
+            alpha**2 + alpha + 1.0e0 - alpha**2 / epsilon**2 - epsilon * (alpha + 1.0e0)
+        )
+
+        hp_term_3 = m * ((3.0e0 + f_poisson_cs_structure) / (8.0e0))
+
+        hp_term_4 = alpha**2 + 1.0e0 - alpha**2 / epsilon**2 - epsilon**2
+
+        return hp_term_1 * hp_term_2 - hp_term_3 * hp_term_4
+
     @staticmethod
     def plot_stress_time_profile(axis: plt.Axes, mfile: mf.MFile, scan: int):
         t_plant_pulse_coil_precharge = mfile.get(
@@ -4007,6 +4090,7 @@ def plot_cs_radial_hoop_stress_profile(
                 r_cs_outer=r_cs_outer,
                 j_cs=j_cs,
                 b_cs_inner=b_cs_inner,
+                f_poisson_cs_structure=tfv.poisson_steel,
             )
             for radius in radii
         ])
@@ -4026,6 +4110,44 @@ def plot_cs_radial_hoop_stress_profile(
         axis.set_title("Central Solenoid Radial Hoop Stress Profile")
         axis.grid(True, alpha=0.3)
 
+    def plot_cs_radial_stress_profile(
+        self,
+        axis: plt.Axes,
+        mfile: mf.MFile,
+        scan: int,
+        j_cs: float,
+        b_cs_inner: float,
+    ):
+        r_cs_inner = mfile.get("r_cs_inner", scan=scan)
+        r_cs_outer = mfile.get("r_cs_outer", scan=scan)
+
+        radii = np.linspace(r_cs_inner, r_cs_outer, num=10)
+        stress_values = np.array([
+            self.calculate_cs_radial_stress(
+                r_stress_point=radius,
+                r_cs_inner=r_cs_inner,
+                r_cs_outer=r_cs_outer,
+                j_cs=j_cs,
+                b_cs_inner=b_cs_inner,
+                f_poisson_cs_structure=tfv.poisson_steel,
+            )
+            for radius in radii
+        ])
+
+        axis.plot(
+            radii,
+            stress_values / 1e6,
+            "o-",
+            linewidth=2,
+            markersize=8,
+            label="$\\sigma_{r}$,Radial Stress",
+        )
+        axis.set_xlabel("Radial Position (m)")
+        axis.set_ylabel("Radial Stress (MPa)")
+        axis.minorticks_on()
+        axis.grid(True, alpha=0.3)
+        axis.legend(loc="best")
+
 
 def peak_b_field_at_pf_coil(
     n_coil: int, n_coil_group: int, t_b_field_peak: int, data: DataStructure

From d8d937b85a024f3ce12ad072f77aa68f619b52b2 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 15:11:01 +0100
Subject: [PATCH 16/26] Refactor CSCoil class to rename and enhance magnetic
 field calculation methods

---
 process/models/pfcoil.py | 57 ++++++++++++++++++++++++++++++++++------
 1 file changed, 49 insertions(+), 8 deletions(-)

diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 35c7a44752..ddd2d56506 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -3606,14 +3606,15 @@ def ohcalc(self):
             )
         self.calculate_cs_self_midplane_axial_stress_time_profile()
 
-    def calculate_cs_self_peak_magnetic_field(
-        self,
+    @staticmethod
+    @numba.njit(cache=True)
+    def calculate_cs_bore_magnetic_field(
         j_cs: float,
         r_cs_inner: float,
         r_cs_outer: float,
         dz_cs_half: float,
     ) -> float:
-        """Calculates the maximum field of a solenoid of circular winding and rectangular cross-section.
+        """Calculates the magnetic field at the centre of the bore of a solenoid of circular winding and rectangular cross-section.
 
         Parameters
         ----------
@@ -3629,19 +3630,17 @@ def calculate_cs_self_peak_magnetic_field(
         Returns
         -------
         float
-            Maximum field of solenoid (T)
+            Magnetic field at the centre of the bore of the solenoid (T)
 
         References
         ----------
-            - Fits are taken from the figure on p.22 of M. Wilson's book
-            "Superconducting Magnets", Clarendon Press, Oxford, N.Y., 1983,
+            -"Superconducting Magnets", Clarendon Press, Oxford, N.Y., 1983,
             ISBN 13: 9780198548102
         """
         beta = dz_cs_half / r_cs_inner
         alpha = r_cs_outer / r_cs_inner
 
-        # Field at the centre of the bore R=0, Z=0 of the solenoid
-        b_cs_bore_centre = (
+        return (
             j_cs
             * constants.RMU0
             * dz_cs_half
@@ -3651,6 +3650,48 @@ def calculate_cs_self_peak_magnetic_field(
             )
         )
 
+    def calculate_cs_self_peak_magnetic_field(
+        self,
+        j_cs: float,
+        r_cs_inner: float,
+        r_cs_outer: float,
+        dz_cs_half: float,
+    ) -> float:
+        """Calculates the maximum field of a solenoid of circular winding and rectangular cross-section.
+
+        Parameters
+        ----------
+        j_cs : float
+            Overall current density (A/m²)
+        r_cs_inner : float
+            Solenoid inner radius (m)
+        r_cs_outer : float
+            Solenoid outer radius (m)
+        dz_cs_half : float
+            Solenoid half height (m)
+
+        Returns
+        -------
+        float
+            Maximum field of solenoid (T)
+
+        References
+        ----------
+            - Fits are taken from the figure on p.22 of M. Wilson's book
+            "Superconducting Magnets", Clarendon Press, Oxford, N.Y., 1983,
+            ISBN 13: 9780198548102
+        """
+        beta = dz_cs_half / r_cs_inner
+        alpha = r_cs_outer / r_cs_inner
+
+        # Field at the centre of the bore R=0, Z=0 of the solenoid
+        b_cs_bore_centre = self.calculate_cs_bore_magnetic_field(
+            j_cs=j_cs,
+            r_cs_inner=r_cs_inner,
+            r_cs_outer=r_cs_outer,
+            dz_cs_half=dz_cs_half,
+        )
+
         # Fits are for 1 < alpha < 2 , and 0.5 < beta < very large
         if beta > 3.0:
             b1 = constants.RMU0 * j_cs * (r_cs_outer - r_cs_inner)

From 85fbe8ce36e1e9bf07b246c9c1f19007bff99b87 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 15:20:47 +0100
Subject: [PATCH 17/26] Add peak field variable for central solenoid bore in PF
 coil module

---
 process/data_structure/pfcoil_variables.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
index 6713426f59..3167f12b68 100644
--- a/process/data_structure/pfcoil_variables.py
+++ b/process/data_structure/pfcoil_variables.py
@@ -122,6 +122,8 @@ class PFCoilData:
     b_cs_peak_pulse_start: float = 0.0
     """maximum field in central solenoid at beginning of pulse (T)"""
 
+b_cs_bore_peak: float = None
+"""peak field in central solenoid bore (T)"""
     b_pf_coil_peak: list[float] = field(default_factory=lambda: np.zeros(NGC2))
     """peak field at coil i (T)"""
 

From 5302620c32e9c8c782a687ed3268cd1839007a04 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 15:34:04 +0100
Subject: [PATCH 18/26] Add magnetic field variable for outer edge of central
 solenoid at midplane

---
 process/data_structure/pfcoil_variables.py | 2 ++
 process/models/pfcoil.py                   | 3 +++
 2 files changed, 5 insertions(+)

diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
index 3167f12b68..030834c17a 100644
--- a/process/data_structure/pfcoil_variables.py
+++ b/process/data_structure/pfcoil_variables.py
@@ -116,6 +116,8 @@ class PFCoilData:
     a_cs_cable_space: float = 0.0
     """central solenoid conductor+void area with area of steel subtracted (m2)"""
 
+b_cs_self_outer_midplane: float = None
+"""magnetic field at outer edge of central solenoid at midplane due to its own current (T)"""
     b_cs_peak_flat_top_end: float = 0.0
     """maximum field in central solenoid at end of flat-top (EoF) (T)"""
 
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index ddd2d56506..0d5df4a64c 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -3318,6 +3318,9 @@ def ohcalc(self):
 
         # Peak field on outboard side of central Solenoid
         # (self-field is assumed to be zero - long solenoid approximation)
+
+        pfcoil_variables.b_cs_self_outer_midplane = 0.0
+
         bohco = abs(bzo)
 
         # Peak field at the Beginning-Of-Pulse (BOP)

From 031be2d637da8e067aa33541f553858e0b8d48eb Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 17:08:17 +0100
Subject: [PATCH 19/26] Add operating temperature variable for central solenoid
 and update initialization

---
 process/data_structure/pfcoil_variables.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
index 030834c17a..3704179a0a 100644
--- a/process/data_structure/pfcoil_variables.py
+++ b/process/data_structure/pfcoil_variables.py
@@ -412,6 +412,9 @@ class PFCoilData:
     temp_cs_superconductor_margin: float = 0.0
     """Central solenoid temperature margin (K)"""
 
+temp_cs_superconductor_operating: float = None
+"""Central solenoid operating temperature (K)"""
+
     n_pf_coil_turns: list[float] = field(default_factory=lambda: np.zeros(NGC2))
     """number of turns in PF coil i"""
 

From 5661812cf9e01e027c8d109fd3e7626a53055ef4 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 7 Apr 2026 17:33:57 +0100
Subject: [PATCH 20/26] Update operating temperature for central solenoid
 superconductor

---
 process/models/pfcoil.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 0d5df4a64c..4a550d1612 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -3523,7 +3523,7 @@ def ohcalc(self):
                 self.data.pf_coil.i_cs_superconductor,
                 tfv.fhts,
                 tfv.str_cs_con_res,
-                tfv.tftmp,
+                pfcoil_variables.temp_cs_superconductor_operating,
                 tfv.bcritsc,
                 tfv.tcritsc,
             )
@@ -3568,7 +3568,7 @@ def ohcalc(self):
                 self.data.pf_coil.i_cs_superconductor,
                 tfv.fhts,
                 tfv.str_cs_con_res,
-                tfv.tftmp,
+                pfcoil_variables.temp_cs_superconductor_operating,
                 tfv.bcritsc,
                 tfv.tcritsc,
             )

From fc93fd808abea5c0cb8bd8fd225a4b1c690a5525 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 19 May 2026 15:39:08 +0100
Subject: [PATCH 21/26] Post rebase fixes

---
 process/core/input.py                      |   3 +
 process/core/io/plot/summary.py            |   5 +-
 process/data_structure/pfcoil_variables.py |  34 ++--
 process/models/costs/costs.py              |  34 ++--
 process/models/pfcoil.py                   | 187 ++++++++++-----------
 tests/unit/models/test_costs_1990.py       |   2 +-
 tests/unit/models/test_pfcoil.py           |   2 +-
 7 files changed, 133 insertions(+), 134 deletions(-)

diff --git a/process/core/input.py b/process/core/input.py
index 565a130876..391f7e1916 100644
--- a/process/core/input.py
+++ b/process/core/input.py
@@ -972,6 +972,9 @@ def __post_init__(self):
     ),
     "f_dr_tf_outboard_inboard": InputVariable("build", float, range=(0.2, 5.0)),
     "tftmp": InputVariable(data_structure.tfcoil_variables, float, range=(0.01, 293.0)),
+    "temp_cs_superconductor_operating": InputVariable(
+        data_structure.pfcoil_variables, float, range=(0.01, 293.0)
+    ),
     "tgain": InputVariable("ife", float, range=(1.0, 500.0)),
     "th_joint_contact": InputVariable(
         data_structure.tfcoil_variables, float, range=(0.0, 1.0)
diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 922b0c9686..780809c265 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -48,6 +48,7 @@
     vacuum_vessel_geometry_double_null,
     vacuum_vessel_geometry_single_null,
 )
+from process.models.pfcoil import CSCoil
 from process.models.physics.bootstrap_current import BootstrapCurrentFractionModel
 from process.models.physics.confinement_time import (
     ConfinementTimeModel,
@@ -14497,9 +14498,7 @@ def main_plot(
 
     plot_current_profiles_over_time(figs[29].add_subplot(111), m_file, scan)
 
-    CSCoil.plot_stress_time_profile(
-        axis=figs[29].add_subplot(322), mfile=m_file, scan=scan
-    )
+    plot_cs_stress_time_profile(axis=figs[29].add_subplot(322), mfile=m_file, scan=scan)
 
     cs_coil = CSCoil(cs_fatigue=CsFatigue())
     cs_coil.plot_cs_radial_hoop_stress_profile(
diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
index 3704179a0a..3599b85f7e 100644
--- a/process/data_structure/pfcoil_variables.py
+++ b/process/data_structure/pfcoil_variables.py
@@ -43,10 +43,12 @@ class PFCoilData:
     """Axial stress (z) in central solenoid at midplane due to its own field at each time point (Pa)"""
 
     stress_hoop_cs_inner: float = 0.0
-"""Hoop stress in central solenoid at inboard edge due to its own field (when at peak current) (Pa)"""
+    """Hoop stress in central solenoid at inboard edge due to its own field (when at peak current) (Pa)"""
 
-stress_hoop_cs_inner_profile: list[float] = None
-"""Hoop stress in central solenoid at inboard edge due to its own field at each time point (Pa)"""
+    stress_hoop_cs_inner_profile: list[float] = field(
+        default_factory=lambda: np.zeros(6)
+    )
+    """Hoop stress in central solenoid at inboard edge due to its own field at each time point (Pa)"""
 
     forc_z_cs_self_peak_midplane: float = 0.0
     """Axial force (z) on central solenoid at midplane due to its own field (when at peak current) (N)"""
@@ -107,8 +109,8 @@ class PFCoilData:
     a_cs_poloidal: float = 0.0
     """Central solenoid vertical cross-sectional area (m2)"""
 
-a_cs_steel_poloidal: float = None
-"""Central solenoid vertical cross-sectional area of steel (m2)"""
+    a_cs_steel_poloidal: float = 0.0
+    """Central solenoid vertical cross-sectional area of steel (m2)"""
 
     a_cs_turn: float = 0.0
     """Central solenoid (OH) trun cross-sectional area (m2)"""
@@ -116,16 +118,18 @@ class PFCoilData:
     a_cs_cable_space: float = 0.0
     """central solenoid conductor+void area with area of steel subtracted (m2)"""
 
-b_cs_self_outer_midplane: float = None
-"""magnetic field at outer edge of central solenoid at midplane due to its own current (T)"""
+    b_cs_self_outer_midplane: float = 0.0
+    """magnetic field at outer edge of central solenoid at midplane due to its own current (T)"""
+
     b_cs_peak_flat_top_end: float = 0.0
     """maximum field in central solenoid at end of flat-top (EoF) (T)"""
 
     b_cs_peak_pulse_start: float = 0.0
     """maximum field in central solenoid at beginning of pulse (T)"""
 
-b_cs_bore_peak: float = None
-"""peak field in central solenoid bore (T)"""
+    b_cs_bore_peak: float = 0.0
+    """peak field in central solenoid bore (T)"""
+
     b_pf_coil_peak: list[float] = field(default_factory=lambda: np.zeros(NGC2))
     """peak field at coil i (T)"""
 
@@ -358,11 +362,11 @@ class PFCoilData:
     r_cs_middle: float = 0.0
     """radius to the centre of the central solenoid (m)"""
 
-r_cs_inner: float = None
-"""inner radius of the central solenoid (m)"""
+    r_cs_inner: float = 0.0
+    """inner radius of the central solenoid (m)"""
 
-r_cs_outer: float = None
-"""outer radius of the central solenoid (m)"""
+    r_cs_outer: float = 0.0
+    """outer radius of the central solenoid (m)"""
 
     dz_cs_full: float = 0.0
     """Full height of the central solenoid (m)"""
@@ -412,8 +416,8 @@ class PFCoilData:
     temp_cs_superconductor_margin: float = 0.0
     """Central solenoid temperature margin (K)"""
 
-temp_cs_superconductor_operating: float = None
-"""Central solenoid operating temperature (K)"""
+    temp_cs_superconductor_operating: float = 4.75
+    """Central solenoid operating temperature (K)"""
 
     n_pf_coil_turns: list[float] = field(default_factory=lambda: np.zeros(NGC2))
     """number of turns in PF coil i"""
diff --git a/process/models/costs/costs.py b/process/models/costs/costs.py
index f399887d3f..2b08ff3071 100644
--- a/process/models/costs/costs.py
+++ b/process/models/costs/costs.py
@@ -1715,12 +1715,12 @@ def acc2222(self):
                 #  Issue #328  Use CS conductor cross-sectional area (m2)
                 if self.data.pf_coil.i_pf_conductor == 0:
                     costpfsc = (
-                        cost_variables.ucsc[pfcoil_variables.i_cs_superconductor - 1]
-                        * pfcoil_variables.a_cs_cable_space
-                        * (1 - pfcoil_variables.f_a_cs_void)
-                        * (1 - pfcoil_variables.fcuohsu)
-                        / pfcoil_variables.n_pf_coil_turns[
-                            pfcoil_variables.n_cs_pf_coils - 1
+                        self.data.costs.ucsc[self.data.pf_coil.i_cs_superconductor - 1]
+                        * self.data.pf_coil.a_cs_cable_space
+                        * (1 - self.data.pf_coil.f_a_cs_void)
+                        * (1 - self.data.pf_coil.fcuohsu)
+                        / self.data.pf_coil.n_pf_coil_turns[
+                            self.data.pf_coil.n_cs_pf_coils - 1
                         ]
                         * tfcoil_variables.dcond[
                             self.data.pf_coil.i_cs_superconductor - 1
@@ -1745,23 +1745,23 @@ def acc2222(self):
 
             if self.data.pf_coil.i_pf_conductor == 0:
                 costpfcu = (
-                    cost_variables.uccu
-                    * pfcoil_variables.a_cs_cable_space
-                    * (1 - pfcoil_variables.f_a_cs_void)
-                    * pfcoil_variables.fcuohsu
-                    / pfcoil_variables.n_pf_coil_turns[
-                        pfcoil_variables.n_cs_pf_coils - 1
+                    self.data.costs.uccu
+                    * self.data.pf_coil.a_cs_cable_space
+                    * (1 - self.data.pf_coil.f_a_cs_void)
+                    * self.data.pf_coil.fcuohsu
+                    / self.data.pf_coil.n_pf_coil_turns[
+                        self.data.pf_coil.n_cs_pf_coils - 1
                     ]
                     * constants.den_copper
                 )
             else:
                 # MDK I don't know if this is ccorrect as we never use the resistive model
                 costpfcu = (
-                    cost_variables.uccu
-                    * pfcoil_variables.a_cs_cable_space
-                    * (1 - pfcoil_variables.f_a_cs_void)
-                    / pfcoil_variables.n_pf_coil_turns[
-                        pfcoil_variables.n_cs_pf_coils - 1
+                    self.data.costs.uccu
+                    * self.data.pf_coil.a_cs_cable_space
+                    * (1 - self.data.pf_coil.f_a_cs_void)
+                    / self.data.pf_coil.n_pf_coil_turns[
+                        self.data.pf_coil.n_cs_pf_coils - 1
                     ]
                     * constants.den_copper
                 )
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 4a550d1612..0df1d0eecc 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -1,6 +1,7 @@
 import logging
 import math
 
+import matplotlib.pyplot as plt
 import numba
 import numpy as np
 from scipy import optimize
@@ -10,6 +11,7 @@
 from process.core import constants
 from process.core import process_output as op
 from process.core.exceptions import ProcessValueError
+from process.core.io.mfile import MFile
 from process.core.model import DataStructure, Model
 from process.data_structure import physics_variables as pv
 from process.data_structure import rebco_variables as rcv
@@ -2145,47 +2147,47 @@ def outpf(self):
                 self.data.pf_coil.r_cs_middle,
                 "OP ",
             )
-                op.ovarre(
-                    self.outfile,
-                    "CS radial inner (m)",
-                    "(r_cs_inner)",
-                    pfcoil_variables.r_cs_inner,
-                    "OP ",
-                )
-                op.ovarre(
-                    self.outfile,
-                    "CS radial outer (m)",
-                    "(r_cs_outer)",
-                    pfcoil_variables.r_cs_outer,
-                    "OP ",
-                )
+            op.ovarre(
+                self.outfile,
+                "CS radial inner (m)",
+                "(r_cs_inner)",
+                self.data.pf_coil.r_cs_inner,
+                "OP ",
+            )
+            op.ovarre(
+                self.outfile,
+                "CS radial outer (m)",
+                "(r_cs_outer)",
+                self.data.pf_coil.r_cs_outer,
+                "OP ",
+            )
             op.ovarre(
                 self.outfile,
                 "CS conductor+void cross-sectional area (m2)",
                 "(a_cs_cable_space)",
-                pfcoil_variables.a_cs_cable_space,
+                self.data.pf_coil.a_cs_cable_space,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
                 "   CS conductor cross-sectional area (m2)",
                 "(a_cs_cable_space*(1-f_a_cs_void))",
-                pfcoil_variables.a_cs_cable_space
-                    * (1.0e0 - pfcoil_variables.f_a_cs_void),
+                self.data.pf_coil.a_cs_cable_space
+                * (1.0e0 - self.data.pf_coil.f_a_cs_void),
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
                 "   CS void cross-sectional area (m2)",
                 "(a_cs_cable_space*f_a_cs_void)",
-                pfcoil_variables.a_cs_cable_space * pfcoil_variables.f_a_cs_void,
+                self.data.pf_coil.a_cs_cable_space * self.data.pf_coil.f_a_cs_void,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
                 "CS steel cross-sectional area (m2)",
                 "(a_cs_steel_poloidal)",
-                pfcoil_variables.a_cs_steel_poloidal,
+                self.data.pf_coil.a_cs_steel_poloidal,
                 "OP ",
             )
             op.ovarre(
@@ -2215,7 +2217,7 @@ def outpf(self):
                 self.outfile,
                 "Hoop stress in CS steel (Pa)",
                 "(stress_hoop_cs_inner)",
-                pfcoil_variables.stress_hoop_cs_inner,
+                self.data.pf_coil.stress_hoop_cs_inner,
                 "OP ",
             )
             op.ovarre(
@@ -3222,11 +3224,11 @@ def ohcalc(self):
             dr_bore=self.data.build.dr_bore,
         )
 
-        pfcoil_variables.r_cs_inner = pfcoil_variables.r_pf_coil_inner[
-            pfcoil_variables.n_cs_pf_coils - 1
+        self.data.pf_coil.r_cs_inner = self.data.pf_coil.r_pf_coil_inner[
+            self.data.pf_coil.n_cs_pf_coils - 1
         ]
-        pfcoil_variables.r_cs_outer = pfcoil_variables.r_pf_coil_outer[
-            pfcoil_variables.n_cs_pf_coils - 1
+        self.data.pf_coil.r_cs_outer = self.data.pf_coil.r_pf_coil_outer[
+            self.data.pf_coil.n_cs_pf_coils - 1
         ]
 
         # Maximum current (MA-turns) in central Solenoid, at either BOP or EOF
@@ -3293,9 +3295,9 @@ def ohcalc(self):
 
         # Peak field due to central Solenoid itself
         b_cs_self_peak_flat_top_end = self.calculate_cs_self_peak_magnetic_field(
-            j_cs=pfcoil_variables.j_cs_flat_top_end,
-            r_cs_inner=pfcoil_variables.r_pf_coil_inner[
-                pfcoil_variables.n_cs_pf_coils - 1
+            j_cs=self.data.pf_coil.j_cs_flat_top_end,
+            r_cs_inner=self.data.pf_coil.r_pf_coil_inner[
+                self.data.pf_coil.n_cs_pf_coils - 1
             ],
             r_cs_outer=self.data.pf_coil.r_pf_coil_outer[
                 self.data.pf_coil.n_cs_pf_coils - 1
@@ -3308,18 +3310,18 @@ def ohcalc(self):
         # Peak field due to other PF coils plus plasma
         timepoint = 5
         _, _, bzi, bzo = peak_b_field_at_pf_coil(
-            n_coil=pfcoil_variables.n_cs_pf_coils,
+            n_coil=self.data.pf_coil.n_cs_pf_coils,
             n_coil_group=99,
             t_b_field_peak=timepoint,
             data=self.data,
         )
 
-        pfcoil_variables.b_cs_peak_flat_top_end = abs(bzi - b_cs_self_peak_flat_top_end)
+        self.data.pf_coil.b_cs_peak_flat_top_end = abs(bzi - b_cs_self_peak_flat_top_end)
 
         # Peak field on outboard side of central Solenoid
         # (self-field is assumed to be zero - long solenoid approximation)
 
-        pfcoil_variables.b_cs_self_outer_midplane = 0.0
+        self.data.pf_coil.b_cs_self_outer_midplane = 0.0
 
         bohco = abs(bzo)
 
@@ -3341,7 +3343,7 @@ def ohcalc(self):
         )
         timepoint = 2
         _, _, bzi, bzo = peak_b_field_at_pf_coil(
-            n_coil=pfcoil_variables.n_cs_pf_coils,
+            n_coil=self.data.pf_coil.n_cs_pf_coils,
             n_coil_group=99,
             t_b_field_peak=timepoint,
             data=self.data,
@@ -3365,18 +3367,18 @@ def ohcalc(self):
             # Superconducting coil
 
             # New calculation from M. N. Wilson for hoop stress
-            pfcoil_variables.stress_hoop_cs_inner = self.calculate_cs_hoop_stress(
-                r_stress_point=pfcoil_variables.r_pf_coil_inner[
-                    pfcoil_variables.n_cs_pf_coils - 1
+            self.data.pf_coil.stress_hoop_cs_inner = self.calculate_cs_hoop_stress(
+                r_stress_point=self.data.pf_coil.r_pf_coil_inner[
+                    self.data.pf_coil.n_cs_pf_coils - 1
                 ],
-                r_cs_inner=pfcoil_variables.r_pf_coil_inner[
-                    pfcoil_variables.n_cs_pf_coils - 1
+                r_cs_inner=self.data.pf_coil.r_pf_coil_inner[
+                    self.data.pf_coil.n_cs_pf_coils - 1
                 ],
-                r_cs_outer=pfcoil_variables.r_pf_coil_outer[
-                    pfcoil_variables.n_cs_pf_coils - 1
+                r_cs_outer=self.data.pf_coil.r_pf_coil_outer[
+                    self.data.pf_coil.n_cs_pf_coils - 1
                 ],
-                j_cs=pfcoil_variables.j_cs_pulse_start,
-                b_cs_inner=pfcoil_variables.b_cs_peak_pulse_start,
+                j_cs=self.data.pf_coil.j_cs_pulse_start,
+                b_cs_inner=self.data.pf_coil.b_cs_peak_pulse_start,
                 f_poisson_cs_structure=tfv.poisson_steel,
             )
 
@@ -3409,7 +3411,7 @@ def ohcalc(self):
                     self.data.cs_fatigue.n_cycle,
                     self.data.cs_fatigue.t_crack_radial,
                 ) = self.cs_fatigue.ncycle(
-                    pfcoil_variables.stress_hoop_cs_inner,
+                    self.data.pf_coil.stress_hoop_cs_inner,
                     self.data.cs_fatigue.residual_sig_hoop,
                     self.data.cs_fatigue.t_crack_vertical,
                     self.data.cs_fatigue.dz_cs_turn_conduit,
@@ -3420,24 +3422,24 @@ def ohcalc(self):
             # equation is used for Central Solenoid stress
 
             # Area of steel in Central Solenoid
-            pfcoil_variables.a_cs_steel_poloidal = (
-                pfcoil_variables.f_a_cs_turn_steel * pfcoil_variables.a_cs_poloidal
+            self.data.pf_coil.a_cs_steel_poloidal = (
+                self.data.pf_coil.f_a_cs_turn_steel * self.data.pf_coil.a_cs_poloidal
             )
 
             if self.data.pf_coil.i_cs_stress == 1:
                 self.data.pf_coil.s_shear_cs_peak = max(
                     abs(
-                        pfcoil_variables.stress_hoop_cs_inner
-                        - pfcoil_variables.stress_z_cs_self_peak_midplane
+                        self.data.pf_coil.stress_hoop_cs_inner
+                        - self.data.pf_coil.stress_z_cs_self_peak_midplane
                     ),
-                    abs(pfcoil_variables.stress_z_cs_self_peak_midplane - 0.0e0),
-                    abs(0.0e0 - pfcoil_variables.stress_hoop_cs_inner),
+                    abs(self.data.pf_coil.stress_z_cs_self_peak_midplane - 0.0e0),
+                    abs(0.0e0 - self.data.pf_coil.stress_hoop_cs_inner),
                 )
             else:
-                pfcoil_variables.s_shear_cs_peak = max(
-                    abs(pfcoil_variables.stress_hoop_cs_inner - 0.0e0),
+                self.data.pf_coil.s_shear_cs_peak = max(
+                    abs(self.data.pf_coil.stress_hoop_cs_inner - 0.0e0),
                     abs(0.0e0 - 0.0e0),
-                    abs(0.0e0 - pfcoil_variables.stress_hoop_cs_inner),
+                    abs(0.0e0 - self.data.pf_coil.stress_hoop_cs_inner),
                 )
 
             # Thickness of hypothetical steel cylinders assumed to encase the CS along
@@ -3445,19 +3447,19 @@ def ohcalc(self):
             # throughout the conductor
             self.data.pf_coil.pfcaseth[self.data.pf_coil.n_cs_pf_coils - 1] = (
                 0.25e0
-                * pfcoil_variables.a_cs_steel_poloidal
-                / pfcoil_variables.z_pf_coil_upper[pfcoil_variables.n_cs_pf_coils - 1]
+                * self.data.pf_coil.a_cs_steel_poloidal
+                / self.data.pf_coil.z_pf_coil_upper[self.data.pf_coil.n_cs_pf_coils - 1]
             )
 
         else:
-            pfcoil_variables.a_cs_steel_poloidal = (
+            self.data.pf_coil.a_cs_steel_poloidal = (
                 0.0e0  # Resistive Central Solenoid - no steel needed
             )
-            pfcoil_variables.pfcaseth[pfcoil_variables.n_cs_pf_coils - 1] = 0.0e0
+            self.data.pf_coil.pfcaseth[self.data.pf_coil.n_cs_pf_coils - 1] = 0.0e0
 
         # Weight of steel
-        pfcoil_variables.m_pf_coil_structure[pfcoil_variables.n_cs_pf_coils - 1] = (
-            pfcoil_variables.a_cs_steel_poloidal
+        self.data.pf_coil.m_pf_coil_structure[self.data.pf_coil.n_cs_pf_coils - 1] = (
+            self.data.pf_coil.a_cs_steel_poloidal
             * 2.0e0
             * np.pi
             * self.data.pf_coil.r_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1]
@@ -3465,33 +3467,37 @@ def ohcalc(self):
         )
 
         # Non-steel cross-sectional area
-        pfcoil_variables.a_cs_cable_space = (
-            pfcoil_variables.a_cs_poloidal - pfcoil_variables.a_cs_steel_poloidal
+        self.data.pf_coil.a_cs_cable_space = (
+            self.data.pf_coil.a_cs_poloidal - self.data.pf_coil.a_cs_steel_poloidal
         )
 
         # Issue #97. Fudge to ensure a_cs_cable_space is positive; result is continuous, smooth and
         # monotonically decreases
 
         da = 0.0001e0  # 1 cm^2
-        if pfcoil_variables.a_cs_cable_space < da:
-            pfcoil_variables.a_cs_cable_space = (
-                da * da / (2.0e0 * da - pfcoil_variables.a_cs_cable_space)
+        if self.data.pf_coil.a_cs_cable_space < da:
+            self.data.pf_coil.a_cs_cable_space = (
+                da * da / (2.0e0 * da - self.data.pf_coil.a_cs_cable_space)
             )
 
         # Weight of conductor in central Solenoid
-        if pfcoil_variables.i_pf_conductor == 0:
-            pfcoil_variables.m_pf_coil_conductor[pfcoil_variables.n_cs_pf_coils - 1] = (
-                pfcoil_variables.a_cs_cable_space
-                * (1.0e0 - pfcoil_variables.f_a_cs_void)
+        if self.data.pf_coil.i_pf_conductor == 0:
+            self.data.pf_coil.m_pf_coil_conductor[
+                self.data.pf_coil.n_cs_pf_coils - 1
+            ] = (
+                self.data.pf_coil.a_cs_cable_space
+                * (1.0e0 - self.data.pf_coil.f_a_cs_void)
                 * 2.0e0
                 * np.pi
                 * self.data.pf_coil.r_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1]
                 * tfv.dcond[self.data.pf_coil.i_cs_superconductor - 1]
             )
         else:
-            pfcoil_variables.m_pf_coil_conductor[pfcoil_variables.n_cs_pf_coils - 1] = (
-                pfcoil_variables.a_cs_cable_space
-                * (1.0e0 - pfcoil_variables.f_a_cs_void)
+            self.data.pf_coil.m_pf_coil_conductor[
+                self.data.pf_coil.n_cs_pf_coils - 1
+            ] = (
+                self.data.pf_coil.a_cs_cable_space
+                * (1.0e0 - self.data.pf_coil.f_a_cs_void)
                 * 2.0e0
                 * np.pi
                 * self.data.pf_coil.r_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1]
@@ -3517,13 +3523,13 @@ def ohcalc(self):
                             self.data.pf_coil.n_cs_pf_coils - 1
                         ]
                     )
-                    / pfcoil_variables.a_cs_cable_space
+                    / self.data.pf_coil.a_cs_cable_space
                 )
                 * 1.0e6,
                 self.data.pf_coil.i_cs_superconductor,
                 tfv.fhts,
                 tfv.str_cs_con_res,
-                pfcoil_variables.temp_cs_superconductor_operating,
+                self.data.pf_coil.temp_cs_superconductor_operating,
                 tfv.bcritsc,
                 tfv.tcritsc,
             )
@@ -3539,10 +3545,10 @@ def ohcalc(self):
                     * (1 - self.data.pf_coil.fcuohsu)
                 )
 
-            pfcoil_variables.j_cs_critical_flat_top_end = (
+            self.data.pf_coil.j_cs_critical_flat_top_end = (
                 jcritwp
-                * pfcoil_variables.a_cs_cable_space
-                / pfcoil_variables.a_cs_poloidal
+                * self.data.pf_coil.a_cs_cable_space
+                / self.data.pf_coil.a_cs_poloidal
             )
 
             # Allowable coil overall current density at BOP
@@ -3562,21 +3568,21 @@ def ohcalc(self):
                             self.data.pf_coil.n_cs_pf_coils - 1
                         ]
                     )
-                    / pfcoil_variables.a_cs_cable_space
+                    / self.data.pf_coil.a_cs_cable_space
                 )
                 * 1.0e6,
                 self.data.pf_coil.i_cs_superconductor,
                 tfv.fhts,
                 tfv.str_cs_con_res,
-                pfcoil_variables.temp_cs_superconductor_operating,
+                self.data.pf_coil.temp_cs_superconductor_operating,
                 tfv.bcritsc,
                 tfv.tcritsc,
             )
 
-            pfcoil_variables.j_pf_wp_critical[pfcoil_variables.n_cs_pf_coils - 1] = (
+            self.data.pf_coil.j_pf_wp_critical[self.data.pf_coil.n_cs_pf_coils - 1] = (
                 jcritwp
-                * pfcoil_variables.a_cs_cable_space
-                / pfcoil_variables.a_cs_poloidal
+                * self.data.pf_coil.a_cs_cable_space
+                / self.data.pf_coil.a_cs_poloidal
             )
             self.data.pf_coil.j_cs_critical_pulse_start = (
                 self.data.pf_coil.j_pf_wp_critical[self.data.pf_coil.n_cs_pf_coils - 1]
@@ -3610,7 +3616,6 @@ def ohcalc(self):
         self.calculate_cs_self_midplane_axial_stress_time_profile()
 
     @staticmethod
-    @numba.njit(cache=True)
     def calculate_cs_bore_magnetic_field(
         j_cs: float,
         r_cs_inner: float,
@@ -3896,9 +3901,8 @@ def calculate_cs_self_midplane_axial_stress_time_profile(
             )
             self.data.pf_coil.stress_z_cs_self_midplane_profile[time] = stress_value
 
-    @staticmethod
-    @numba.njit(cache=True)
     def calculate_cs_hoop_stress(
+        self,
         r_stress_point: float | np.ndarray,
         r_cs_inner: float,
         r_cs_outer: float,
@@ -3936,7 +3940,6 @@ def calculate_cs_hoop_stress(
         - M. N. Wilson, Superconducting Magnets. Oxford University Press, USA, 1983.
         ‌
         """
-
         alpha = r_cs_outer / r_cs_inner
 
         # alpha
@@ -3949,9 +3952,6 @@ def calculate_cs_hoop_stress(
         # Assume to be 0 for now
         b_b = 0.0e0
 
-        # current density [A/m^2]
-        j_cs = pfcoil_variables.j_cs_pulse_start
-
         # K term
         k = ((alpha * b_cs_inner - b_b) * j_cs * r_cs_inner) / (alpha - 1.0e0)
 
@@ -3988,11 +3988,10 @@ def calculate_cs_hoop_stress(
 
         s_hoop_nom = hp_term_1 * hp_term_2 - hp_term_3 * hp_term_4
 
-        return s_hoop_nom / pfcoil_variables.f_a_cs_turn_steel
+        return s_hoop_nom / self.data.pf_coil.f_a_cs_turn_steel
 
-    @staticmethod
-    @numba.njit(cache=True)
     def calculate_cs_radial_stress(
+        self,
         r_stress_point: float | np.ndarray,
         r_cs_inner: float,
         r_cs_outer: float,
@@ -4030,9 +4029,6 @@ def calculate_cs_radial_stress(
         - M. N. Wilson, Superconducting Magnets. Oxford University Press, USA, 1983.
         ‌
         """
-
-        alpha = r_cs_outer / r_cs_inner
-
         # alpha
         alpha = r_cs_outer / r_cs_inner
 
@@ -4043,9 +4039,6 @@ def calculate_cs_radial_stress(
         # Assume to be 0 for now
         b_b = 0.0e0
 
-        # current density [A/m^2]
-        j_cs = pfcoil_variables.j_cs_pulse_start
-
         # K term
         k = ((alpha * b_cs_inner - b_b) * j_cs * r_cs_inner) / (alpha - 1.0e0)
 
@@ -4066,7 +4059,7 @@ def calculate_cs_radial_stress(
         return hp_term_1 * hp_term_2 - hp_term_3 * hp_term_4
 
     @staticmethod
-    def plot_stress_time_profile(axis: plt.Axes, mfile: mf.MFile, scan: int):
+    def plot_stress_time_profile(axis: plt.Axes, mfile: MFile, scan: int):
         t_plant_pulse_coil_precharge = mfile.get(
             "t_plant_pulse_coil_precharge", scan=scan
         )
@@ -4118,7 +4111,7 @@ def plot_stress_time_profile(axis: plt.Axes, mfile: mf.MFile, scan: int):
     def plot_cs_radial_hoop_stress_profile(
         self,
         axis: plt.Axes,
-        mfile: mf.MFile,
+        mfile: MFile,
         scan: int,
         j_cs: float,
         b_cs_inner: float,
@@ -4157,7 +4150,7 @@ def plot_cs_radial_hoop_stress_profile(
     def plot_cs_radial_stress_profile(
         self,
         axis: plt.Axes,
-        mfile: mf.MFile,
+        mfile: MFile,
         scan: int,
         j_cs: float,
         b_cs_inner: float,
diff --git a/tests/unit/models/test_costs_1990.py b/tests/unit/models/test_costs_1990.py
index e11834ce01..cbbcd339ae 100644
--- a/tests/unit/models/test_costs_1990.py
+++ b/tests/unit/models/test_costs_1990.py
@@ -2768,7 +2768,7 @@ def test_acc2222(acc2222param, monkeypatch, costs):
     )
 
     monkeypatch.setattr(
-        pfcoil_variables, "a_cs_cable_space", acc2222param.a_cs_cable_space
+        costs.data.pf_coil, "a_cs_cable_space", acc2222param.a_cs_cable_space
     )
 
     monkeypatch.setattr(costs.data.structure, "fncmass", acc2222param.fncmass)
diff --git a/tests/unit/models/test_pfcoil.py b/tests/unit/models/test_pfcoil.py
index 47dd35bbbd..44f982e480 100644
--- a/tests/unit/models/test_pfcoil.py
+++ b/tests/unit/models/test_pfcoil.py
@@ -2611,7 +2611,7 @@ def test_ohcalc(monkeypatch, reinitialise_error_module, cs_coil):
     monkeypatch.setattr(cs_coil.data.pf_coil, "p_cs_resistive_flat_top", 0.0)
     monkeypatch.setattr(cs_coil.data.pf_coil, "j_cs_critical_pulse_start", 3.048e7)
     monkeypatch.setattr(cs_coil.data.pf_coil, "s_shear_cs_peak", 5.718e8)
-    monkeypatch.setattr(cs_coil.data.pf_coil, "awpoh", 4.232)
+    monkeypatch.setattr(cs_coil.data.pf_coil, "a_cs_cable_space", 4.232)
     monkeypatch.setattr(cs_coil.data.pf_coil, "f_a_cs_turn_steel", 5.926e-1)
     monkeypatch.setattr(cs_coil.data.pf_coil, "b_cs_peak_pulse_start", 1.4e1)
     monkeypatch.setattr(cs_coil.data.pf_coil, "j_cs_critical_flat_top_end", 4.070e7)

From 49c893bf5654f414cc38b8251fabf49014e078de Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 20 May 2026 08:57:11 +0100
Subject: [PATCH 22/26] Update unit test

---
 process/models/pfcoil.py         |  8 ++-
 tests/unit/models/test_pfcoil.py | 86 +++++++-------------------------
 2 files changed, 26 insertions(+), 68 deletions(-)

diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 0df1d0eecc..8cf62c9b4c 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -3380,6 +3380,7 @@ def ohcalc(self):
                 j_cs=self.data.pf_coil.j_cs_pulse_start,
                 b_cs_inner=self.data.pf_coil.b_cs_peak_pulse_start,
                 f_poisson_cs_structure=tfv.poisson_steel,
+                f_a_cs_turn_steel=self.data.pf_coil.f_a_cs_turn_steel,
             )
 
             # New calculation from Y. Iwasa for axial stress
@@ -3909,6 +3910,7 @@ def calculate_cs_hoop_stress(
         j_cs: float,
         b_cs_inner: float,
         f_poisson_cs_structure: float,
+        f_a_cs_turn_steel: float,
     ) -> float:
         """Calculation of hoop stress of central solenoid.
 
@@ -3929,6 +3931,9 @@ def calculate_cs_hoop_stress(
             Magnetic field at the inner radius of the central solenoid (T)
         f_poisson_cs_structure : float
             Poisson's ratio of the central solenoid structure (dimensionless)
+        f_a_cs_turn_steel : float
+            Steel area fraction of the central solenoid turn cross-section
+            (dimensionless)
 
         Returns
         -------
@@ -3988,7 +3993,7 @@ def calculate_cs_hoop_stress(
 
         s_hoop_nom = hp_term_1 * hp_term_2 - hp_term_3 * hp_term_4
 
-        return s_hoop_nom / self.data.pf_coil.f_a_cs_turn_steel
+        return s_hoop_nom / f_a_cs_turn_steel
 
     def calculate_cs_radial_stress(
         self,
@@ -4128,6 +4133,7 @@ def plot_cs_radial_hoop_stress_profile(
                 j_cs=j_cs,
                 b_cs_inner=b_cs_inner,
                 f_poisson_cs_structure=tfv.poisson_steel,
+                f_a_cs_turn_steel=mfile.get("f_a_cs_turn_steel", scan=scan),
             )
             for radius in radii
         ])
diff --git a/tests/unit/models/test_pfcoil.py b/tests/unit/models/test_pfcoil.py
index 44f982e480..07f268d585 100644
--- a/tests/unit/models/test_pfcoil.py
+++ b/tests/unit/models/test_pfcoil.py
@@ -1974,78 +1974,30 @@ def test_hoop_stress(cs_coil, monkeypatch):
     discovered using gdb to break on the first subroutine call when running the
     baseline 2018 IN.DAT.
 
-    :param pfcoil: PFCoil object
-    :type pfcoil: process.pfcoil.PFCoil
+    :param cs_coil: CSCoil object
+    :type cs_coil: process.pfcoil.CSCoil
     :param monkeypatch: mocking fixture
     :type monkeypatch: _pytest.monkeypatch.MonkeyPatch
     """
-    monkeypatch.setattr(cs_coil.data.pf_coil, "f_a_cs_turn_steel", 0.57874999999999999)
-    monkeypatch.setattr(
-        cs_coil.data.pf_coil, "b_cs_peak_pulse_start", 13.522197474024983
-    )
-    monkeypatch.setattr(cs_coil.data.pf_coil, "j_cs_pulse_start", 19311657.760000002)
-    monkeypatch.setattr(cs_coil.data.pf_coil, "n_cs_pf_coils", 7)
-    monkeypatch.setattr(
-        cs_coil.data.pf_coil,
-        "r_pf_coil_outer",
-        np.array([
-            6.8520884119768697,
-            6.9480065348448967,
-            18.98258241468535,
-            18.98258241468535,
-            17.22166645654087,
-            17.22166645654087,
-            2.88462,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-        ]),
-    )
-    monkeypatch.setattr(
-        cs_coil.data.pf_coil,
-        "r_pf_coil_inner",
-        np.array([
-            5.6944236847973242,
-            5.5985055619292972,
-            17.819978201682968,
-            17.819978201682968,
-            16.385123084628962,
-            16.385123084628962,
-            2.3321999999999998,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-            0.0,
-        ]),
-    )
-    monkeypatch.setattr(tfv, "poisson_steel", 0.29999999999999999)
 
-    r = 2.3
+    r_stress_point = 2.3
+    r_cs_inner = 2.3321999999999998
+    r_cs_outer = 2.88462
+    j_cs = 19311657.760000002
+    b_cs_inner = 13.522197474024983
+    f_poisson_cs_structure = 0.29999999999999999
+    f_a_cs_turn_steel = 0.57874999999999999
+
     s_hoop_exp = 6.737108e8
-    s_hoop = cs_coil.calculate_cs_hoop_stress(r)
+    s_hoop = cs_coil.calculate_cs_hoop_stress(
+        r_stress_point,
+        r_cs_inner,
+        r_cs_outer,
+        j_cs,
+        b_cs_inner,
+        f_poisson_cs_structure,
+        f_a_cs_turn_steel,
+    )
 
     assert pytest.approx(s_hoop) == s_hoop_exp
 

From 435f3a38427debb5179f321bf662abea5cc2ac60 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 20 May 2026 09:27:30 +0100
Subject: [PATCH 23/26] Add CSGeometry dataclass

---
 process/models/pfcoil.py         | 146 ++++++++++++++++++++-----------
 tests/unit/models/test_pfcoil.py |  49 ++++++++++-
 2 files changed, 142 insertions(+), 53 deletions(-)

diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 8cf62c9b4c..181c9401a7 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -1,5 +1,6 @@
 import logging
 import math
+from dataclasses import dataclass
 
 import matplotlib.pyplot as plt
 import numba
@@ -154,24 +155,36 @@ def pfcoil(self):
         # Set up call to MHD scaling routine for coil currents.
         # First break up Central Solenoid solenoid into 'filaments'
 
-        (
-            self.data.pf_coil.z_pf_coil_upper[self.data.pf_coil.n_cs_pf_coils - 1],
-            self.data.pf_coil.z_pf_coil_lower[self.data.pf_coil.n_cs_pf_coils - 1],
-            self.data.pf_coil.r_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1],
-            self.data.pf_coil.r_cs_middle,
-            self.data.pf_coil.z_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1],
-            self.data.pf_coil.r_pf_coil_outer[self.data.pf_coil.n_cs_pf_coils - 1],
-            self.data.pf_coil.r_pf_coil_inner[self.data.pf_coil.n_cs_pf_coils - 1],
-            self.data.pf_coil.a_cs_poloidal,
-            self.data.pf_coil.dz_cs_full,
-            self.data.pf_coil.dr_cs_full,
-        ) = self.cs_coil.calculate_cs_geometry(
+        cs_geometry: CSGeometry = self.cs_coil.calculate_cs_geometry(
             z_tf_inside_half=self.data.build.z_tf_inside_half,
             f_z_cs_tf_internal=self.data.pf_coil.f_z_cs_tf_internal,
             dr_cs=self.data.build.dr_cs,
             dr_bore=self.data.build.dr_bore,
         )
 
+        self.data.pf_coil.z_pf_coil_upper[self.data.pf_coil.n_cs_pf_coils - 1] = (
+            cs_geometry.z_cs_coil_upper
+        )
+        self.data.pf_coil.z_pf_coil_lower[self.data.pf_coil.n_cs_pf_coils - 1] = (
+            cs_geometry.z_cs_coil_lower
+        )
+        self.data.pf_coil.r_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1] = (
+            cs_geometry.r_cs_coil_middle
+        )
+        self.data.pf_coil.r_cs_middle = cs_geometry.r_cs_middle
+        self.data.pf_coil.z_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1] = (
+            cs_geometry.z_cs_coil_middle
+        )
+        self.data.pf_coil.r_pf_coil_outer[self.data.pf_coil.n_cs_pf_coils - 1] = (
+            cs_geometry.r_cs_coil_outer
+        )
+        self.data.pf_coil.r_pf_coil_inner[self.data.pf_coil.n_cs_pf_coils - 1] = (
+            cs_geometry.r_cs_coil_inner
+        )
+        self.data.pf_coil.a_cs_poloidal = cs_geometry.a_cs_poloidal
+        self.data.pf_coil.dz_cs_full = cs_geometry.dz_cs_full
+        self.data.pf_coil.dr_cs_full = cs_geometry.dr_cs_full
+
         # nfxf is the total no of filaments into which the Central Solenoid is split,
         # if present
         if self.data.build.iohcl == 0:
@@ -2966,6 +2979,32 @@ def waveform(self):
             self.data.pf_coil.f_c_pf_cs_peak_time_array[ic, 5] = 0.0e0
 
 
+@dataclass
+class CSGeometry:
+    """Data class to hold the geometry parameters of the Central Solenoid (CS) coil."""
+
+    z_cs_coil_upper: float
+    """Upper Z coordinate of CS coil (m)"""
+    z_cs_coil_lower: float
+    """Lower Z coordinate of CS coil (m)"""
+    r_cs_coil_middle: float
+    """Radial coordinate of CS coil centre (m)"""
+    r_cs_middle: float
+    """Mean radius of CS coil (m)"""
+    z_cs_coil_middle: float
+    """Z coordinate of CS coil centre (m)"""
+    r_cs_coil_outer: float
+    """Outer radius of CS coil (m)"""
+    r_cs_coil_inner: float
+    """Inner radius of CS coil (m)"""
+    a_cs_poloidal: float
+    """Total poloidal cross-sectional area of CS coil (m²)"""
+    dz_cs_full: float
+    """Full height of CS coil (m)"""
+    dr_cs_full: float
+    """Full radial thickness of CS coil (m)"""
+
+
 class CSCoil(Model):
     """Calculate central solenoid coil system parameters."""
 
@@ -2987,7 +3026,7 @@ def calculate_cs_geometry(
         f_z_cs_tf_internal: float,
         dr_cs: float,
         dr_bore: float,
-    ) -> tuple[float, float, float, float, float, float, float, float, float]:
+    ) -> CSGeometry:
         """Calculate the geometry of the Central Solenoid (CS) coil.
 
         Parameters
@@ -3003,16 +3042,18 @@ def calculate_cs_geometry(
 
         Returns
         -------
-        tuple[float, float, float, float, float, float, float, float]
-            Tuple containing:
+        CSGeometry
+            Data class containing the geometry parameters of the CS coil, including:
             - z_cs_coil_upper: Upper Z coordinate of CS coil (m)
             - z_cs_coil_lower: Lower Z coordinate of CS coil (m)
             - r_cs_coil_middle: Radial coordinate of CS coil centre (m)
+            - r_cs_middle: Mean radius of CS coil (m)
             - z_cs_coil_middle: Z coordinate of CS coil centre (m)
             - r_cs_coil_outer: Outer radius of CS coil (m)
             - r_cs_coil_inner: Inner radius of CS coil (m)
             - a_cs_poloidal: Total poloidal cross-sectional area of CS coil (m²)
             - dz_cs_full: Full height of CS coil (m)
+            - dr_cs_full: Full radial thickness of CS coil (m)
         """
         # Central Solenoid mean radius
         r_cs_middle = dr_bore + (0.5e0 * dr_cs)
@@ -3042,17 +3083,17 @@ def calculate_cs_geometry(
         # Total cross-sectional area
         a_cs_poloidal = 2.0e0 * z_cs_inside_half * dr_cs
 
-        return (
-            z_cs_coil_upper,
-            z_cs_coil_lower,
-            r_cs_coil_middle,
-            r_cs_middle,
-            z_cs_coil_middle,
-            r_cs_coil_outer,
-            r_cs_coil_inner,
-            a_cs_poloidal,
-            dz_cs_full,
-            dr_cs_full,
+        return CSGeometry(
+            z_cs_coil_upper=z_cs_coil_upper,
+            z_cs_coil_lower=z_cs_coil_lower,
+            r_cs_coil_middle=r_cs_coil_middle,
+            r_cs_middle=r_cs_middle,
+            z_cs_coil_middle=z_cs_coil_middle,
+            r_cs_coil_outer=r_cs_coil_outer,
+            r_cs_coil_inner=r_cs_coil_inner,
+            a_cs_poloidal=a_cs_poloidal,
+            dz_cs_full=dz_cs_full,
+            dr_cs_full=dr_cs_full,
         )
 
     def calculate_cs_turn_geometry_eu_demo(
@@ -3062,12 +3103,13 @@ def calculate_cs_turn_geometry_eu_demo(
         radius_cs_turn_corners: float,
         f_a_cs_turn_steel: float,
     ) -> tuple[float, float, float, float, float]:
-        """Calculate the geometry of a CS (Central Solenoid) turn using the EU DEMO stadium-shaped model.
+        """Calculate the geometry of a CS (Central Solenoid) turn using the
+        EU DEMO stadium-shaped model.
 
         Parameters
         ----------
          a_cs_turn : float
-             Poloidal area of a CS turn (m^2)
+             Poloidal area of a CS turn (m²)
          f_dr_dz_cs_turn : float
              Length-to-height ratio of the CS turn
          radius_cs_turn_corners : float
@@ -3087,14 +3129,15 @@ def calculate_cs_turn_geometry_eu_demo(
 
         Notes
         -----
-            - The calculation assumes a stadium-shaped cross-section for the CS turn.
-            - If the calculated conduit thickness is negative or too small, it is set to a minimum value of 1 mm.
+        - The calculation assumes a stadium-shaped cross-section for the CS turn.
+        - If the calculated conduit thickness is negative or too small, it is set
+            to a minimum value of 1 mm.
 
         References
         ----------
-            - R. Wesche et al., “Central solenoid winding pack design for DEMO,”
-            Fusion Engineering and Design, vol. 124, pp. 82-85, Apr. 2017,
-            doi: https://doi.org/10.1016/j.fusengdes.2017.04.052.
+        [1] R. Wesche et al., “Central solenoid winding pack design for DEMO,”
+        Fusion Engineering and Design, vol. 124, pp. 82-85, Apr. 2017,
+        doi: https://doi.org/10.1016/j.fusengdes.2017.04.052.
         """
         # Vertical height of CS turn conduit/turn
         dz_cs_turn = (a_cs_turn / f_dr_dz_cs_turn) ** 0.5
@@ -3206,30 +3249,35 @@ def place_cs_filaments(
 
     def ohcalc(self):
         """Routine to perform calculations for the Central Solenoid."""
-        (
-            self.data.pf_coil.z_pf_coil_upper[self.data.pf_coil.n_cs_pf_coils - 1],
-            self.data.pf_coil.z_pf_coil_lower[self.data.pf_coil.n_cs_pf_coils - 1],
-            self.data.pf_coil.r_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1],
-            self.data.pf_coil.r_cs_middle,
-            self.data.pf_coil.z_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1],
-            self.data.pf_coil.r_pf_coil_outer[self.data.pf_coil.n_cs_pf_coils - 1],
-            self.data.pf_coil.r_pf_coil_inner[self.data.pf_coil.n_cs_pf_coils - 1],
-            self.data.pf_coil.a_cs_poloidal,
-            self.data.pf_coil.dz_cs_full,
-            self.data.pf_coil.dr_cs_full,
-        ) = self.calculate_cs_geometry(
+        cs_geometry: CSGeometry = self.calculate_cs_geometry(
             z_tf_inside_half=self.data.build.z_tf_inside_half,
             f_z_cs_tf_internal=self.data.pf_coil.f_z_cs_tf_internal,
             dr_cs=self.data.build.dr_cs,
             dr_bore=self.data.build.dr_bore,
         )
 
-        self.data.pf_coil.r_cs_inner = self.data.pf_coil.r_pf_coil_inner[
-            self.data.pf_coil.n_cs_pf_coils - 1
-        ]
+        self.data.pf_coil.z_pf_coil_upper[self.data.pf_coil.n_cs_pf_coils - 1] = (
+            cs_geometry.z_cs_coil_upper
+        )
+        self.data.pf_coil.z_pf_coil_lower[self.data.pf_coil.n_cs_pf_coils - 1] = (
+            cs_geometry.z_cs_coil_lower
+        )
+        self.data.pf_coil.r_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1] = (
+            cs_geometry.r_cs_coil_middle
+        )
+        self.data.pf_coil.r_cs_middle = cs_geometry.r_cs_middle
+        self.data.pf_coil.z_pf_coil_middle[self.data.pf_coil.n_cs_pf_coils - 1] = (
+            cs_geometry.z_cs_coil_middle
+        )
         self.data.pf_coil.r_cs_outer = self.data.pf_coil.r_pf_coil_outer[
             self.data.pf_coil.n_cs_pf_coils - 1
-        ]
+        ] = cs_geometry.r_cs_coil_outer
+        self.data.pf_coil.r_cs_inner = self.data.pf_coil.r_pf_coil_inner[
+            self.data.pf_coil.n_cs_pf_coils - 1
+        ] = cs_geometry.r_cs_coil_inner
+        self.data.pf_coil.a_cs_poloidal = cs_geometry.a_cs_poloidal
+        self.data.pf_coil.dz_cs_full = cs_geometry.dz_cs_full
+        self.data.pf_coil.dr_cs_full = cs_geometry.dr_cs_full
 
         # Maximum current (MA-turns) in central Solenoid, at either BOP or EOF
         if self.data.pf_coil.j_cs_pulse_start > self.data.pf_coil.j_cs_flat_top_end:
diff --git a/tests/unit/models/test_pfcoil.py b/tests/unit/models/test_pfcoil.py
index 07f268d585..ef77bbbb3d 100644
--- a/tests/unit/models/test_pfcoil.py
+++ b/tests/unit/models/test_pfcoil.py
@@ -21,6 +21,7 @@
 from process.data_structure import tfcoil_variables as tfv
 from process.data_structure.pfcoil_variables import N_PF_COILS_IN_GROUP_MAX
 from process.models.pfcoil import (
+    CSGeometry,
     PFCoil,
     calculate_b_field_at_point,
     fixb,
@@ -2045,17 +2046,57 @@ def test_brookscoil(pfcoil):
     ("z_tf_inside_half", "f_z_cs_tf_internal", "dr_cs", "dr_bore", "expected"),
     [
         # Typical values
-        (2.0, 0.5, 0.3, 0.15, (1.0, -1.0, 0.3, 0.3, 0.0, 0.45, 0.15, 0.6, 2.0, 0.9)),
+        (
+            2.0,
+            0.5,
+            0.3,
+            0.15,
+            CSGeometry(
+                1.0,
+                -1.0,
+                0.3,
+                0.3,
+                0.0,
+                0.44999999999999996,
+                0.14999999999999997,
+                0.6,
+                2.0,
+                0.8999999999999999,
+            ),
+        ),
         # Zero thickness
-        (2.0, 0.5, 0.0, 0.0, (1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0)),
+        (
+            2.0,
+            0.5,
+            0.0,
+            0.0,
+            CSGeometry(1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0),
+        ),
         # Zero fractional height
-        (2.0, 0.0, 0.3, 1.5, (0.0, -0.0, 1.65, 1.65, 0.0, 1.8, 1.5, 0.0, 0.0, 3.6)),
+        (
+            2.0,
+            0.0,
+            0.3,
+            1.5,
+            CSGeometry(
+                0.0,
+                -0.0,
+                1.65,
+                1.65,
+                0.0,
+                1.7999999999999998,
+                1.4999999999999998,
+                0.0,
+                0.0,
+                3.5999999999999996,
+            ),
+        ),
     ],
 )
 def test_calculate_cs_geometry(
     cs_coil, z_tf_inside_half, f_z_cs_tf_internal, dr_cs, dr_bore, expected
 ):
-    result = cs_coil.calculate_cs_geometry(
+    result: CSGeometry = cs_coil.calculate_cs_geometry(
         z_tf_inside_half=z_tf_inside_half,
         f_z_cs_tf_internal=f_z_cs_tf_internal,
         dr_cs=dr_cs,

From a3ffa8302b563d2a660ec2cc09cc2ac9f2ef706f Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 20 May 2026 11:06:59 +0100
Subject: [PATCH 24/26] Refactor CS stress plotting functions for clarity and
 consistency

---
 process/core/io/plot/summary.py | 17 +++++-----
 process/models/pfcoil.py        | 57 ++-------------------------------
 2 files changed, 11 insertions(+), 63 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 780809c265..ed412ef50e 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9837,14 +9837,14 @@ def plot_cs_stress_time_profile(axis: plt.Axes, mfile: MFile, scan: int) -> None
         stress_z_cs_self_midplane_profile / 1e6,
         "o-",
         linewidth=2,
-        markersize=8,
-        label="Midplane Axial Stress",
+        markersize=4,
+        label="$\\sigma_{z}$,Axial Stress",
     )
     axis.set_xlabel("Pulse Time (s)")
-    axis.set_ylabel("Stress (MPa)")
+    axis.set_ylabel("Midplane Axial Stress (MPa)")
     axis.minorticks_on()
     axis.legend(loc="best")
-    axis.set_title("Central Solenoid Stress")
+    axis.set_title("CS Midplane Axial Stress Time Profile")
     axis.grid(True, alpha=0.3)
 
 
@@ -14498,11 +14498,11 @@ def main_plot(
 
     plot_current_profiles_over_time(figs[29].add_subplot(111), m_file, scan)
 
-    plot_cs_stress_time_profile(axis=figs[29].add_subplot(322), mfile=m_file, scan=scan)
+    plot_cs_stress_time_profile(axis=figs[30].add_subplot(431), mfile=m_file, scan=scan)
 
     cs_coil = CSCoil(cs_fatigue=CsFatigue())
     cs_coil.plot_cs_radial_hoop_stress_profile(
-        axis=figs[29].add_subplot(324),
+        axis=figs[30].add_subplot(432),
         mfile=m_file,
         scan=scan,
         j_cs=m_file.get("j_cs_pulse_start", scan=scan),
@@ -14510,7 +14510,7 @@ def main_plot(
     )
 
     cs_coil.plot_cs_radial_stress_profile(
-        axis=figs[29].add_subplot(326),
+        axis=figs[30].add_subplot(433),
         mfile=m_file,
         scan=scan,
         j_cs=m_file.get("j_cs_pulse_start", scan=scan),
@@ -14518,11 +14518,12 @@ def main_plot(
     )
 
     plot_cs_coil_structure(
-        figs[30].add_subplot(121, aspect="equal"), figs[30], m_file, scan
+        figs[30].add_subplot(223, aspect="equal"), figs[30], m_file, scan
     )
     plot_cs_turn_structure(
         figs[30].add_subplot(326, aspect="equal"), figs[30], m_file, scan
     )
+    figs[30].subplots_adjust(wspace=0.3)
 
     plot_first_wall_top_down_cross_section(
         figs[31].add_subplot(221, aspect="equal"), m_file, scan
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 181c9401a7..730dba0bdc 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -4111,56 +4111,6 @@ def calculate_cs_radial_stress(
 
         return hp_term_1 * hp_term_2 - hp_term_3 * hp_term_4
 
-    @staticmethod
-    def plot_stress_time_profile(axis: plt.Axes, mfile: MFile, scan: int):
-        t_plant_pulse_coil_precharge = mfile.get(
-            "t_plant_pulse_coil_precharge", scan=scan
-        )
-        t_plant_pulse_plasma_current_ramp_up = mfile.get(
-            "t_plant_pulse_plasma_current_ramp_up", scan=scan
-        )
-        t_plant_pulse_fusion_ramp = mfile.get("t_plant_pulse_fusion_ramp", scan=scan)
-        t_plant_pulse_burn = mfile.get("t_plant_pulse_burn", scan=scan)
-        t_plant_pulse_plasma_current_ramp_down = mfile.get(
-            "t_plant_pulse_plasma_current_ramp_down", scan=scan
-        )
-
-        # Define a cumulative sum list for each point in the pulse
-        t_steps = np.cumsum([
-            0,
-            t_plant_pulse_coil_precharge,
-            t_plant_pulse_plasma_current_ramp_up,
-            t_plant_pulse_fusion_ramp,
-            t_plant_pulse_burn,
-            t_plant_pulse_plasma_current_ramp_down,
-        ])
-
-        stress_times = t_steps[
-            :6
-        ]  # Get the first 6 time points corresponding to the stress profile
-
-        stress_z_cs_self_midplane_profile = np.zeros(6)
-        for i in range(6):
-            stress_z_cs_self_midplane_profile[i] = mfile.get(
-                f"stress_z_cs_self_midplane_profile[{i}]", scan=scan
-            )
-
-        # Plot stress vs time
-        axis.plot(
-            stress_times,
-            stress_z_cs_self_midplane_profile / 1e6,
-            "o-",
-            linewidth=2,
-            markersize=8,
-            label=r"$\sigma_{{z}}$,Midplane Axial Stress",
-        )
-        axis.set_xlabel("Pulse Time (s)")
-        axis.set_ylabel("Stress (MPa)")
-        axis.minorticks_on()
-        axis.legend(loc="best")
-        axis.set_title("Central Solenoid Stress")
-        axis.grid(True, alpha=0.3)
-
     def plot_cs_radial_hoop_stress_profile(
         self,
         axis: plt.Axes,
@@ -4189,16 +4139,14 @@ def plot_cs_radial_hoop_stress_profile(
         axis.plot(
             radii,
             stress_values / 1e6,
-            "o-",
             linewidth=2,
-            markersize=8,
             label="$\\sigma_{\\theta}$,Hoop Stress",
         )
         axis.set_xlabel("Radial Position (m)")
         axis.set_ylabel("Hoop Stress (MPa)")
         axis.minorticks_on()
         axis.legend(loc="best")
-        axis.set_title("Central Solenoid Radial Hoop Stress Profile")
+        axis.set_title("CS Hoop Stress at BOP")
         axis.grid(True, alpha=0.3)
 
     def plot_cs_radial_stress_profile(
@@ -4228,15 +4176,14 @@ def plot_cs_radial_stress_profile(
         axis.plot(
             radii,
             stress_values / 1e6,
-            "o-",
             linewidth=2,
-            markersize=8,
             label="$\\sigma_{r}$,Radial Stress",
         )
         axis.set_xlabel("Radial Position (m)")
         axis.set_ylabel("Radial Stress (MPa)")
         axis.minorticks_on()
         axis.grid(True, alpha=0.3)
+        axis.set_title("CS Radial Stress at BOP")
         axis.legend(loc="best")
 
 

From 2923f211fff77e7bcdcb6a7ecd53fabd8432b74b Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 21 May 2026 16:17:52 +0100
Subject: [PATCH 25/26] Update docs

---
 .../source/eng-models/central-solenoid.md     | 98 +++++++++++++------
 process/models/pfcoil.py                      | 12 +--
 2 files changed, 76 insertions(+), 34 deletions(-)

diff --git a/documentation/source/eng-models/central-solenoid.md b/documentation/source/eng-models/central-solenoid.md
index 6ab0b083e6..b024fa6278 100644
--- a/documentation/source/eng-models/central-solenoid.md
+++ b/documentation/source/eng-models/central-solenoid.md
@@ -66,21 +66,27 @@ This method calculates the CS geometry parameters. The CS is assumed to be a per
 
 -----------
 
-### Self peak magnetic field | `calculate_cs_self_peak_magnetic_field()`
+### Self peak bore magnetic field | `calculate_cs_bore_magnetic_field()`
 
-The general form for the field at the very centre of the central solenoid bore with uniform current density and rectangular cross-section is given by:
+The self induced peak field at the central bore of a homogenous current density rectangular cross-section solenoid is given by:
 
 $$
-B_0 = J_{\text{CS}}aF(\alpha,\beta)
+B_0 = J_{\text{CS}}r_{\text{CS,inner}}F(\alpha,\beta)
 $$
 
 $$
 F(\alpha,\beta) = \mu_0\beta \ln{\left[\frac{\alpha+\sqrt{\alpha^2+\beta^2}}{1+\sqrt{1+\beta^2}}\right]}
 $$
 
-where $\alpha = \frac{r_{\text{CS,outer}}}{r_{\text{CS,inner}}}$, is the ratio of the outer and inner radii of the solenoid and $\beta = \frac{z_{\text{CS,half}}}{r_{\text{CS,outer}}}$, is the ratio of the solenoid half height to its inboard radius.
+where $\alpha = \frac{r_{\text{CS,outer}}}{r_{\text{CS,inner}}}$, is the ratio of the outer and inner radii of the solenoid and $\beta = \frac{z_{\text{CS,half}}}{r_{\text{CS,inner}}}$, is the ratio of the solenoid half height to its inboard radius.
 
-The peak field at the bore of the central solenoid will not be the same as that felt by the conductors inside the structures. We require to know the peak field on the conductor if we are to design a superconducting central solenoid that has enough margin. Fits to data[^1] for different ranges of $\beta$ have been calculated as follows:
+This is Equation 3.13 from "Case Studies in Superconducting Magnets"[^2].
+
+-------------
+
+### Self peak on coil magnetic field | `calculate_cs_self_peak_magnetic_field()`
+
+The peak field at the bore of the central solenoid will not be the same as that felt by the conductors inside the structures.So wecannot use the bore value directly calculated by [`calculate_cs_bore_magnetic_field()`](#self-peak-bore-magnetic-field--calculate_cs_bore_magnetic_field). We require to know the peak field on the conductor if we are to design a superconducting central solenoid that has enough margin. Fits to data[^1] for different ranges of $\beta$ have been calculated as follows to scale the bore field value by:
 
 - $\beta > 3.0$
 
@@ -123,10 +129,10 @@ The peak field at the bore of the central solenoid will not be the same as that
 
 -----------
 
-### Axial stresses | `calculate_cs_self_peak_midplane_axial_stress()`
+### Peak Axial Stresses | `calculate_cs_self_peak_midplane_axial_stress()`
 
 
-The vertical (axial) force for a "thin-walled" solenoid ($\alpha = 1$) at the midplane is given by[^2]:
+The vertical (axial) force for a "thin-walled" solenoid ($\alpha = 1$) at the midplane where the force is maximum is given by Equation 3.41 in "Case studies in superconducting magnets" [^2]:
 
 $$
 F_{z}(0)=\frac{\mu_0}{2}\left(\frac{N I}{2 \times dz_{\text{half}}}\right) \times \\
@@ -148,64 +154,100 @@ $$
 
 Here $K(k)$ and $E(k)$ are the complete elliptic integrals, respectively of the first and second kinds.
 
+!!! info "Non thin-walled solenoids"
+
+    For solenoids that can be classed as "thick-walled" ($\alpha > 1$), the coil has to be divided radially into several thin-walled solenoids. This is not currently performed though more info can be found in Section 3.5.5 of "Case studies in superconducting magnets" [^2].
 
 The axial compressive force at $z$ in an isolated solenoid increases from 0 at $z  = dz_{\text{half}}$
 to the maximum at the midplane, $F_{z}(0)$.
 
+The axial stress in the steel is given by:
 
+$$
+\sigma_z = \frac{F_z}{f_z A_z}
+$$
 
---------------------------
+where $F_z$ is the axial force, $f_z$ is the fraction of the horizontal cross-section occupied by 
+steel, and $A_z$ is the area of the horizontal cross-section.
 
+The fraction of the horizontal cross-section occupied by steel is calculated assuming that the 
+conductor is square and has a steel jacket with the same thickness on all four sides, giving:
 
-### Hoop stress | `calculate_cs_hoop_stress()`
+$$
+f_z = 0.5.
+$$
+
+
+
+--------------------------
 
 
+### Radial stress | `calculate_cs_radial_stress()`
 
 
-The hoop stress is calculated using equations 4.10 and 4.11 from "Superconducting magnets", Martin N. 
-Wilson (1983).  This is divided by the fraction of the area occupied by steel to obtain the hoop 
-stress in the steel, $\sigma_{hoop}$.
+The self induced radial stress is calculated using Equation 4.11 from "Superconducting magnets" [^1].
 
 $$
-\sigma_{\theta} = \frac{K(2+v)}{3(\alpha+1)}\times \left(\alpha^2+\alpha+1+\frac{\alpha^2}{\epsilon^2}-\epsilon \frac{(1+2v)(\alpha+1)}{(2+v)}\right) \\
-- \frac{M(3+v)}{8}\left(\alpha^2+1+\frac{\alpha^2}{\epsilon^2}-\frac{(1+3v)}{(3+v)}\epsilon^2\right)
+\sigma_{r} = \frac{K(2+v)}{3(\alpha+1)}\times \left(\alpha^2+\alpha+1-\frac{\alpha^2}{\epsilon^2}-\epsilon(\alpha+1)\right) \\
+- \frac{M(3+v)}{8}\left(\alpha^2+1-\frac{\alpha^2}{\epsilon^2}-\epsilon^2\right)
 $$
 
 Where:
 
-- $\epsilon = \frac{r}{r_{CS,inner}}$
-- $\alpha = \frac{r_{CS,outer}}{r_{CS,inner}}$
+- $\epsilon = \frac{r}{r_{\text{CS,inner}}}$
+- $\alpha = \frac{r_{\text{CS,outer}}}{r_{\text{CS,inner}}}$
 
 The terms $K$ and $M$ are given by:
 
 $$
-K = \frac{Ja\left(\alpha B_{\text{a}} - B_{\text{b}}\right)}{(\alpha-1)}
+K = \frac{J r_{\text{CS,inner}}\left(\alpha B_{\text{CS,inner}} - B_{\text{CS,outer}}\right)}{(\alpha-1)}
 $$
 
 $$
-M = \frac{Ja\left(B_{\text{a}} - B_{\text{b}}\right)}{(\alpha-1)}
+M = \frac{J r_{\text{CS,inner}}\left(B_{\text{CS,inner}} - B_{\text{CS,outer}}\right)}{(\alpha-1)}
 $$
 
-For a special infinite solenoid case $B$
+-------------
 
-The axial stress can be calculated using "Case studies in superconducting magnets", Y. Iwasa, p. 
-86, 3.5.2, Special Case 4: Midplane force.  This applies exactly only to a thin-walled solenoid. 
-The axial stress in the steel is given by:
+
+### Hoop stress | `calculate_cs_hoop_stress()`
+
+
+
+
+The hoop stress is calculated using Equation 4.10 from "Superconducting magnets" [^1]. This is divided by the fraction of the area occupied by steel to obtain the hoop 
+stress in the steel, $\sigma_{\theta}$.
 
 $$
-\sigma_z = \frac{F_z}{f_z A_z}
+\sigma_{\theta} = \frac{K(2+v)}{3(\alpha+1)}\times \left(\alpha^2+\alpha+1+\frac{\alpha^2}{\epsilon^2}-\epsilon \frac{(1+2v)(\alpha+1)}{(2+v)}\right) \\
+- \frac{M(3+v)}{8}\left(\alpha^2+1+\frac{\alpha^2}{\epsilon^2}-\frac{(1+3v)}{(3+v)}\epsilon^2\right)
 $$
 
-where $F_z$ is the axial force, $f_z$ is the fraction of the horizontal cross-section occupied by 
-steel, and $A_z$ is the area of the horizontal cross-section.
+Where:
 
-The fraction of the horizontal cross-section occupied by steel is calculated assuming that the 
-conductor is square and has a steel jacket with the same thickness on all four sides, giving:
+- $\epsilon = \frac{r}{r_{\text{CS,inner}}}$
+- $\alpha = \frac{r_{\text{CS,outer}}}{r_{\text{CS,inner}}}$
+
+The terms $K$ and $M$ are given by:
+
+$$
+K = \frac{J r_{\text{CS,inner}}\left(\alpha B_{\text{CS,inner}} - B_{\text{CS,outer}}\right)}{(\alpha-1)}
+$$
 
 $$
-f_z = \frac{f_V}{2}.
+M = \frac{J r_{\text{CS,inner}}\left(B_{\text{CS,inner}} - B_{\text{CS,outer}}\right)}{(\alpha-1)}
 $$
 
+!!! warning "Assumption of outer field, $B_{\text{CS,outer}}$" 
+
+    In this case we currently assume that $B_{\text{CS,outer}} = 0$. This is the same as that for an infinite solenoid. Approximation of the outboard field is currently not performed.
+
+--------------------------
+
+
+
+
+
 The radial stress is neglected. The hoop and axial stresses are combined to give the maximum shear 
 stress, as required by the Tresca stress criterion:
 
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 730dba0bdc..4862b3feec 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -3691,8 +3691,8 @@ def calculate_cs_bore_magnetic_field(
 
         References
         ----------
-            -"Superconducting Magnets", Clarendon Press, Oxford, N.Y., 1983,
-            ISBN 13: 9780198548102
+        [1] "Superconducting Magnets", Clarendon Press, Oxford, N.Y., 1983,
+        ISBN 13: 9780198548102
         """
         beta = dz_cs_half / r_cs_inner
         alpha = r_cs_outer / r_cs_inner
@@ -4089,14 +4089,14 @@ def calculate_cs_radial_stress(
         epsilon = r_stress_point / r_cs_inner
 
         # Field at outer radius of coil [T]
-        # Assume to be 0 for now
-        b_b = 0.0e0
+        # Assume to be 0 for now same as for an infinite solenoid
+        b_cs_outer = 0.0e0
 
         # K term
-        k = ((alpha * b_cs_inner - b_b) * j_cs * r_cs_inner) / (alpha - 1.0e0)
+        k = ((alpha * b_cs_inner - b_cs_outer) * j_cs * r_cs_inner) / (alpha - 1.0e0)
 
         # M term
-        m = ((b_cs_inner - b_b) * j_cs * r_cs_inner) / (alpha - 1.0e0)
+        m = ((b_cs_inner - b_cs_outer) * j_cs * r_cs_inner) / (alpha - 1.0e0)
 
         # calculate hoop stress terms
         hp_term_1 = k * ((2.0e0 + f_poisson_cs_structure) / (3.0e0 * (alpha + 1.0e0)))

From 4b2a42c2abc0f1d34861d88d215439dd286c7df7 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 21 May 2026 18:01:07 +0100
Subject: [PATCH 26/26] :sparkle: Add `a_cs_toroidal` to store the top down
 cross section

---
 process/core/io/plot/summary.py            |   1 +
 process/data_structure/pfcoil_variables.py |   3 +
 process/models/pfcoil.py                   | 104 +++++++++++----------
 tests/unit/models/test_pfcoil.py           |  20 ++--
 4 files changed, 71 insertions(+), 57 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index ed412ef50e..7857209898 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9756,6 +9756,7 @@ def plot_cs_coil_structure(
         f"CS full height: {mfile.get('dz_cs_full', scan=scan):.4f} m\n"
         f"CS full width: {mfile.get('dr_cs_full', scan=scan):.4f} m\n"
         f"CS poloidal area: {mfile.get('a_cs_poloidal', scan=scan):.4f} m$^2$\n"
+        f"CS top-down toroidal area: {mfile.get('a_cs_toroidal', scan=scan):.4f} m$^2$\n"
         f"$N_{{\\text{{turns}}}}:$ {mfile.get('n_pf_coil_turns[n_cs_pf_coils-1]', scan=scan):,.2f}\n"
         f"$I_{{\\text{{peak}}}}:$ {mfile.get('c_pf_cs_coils_peak_ma[n_cs_pf_coils-1]', scan=scan):.3f}$ \\ MA$\n"
         f"$B_{{\\text{{peak}}}}:$ {mfile.get('b_pf_coil_peak[n_cs_pf_coils-1]', scan=scan):.3f}$ \\ T$\n"
diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
index 3599b85f7e..6d3a30f1ea 100644
--- a/process/data_structure/pfcoil_variables.py
+++ b/process/data_structure/pfcoil_variables.py
@@ -106,6 +106,9 @@ class PFCoilData:
     - =1 Hoop + Axial stress
     """
 
+    a_cs_toroidal: float = 0.0
+    """Central solenoid top-down toroidal cross-sectional area (m2)"""
+
     a_cs_poloidal: float = 0.0
     """Central solenoid vertical cross-sectional area (m2)"""
 
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 4862b3feec..223999abb4 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -182,6 +182,7 @@ def pfcoil(self):
             cs_geometry.r_cs_coil_inner
         )
         self.data.pf_coil.a_cs_poloidal = cs_geometry.a_cs_poloidal
+        self.data.pf_coil.a_cs_toroidal = cs_geometry.a_cs_toroidal
         self.data.pf_coil.dz_cs_full = cs_geometry.dz_cs_full
         self.data.pf_coil.dr_cs_full = cs_geometry.dr_cs_full
 
@@ -2056,28 +2057,28 @@ def outpf(self):
             )
             op.ovarre(
                 self.outfile,
-                "Critical superconductor current density at BOP (A/m2)",
+                "Critical superconductor current density at BOP [A/m²]",
                 "(j_cs_conductor_critical_pulse_start)",
                 self.data.pf_coil.j_cs_conductor_critical_pulse_start,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
-                "Critical cable current density at BOP (A/m2)",
+                "Critical cable current density at BOP [A/m²]",
                 "(jcableoh_bop)",
                 self.data.pf_coil.jcableoh_bop,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
-                "Allowable overall current density at BOP (A/m2)",
+                "Allowable overall current density at BOP [A/m²]",
                 "(j_cs_critical_pulse_start)",
                 self.data.pf_coil.j_cs_critical_pulse_start,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
-                "Actual overall current density at BOP (A/m2)",
+                "Actual overall current density at BOP [A/m²]",
                 "(j_cs_pulse_start)",
                 self.data.pf_coil.j_cs_pulse_start,
                 "OP ",
@@ -2085,35 +2086,35 @@ def outpf(self):
             op.oblnkl(self.outfile)
             op.ovarre(
                 self.outfile,
-                "Maximum field at End Of Flattop (T)",
+                "Maximum field at End Of Flattop [T]",
                 "(b_cs_peak_flat_top_end)",
                 self.data.pf_coil.b_cs_peak_flat_top_end,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
-                "Critical superconductor current density at EOF (A/m2)",
+                "Critical superconductor current density at EOF [A/m²]",
                 "(j_cs_conductor_critical_flat_top_end)",
                 self.data.pf_coil.j_cs_conductor_critical_flat_top_end,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
-                "Critical cable current density at EOF (A/m2)",
+                "Critical cable current density at EOF [A/m²]",
                 "(jcableoh_eof)",
                 self.data.pf_coil.jcableoh_eof,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
-                "Allowable overall current density at EOF (A/m2)",
+                "Allowable overall current density at EOF [A/m²]",
                 "(j_cs_critical_flat_top_end)",
                 self.data.pf_coil.j_cs_critical_flat_top_end,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
-                "Actual overall current density at EOF (A/m2)",
+                "Actual overall current density at EOF [A/m²]",
                 "(j_cs_flat_top_end)",
                 self.data.pf_coil.j_cs_flat_top_end,
             )
@@ -2135,55 +2136,62 @@ def outpf(self):
             # MDK add self.data.build.dr_cs, self.data.build.dr_bore and self.data.build.dr_cs_tf_gap as they can be iteration variables
             op.ovarre(
                 self.outfile,
-                "CS inside radius (m)",
+                "CS inside radius [m]",
                 "(dr_bore)",
                 self.data.build.dr_bore,
             )
-            op.ovarre(self.outfile, "CS thickness (m)", "(dr_cs)", self.data.build.dr_cs)
+            op.ovarre(self.outfile, "CS thickness [m]", "(dr_cs)", self.data.build.dr_cs)
             op.ovarre(
                 self.outfile,
-                "Gap between central solenoid and TF coil (m)",
+                "Gap between central solenoid and TF coil [m]",
                 "(dr_cs_tf_gap)",
                 self.data.build.dr_cs_tf_gap,
             )
             op.ovarre(
                 self.outfile,
-                "CS overall cross-sectional area (m2)",
+                "CS total poloidal cross-sectional area [m²]",
                 "(a_cs_poloidal)",
                 self.data.pf_coil.a_cs_poloidal,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
-                "CS radial middle (m)",
+                "CS total top-down toroidal cross-sectional area [m²]",
+                "(a_cs_toroidal)",
+                self.data.pf_coil.a_cs_toroidal,
+                "OP ",
+            )
+            op.ovarre(
+                self.outfile,
+                "CS radial middle [m]",
                 "(r_cs_middle)",
                 self.data.pf_coil.r_cs_middle,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
-                "CS radial inner (m)",
+                "CS radial inner [m]",
                 "(r_cs_inner)",
                 self.data.pf_coil.r_cs_inner,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
-                "CS radial outer (m)",
+                "CS radial outer [m]",
                 "(r_cs_outer)",
                 self.data.pf_coil.r_cs_outer,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
-                "CS conductor+void cross-sectional area (m2)",
+                "CS conductor+void cross-sectional area [m²]",
                 "(a_cs_cable_space)",
                 self.data.pf_coil.a_cs_cable_space,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
-                "   CS conductor cross-sectional area (m2)",
+                "   CS conductor cross-sectional area [m²]",
                 "(a_cs_cable_space*(1-f_a_cs_void))",
                 self.data.pf_coil.a_cs_cable_space
                 * (1.0e0 - self.data.pf_coil.f_a_cs_void),
@@ -2191,14 +2199,14 @@ def outpf(self):
             )
             op.ovarre(
                 self.outfile,
-                "   CS void cross-sectional area (m2)",
+                "   CS void cross-sectional area [m²]",
                 "(a_cs_cable_space*f_a_cs_void)",
                 self.data.pf_coil.a_cs_cable_space * self.data.pf_coil.f_a_cs_void,
                 "OP ",
             )
             op.ovarre(
                 self.outfile,
-                "CS steel cross-sectional area (m2)",
+                "CS steel cross-sectional area [m²]",
                 "(a_cs_steel_poloidal)",
                 self.data.pf_coil.a_cs_steel_poloidal,
                 "OP ",
@@ -2999,6 +3007,8 @@ class CSGeometry:
     """Inner radius of CS coil (m)"""
     a_cs_poloidal: float
     """Total poloidal cross-sectional area of CS coil (m²)"""
+    a_cs_toroidal: float
+    """Total top-down toroidal cross-sectional area of CS coil (m²)"""
     dz_cs_full: float
     """Full height of CS coil (m)"""
     dr_cs_full: float
@@ -3044,16 +3054,17 @@ def calculate_cs_geometry(
         -------
         CSGeometry
             Data class containing the geometry parameters of the CS coil, including:
-            - z_cs_coil_upper: Upper Z coordinate of CS coil (m)
-            - z_cs_coil_lower: Lower Z coordinate of CS coil (m)
-            - r_cs_coil_middle: Radial coordinate of CS coil centre (m)
-            - r_cs_middle: Mean radius of CS coil (m)
-            - z_cs_coil_middle: Z coordinate of CS coil centre (m)
-            - r_cs_coil_outer: Outer radius of CS coil (m)
-            - r_cs_coil_inner: Inner radius of CS coil (m)
-            - a_cs_poloidal: Total poloidal cross-sectional area of CS coil (m²)
-            - dz_cs_full: Full height of CS coil (m)
-            - dr_cs_full: Full radial thickness of CS coil (m)
+            - z_cs_coil_upper: Upper Z coordinate of CS coil [m]
+            - z_cs_coil_lower: Lower Z coordinate of CS coil [m]
+            - r_cs_coil_middle: Radial coordinate of CS coil centre [m]
+            - r_cs_middle: Mean radius of CS coil [m]
+            - z_cs_coil_middle: Z coordinate of CS coil centre [m]
+            - r_cs_coil_outer: Outer radius of CS coil [m]
+            - r_cs_coil_inner: Inner radius of CS coil [m]
+            - a_cs_poloidal: Total poloidal cross-sectional area of CS coil [m²]
+            - a_cs_toroidal: Total top-down toroidal cross-sectional area of CS coil [m²]
+            - dz_cs_full: Full height of CS coil [m]
+            - dr_cs_full: Full radial thickness of CS coil [m]
         """
         # Central Solenoid mean radius
         r_cs_middle = dr_bore + (0.5e0 * dr_cs)
@@ -3080,9 +3091,12 @@ def calculate_cs_geometry(
         # Full radial thickness of CS coil
         dr_cs_full = 2 * r_cs_coil_outer
 
-        # Total cross-sectional area
+        # Total poloidal cross-sectional area [m²]
         a_cs_poloidal = 2.0e0 * z_cs_inside_half * dr_cs
 
+        # Total top-down toroidal cross-sectional area [m²]
+        a_cs_toroidal = np.pi * (r_cs_coil_outer**2 - r_cs_coil_inner**2)
+
         return CSGeometry(
             z_cs_coil_upper=z_cs_coil_upper,
             z_cs_coil_lower=z_cs_coil_lower,
@@ -3092,6 +3106,7 @@ def calculate_cs_geometry(
             r_cs_coil_outer=r_cs_coil_outer,
             r_cs_coil_inner=r_cs_coil_inner,
             a_cs_poloidal=a_cs_poloidal,
+            a_cs_toroidal=a_cs_toroidal,
             dz_cs_full=dz_cs_full,
             dr_cs_full=dr_cs_full,
         )
@@ -3276,6 +3291,7 @@ def ohcalc(self):
             self.data.pf_coil.n_cs_pf_coils - 1
         ] = cs_geometry.r_cs_coil_inner
         self.data.pf_coil.a_cs_poloidal = cs_geometry.a_cs_poloidal
+        self.data.pf_coil.a_cs_toroidal = cs_geometry.a_cs_toroidal
         self.data.pf_coil.dz_cs_full = cs_geometry.dz_cs_full
         self.data.pf_coil.dr_cs_full = cs_geometry.dr_cs_full
 
@@ -3431,7 +3447,6 @@ def ohcalc(self):
                 f_a_cs_turn_steel=self.data.pf_coil.f_a_cs_turn_steel,
             )
 
-            # New calculation from Y. Iwasa for axial stress
             (
                 self.data.pf_coil.stress_z_cs_self_peak_midplane,
                 self.data.pf_coil.forc_z_cs_self_peak_midplane,
@@ -3439,14 +3454,12 @@ def ohcalc(self):
                 r_cs_outer=self.data.pf_coil.r_pf_coil_outer[
                     self.data.pf_coil.n_cs_pf_coils - 1
                 ],
-                r_cs_inner=self.data.pf_coil.r_pf_coil_inner[
-                    self.data.pf_coil.n_cs_pf_coils - 1
-                ],
                 dz_cs_half=self.data.pf_coil.dz_cs_full / 2.0,
                 c_cs_peak=self.data.pf_coil.c_pf_cs_coils_peak_ma[
                     self.data.pf_coil.n_cs_pf_coils - 1
                 ]
                 * 1.0e6,
+                a_cs_toroidal=self.data.pf_coil.a_cs_toroidal,
             )
 
             # Allowable (hoop) stress (Pa) alstroh
@@ -3847,22 +3860,22 @@ def output_cs_structure(self):
     def calculate_cs_self_peak_midplane_axial_stress(
         self,
         r_cs_outer: float,
-        r_cs_inner: float,
         dz_cs_half: float,
         c_cs_peak: float,
+        a_cs_toroidal: float,
     ) -> tuple[float, float]:
         """Calculate axial stress and axial force for the central solenoid.
 
         Parameters
         ----------
         r_cs_outer:
-            Outer radius of the central solenoid (m).
-        r_cs_inner:
-            Inner radius of the central solenoid (m).
+            Outer radius of the central solenoid [m].
         dz_cs_half:
-            Half-height of the central solenoid (m).
+            Half-height of the central solenoid [m].
         c_cs_peak:
-            Peak CS coil current (A).
+            Peak CS coil current [A].
+        a_cs_toroidal:
+            Total top-down toroidal area of the CS [m²].
 
         Returns
         -------
@@ -3918,12 +3931,9 @@ def calculate_cs_self_peak_midplane_axial_stress(
         # calculate axial force [N]
         forc_z_cs_self_peak_midplane = axial_term_1 * (axial_term_2 - axial_term_3)
 
-        # axial area [m2]
-        area_ax = np.pi * (r_cs_outer**2 - r_cs_inner**2)
-
         # Calculate unsmeared axial stress
         # Average axial stress at the interface of each half of the coil
-        s_axial = forc_z_cs_self_peak_midplane / (0.5 * area_ax)
+        s_axial = forc_z_cs_self_peak_midplane / (0.5 * a_cs_toroidal)
 
         return s_axial, forc_z_cs_self_peak_midplane
 
@@ -3935,9 +3945,6 @@ def calculate_cs_self_midplane_axial_stress_time_profile(
                 r_cs_outer=self.data.pf_coil.r_pf_coil_outer[
                     self.data.pf_coil.n_cs_pf_coils - 1
                 ],
-                r_cs_inner=self.data.pf_coil.r_pf_coil_inner[
-                    self.data.pf_coil.n_cs_pf_coils - 1
-                ],
                 dz_cs_half=self.data.pf_coil.dz_cs_full / 2.0,
                 c_cs_peak=(
                     self.data.pf_coil.c_pf_coil_turn[
@@ -3947,6 +3954,7 @@ def calculate_cs_self_midplane_axial_stress_time_profile(
                         self.data.pf_coil.n_cs_pf_coils - 1
                     ]
                 ),
+                a_cs_toroidal=self.data.pf_coil.a_cs_toroidal,
             )
             self.data.pf_coil.stress_z_cs_self_midplane_profile[time] = stress_value
 
diff --git a/tests/unit/models/test_pfcoil.py b/tests/unit/models/test_pfcoil.py
index ef77bbbb3d..dda3c02ee2 100644
--- a/tests/unit/models/test_pfcoil.py
+++ b/tests/unit/models/test_pfcoil.py
@@ -2060,6 +2060,7 @@ def test_brookscoil(pfcoil):
                 0.44999999999999996,
                 0.14999999999999997,
                 0.6,
+                0.5654866776461627,
                 2.0,
                 0.8999999999999999,
             ),
@@ -2070,7 +2071,7 @@ def test_brookscoil(pfcoil):
             0.5,
             0.0,
             0.0,
-            CSGeometry(1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0),
+            CSGeometry(1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0),
         ),
         # Zero fractional height
         (
@@ -2087,6 +2088,7 @@ def test_brookscoil(pfcoil):
                 1.7999999999999998,
                 1.4999999999999998,
                 0.0,
+                3.1101767270538945,
                 0.0,
                 3.5999999999999996,
             ),
@@ -2359,12 +2361,12 @@ def test_calculate_cs_turn_geometry_eu_demo_output_consistency(
 def test_calculate_cs_self_peak_midplane_axial_stress_basic(cs_coil):
     # Typical values for a CS coil
     r_cs_outer = 2.0
-    r_cs_inner = 1.5
+    a_cs_toroidal = 5.497787144
     dz_cs_half = 1.0
     c_cs_peak = 1.2e7  # 12 MA
 
     s_axial, force_axial = cs_coil.calculate_cs_self_peak_midplane_axial_stress(
-        r_cs_outer, r_cs_inner, dz_cs_half, c_cs_peak
+        r_cs_outer, dz_cs_half, c_cs_peak, a_cs_toroidal
     )
 
     # Check actual output numbers
@@ -2379,12 +2381,12 @@ def test_calculate_cs_self_peak_midplane_axial_stress_basic(cs_coil):
 def test_calculate_cs_self_peak_midplane_axial_stress_zero_current(cs_coil):
     # Zero current should yield zero force and stress
     r_cs_outer = 2.0
-    r_cs_inner = 1.5
+    a_cs_toroidal = 5.497787144
     dz_cs_half = 1.0
     c_cs_peak = 0.0
 
     s_axial, force_axial = cs_coil.calculate_cs_self_peak_midplane_axial_stress(
-        r_cs_outer, r_cs_inner, dz_cs_half, c_cs_peak
+        r_cs_outer, dz_cs_half, c_cs_peak, a_cs_toroidal
     )
 
     assert pytest.approx(s_axial) == 0.0
@@ -2394,12 +2396,12 @@ def test_calculate_cs_self_peak_midplane_axial_stress_zero_current(cs_coil):
 def test_calculate_cs_self_peak_midplane_axial_stress_extreme_geometry(cs_coil):
     # Very thin coil, large dz_cs_half
     r_cs_outer = 1.01
-    r_cs_inner = 1.0
+    a_cs_toroidal = 0.06314601234
     dz_cs_half = 10.0
     c_cs_peak = 1.0e7
 
     s_axial, force_axial = cs_coil.calculate_cs_self_peak_midplane_axial_stress(
-        r_cs_outer, r_cs_inner, dz_cs_half, c_cs_peak
+        r_cs_outer, dz_cs_half, c_cs_peak, a_cs_toroidal
     )
 
     # Check actual output numbers
@@ -2413,12 +2415,12 @@ def test_calculate_cs_self_peak_midplane_axial_stress_extreme_geometry(cs_coil):
 def test_calculate_cs_self_peak_midplane_axial_stress_invalid(cs_coil):
     # Negative dz_cs_half should still return floats, but may be positive force
     r_cs_outer = 2.0
-    r_cs_inner = 1.5
+    a_cs_toroidal = 5.497787144
     dz_cs_half = 2.0
     c_cs_peak = 1.2e7
 
     s_axial, force_axial = cs_coil.calculate_cs_self_peak_midplane_axial_stress(
-        r_cs_outer, r_cs_inner, dz_cs_half, c_cs_peak
+        r_cs_outer, dz_cs_half, c_cs_peak, a_cs_toroidal
     )
 
     # Check actual output numbers