diff --git a/documentation/source/eng-models/central-solenoid.md b/documentation/source/eng-models/central-solenoid.md
index 354aa3ae10..43095c9b8b 100644
--- a/documentation/source/eng-models/central-solenoid.md
+++ b/documentation/source/eng-models/central-solenoid.md
@@ -33,7 +33,7 @@ This method calculates the CS geometry parameters. The CS is assumed to be a per
 2. The half height of the CS is set relative to that of the inside height of the TF and can be scaled by changing the input value of `f_z_cs_tf_internal`:
 
     $$
-    \overbrace{z_{\text{CS,half}}}^{\texttt{z_cs_inside_half}} = \overbrace{z_{\text{TF,inside-half}}}^{\texttt{z_tf_inside_half}} \times \texttt{f_z_cs_tf_internal}
+    \overbrace{z_{\text{CS,half}}}^{\texttt{z_cs_half}} = \overbrace{z_{\text{TF,inside-half}}}^{\texttt{z_tf_inside_half}} \times \texttt{f_z_cs_tf_internal}
     $$
 
 3. The full height of the CS is thus simply given by
@@ -42,6 +42,17 @@ This method calculates the CS geometry parameters. The CS is assumed to be a per
     \overbrace{dz_{\text{CS}}}^{\texttt{dz_cs_full}} = z_{\text{CS,half}} \times 2
     $$
 
+    Along with the upper and lower dimensions:
+
+    $$
+    z_{\text{CS,upper}} = z_{\text{CS,half}}
+    $$
+
+    $$
+    z_{\text{CS,lower}} = -z_{\text{CS,upper}}
+    $$
+
+
 4. The outboard edge of the CS is given by:
 
 
@@ -49,12 +60,75 @@ This method calculates the CS geometry parameters. The CS is assumed to be a per
     r_{\text{CS,outer}} = r_{\text{CS,middle}} + \frac{dr_{\text{CS}}}{2}
     $$
 
-5. The full poloidal cross-sectional area is given by:
+5. The inboard edge of the CS is given by:
+
+    $$
+    r_{\text{CS,inner}} = r_{\text{CS,outer}} - dr_{\text{CS}}
+    $$
+
+6. The full radial width of the CS is given by:
+
+    $$
+    dr_{\text{CS,full}}  = 2 \times r_{\text{CS,outer}}
+    $$
+
+
+7. The full poloidal cross-sectional area is given by:
 
     $$
     \overbrace{A_{\text{CS,poloidal}}}^{\texttt{a_cs_poloidal}} = 2 \times dr_{\text{CS}} \times dz_{\text{CS}}
     $$
 
+8. The full top-down toroidal cross-sectional area is given by:
+
+    $$
+    \overbrace{A_{\text{CS,toroidal}}}^{\texttt{a_cs_toroidal}} = \pi \left(r_{\text{CS,outer}}^2 - r_{\text{CS,inner}}^2 \right)
+    $$
+
+------------
+
+### EU-DEMO Turn Geometry | `calculate_cs_turn_geometry_eu_demo()`
+
+This turn superconducting turn strucutre for the CS assumes a rectangular turn shape with a "stadium" shaped cable area[^eu_demo_turn].
+
+
+![CS turn layout](../eng-models/images/cs_eu_demo_turn.PNG "CS EU-DEMO like turn")
+
+The turn geometry is calculated as follows:
+
+1. The vertical height of the turn is given by:
+
+    $$
+    dz_{\text{CS,turn}} = \left(\frac{A_{\text{CS,turn}}}{\texttt{f_dr_dz_cs_turn}}\right)^2
+    $$
+
+    $\texttt{f_dr_dz_cs_turn}$ is the intended length to height ratio of the turn
+
+2. The radial width or length of the turn is now given by:
+
+    $$
+    dr_{\text{CS,turn}} = \texttt{f_dr_dz_cs_turn} \times dz_{\text{CS,turn}}
+    $$
+
+3. The radius of the corners of the cable space is given by:
+
+    $$
+    r_{\text{CS,cable space corner}} = - \frac{dr_{\text{CS,turn}}-dz_{\text{CS,turn}}}{\pi} \\
+    + \sqrt{\left(\frac{dr_{\text{CS,turn}}-dz_{\text{CS,turn}}}{\pi}\right)^2+ \frac{dr_{\text{CS,turn}}dz_{\text{CS,turn}}(4-\pi)r_{\text{CS,turn corners}}^2 - (A_{\text{CS,turn}}\times \texttt{f_a_cs_turn_steel})}{\pi}}
+    $$
+
+4. The thickness of the conduit around the cable space is given by:
+
+    $$
+    dz_{\text{CS,turn conduit}} = \frac{dz_{\text{CS,turn}}}{2} - r_{\text{CS,cable space corner}}
+    $$
+
+    In this model the vertical and radial conduit thicknesses are equal so:
+
+    $$
+    dr_{\text{CS,turn conduit}}  = dz_{\text{CS,turn conduit}}
+    $$
+
 
 ------------
 
@@ -66,21 +140,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.
+
+This is Equation 3.13 from "Case Studies in Superconducting Magnets"[^2].
+
+-------------
 
-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:
+### 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 +203,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,41 +228,99 @@ $$
 
 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:
+
+$$
+f_z = 0.5.
+$$
+
 
 
 --------------------------
 
 
-### Hoop stress | `hoop_stress()`
+### Radial stress | `calculate_cs_radial_stress()`
 
 
+The self induced radial stress is calculated using Equation 4.11 from "Superconducting magnets" [^1].
 
+$$
+\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)
+$$
 
-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}$.
+Where:
 
-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:
+- $\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:
 
 $$
-\sigma_z = \frac{F_z}{f_z A_z}
+K = \frac{J r_{\text{CS,inner}}\left(\alpha B_{\text{CS,inner}} - B_{\text{CS,outer}}\right)}{(\alpha-1)}
 $$
 
-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.
+$$
+M = \frac{J r_{\text{CS,inner}}\left(B_{\text{CS,inner}} - B_{\text{CS,outer}}\right)}{(\alpha-1)}
+$$
+
+-------------
+
+
+### Hoop stress | `calculate_cs_hoop_stress()`
+
+
 
-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:
 
+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_{\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)
 $$
-f_z = \frac{f_V}{2}.
+
+Where:
+
+- $\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)}
+$$
+
+$$
+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:
@@ -320,4 +458,5 @@ constraints (26 and 27) are activated.
 
 [^1]: M. N. Wilson, Superconducting Magnets. Oxford University Press, USA, 1983, ISBN 13: 9780198548102
 [^2]: Case Studies in Superconducting Magnets. Boston, MA: Springer US, 2009. doi: https://doi.org/10.1007/b112047.
+[^eu_demo_turn]: 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.
 ‌
\ No newline at end of file
diff --git a/documentation/source/eng-models/images/cs_eu_demo_turn.PNG b/documentation/source/eng-models/images/cs_eu_demo_turn.PNG
new file mode 100644
index 0000000000..12da080491
Binary files /dev/null and b/documentation/source/eng-models/images/cs_eu_demo_turn.PNG differ
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 ac1e0951d5..7857209898 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,
@@ -47,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,
@@ -9754,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"
@@ -9835,14 +9838,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)
 
 
@@ -14496,7 +14499,24 @@ 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)
+    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[30].add_subplot(432),
+        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[30].add_subplot(433),
+        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
@@ -14504,6 +14524,7 @@ def main_plot(
     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/data_structure/build_variables.py b/process/data_structure/build_variables.py
index 75c9c2dece..1b159db5ef 100644
--- a/process/data_structure/build_variables.py
+++ b/process/data_structure/build_variables.py
@@ -138,10 +138,10 @@ class BuildData:
     """
 
     dr_cs: float = 0.811
-    """Central solenoid thickness (m) (`iteration variable 16`)"""
+    """Central solenoid radial thickness (m) (`iteration variable 16`)"""
 
     dr_cs_precomp: float = 0.0
-    """CS coil precompression structure thickness (m)"""
+    """CS coil precompression structure radial thickness (m)"""
 
     rbld: float = 0.0
     """sum of thicknesses to the major radius (m)"""
diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
index 720e813ae5..7cbc9943b3 100644
--- a/process/data_structure/pfcoil_variables.py
+++ b/process/data_structure/pfcoil_variables.py
@@ -42,7 +42,13 @@ 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_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] = 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)"""
@@ -100,21 +106,33 @@ 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)"""
 
+    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)"""
 
-    awpoh: float = 0.0
+    a_cs_cable_space: float = 0.0
     """central solenoid conductor+void area with area of steel subtracted (m2)"""
 
+    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 = 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)"""
 
@@ -347,6 +365,21 @@ class PFCoilData:
     r_cs_middle: float = 0.0
     """radius to the centre of the central solenoid (m)"""
 
+    r_cs_inner: float = 0.0
+    """Inner radius of the central solenoid (m)"""
+
+    r_cs_outer: float = 0.0
+    """Outer radius of the central solenoid (m)"""
+
+    z_cs_upper: float = 0.0
+    """z location of upper end of central solenoid (m)"""
+
+    z_cs_lower: float = 0.0
+    """z location of lower end of central solenoid (m)"""
+
+    z_cs_middle: float = 0.0
+    """z location of middle of central solenoid (m)"""
+
     dz_cs_full: float = 0.0
     """Full height of the central solenoid (m)"""
 
@@ -395,6 +428,9 @@ class PFCoilData:
     temp_cs_superconductor_margin: float = 0.0
     """Central solenoid temperature margin (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 bc68a495e3..2b08ff3071 100644
--- a/process/models/costs/costs.py
+++ b/process/models/costs/costs.py
@@ -1716,7 +1716,7 @@ def acc2222(self):
                 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
+                        * 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[
@@ -1746,7 +1746,7 @@ def acc2222(self):
             if self.data.pf_coil.i_pf_conductor == 0:
                 costpfcu = (
                     self.data.costs.uccu
-                    * self.data.pf_coil.awpoh
+                    * 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[
@@ -1758,7 +1758,7 @@ def acc2222(self):
                 # 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
+                    * 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
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index 71f73afe9c..dc93b37bfd 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -1,6 +1,8 @@
 import logging
 import math
+from dataclasses import dataclass
 
+import matplotlib.pyplot as plt
 import numba
 import numpy as np
 from scipy import optimize
@@ -10,6 +12,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
@@ -152,24 +155,37 @@ 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_cs_upper = 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_cs_lower = 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_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.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
+
         # nfxf is the total no of filaments into which the Central Solenoid is split,
         # if present
         if self.data.build.iohcl == 0:
@@ -2000,36 +2016,34 @@ def outpf(self):
                 self.data.pf_coil.i_cs_superconductor,
             )
 
-            if self.data.pf_coil.i_cs_superconductor == 1:
-                op.ocmmnt(self.outfile, "  (ITER Nb3Sn critical surface model)")
-            elif self.data.pf_coil.i_cs_superconductor == 2:
-                op.ocmmnt(self.outfile, "  (Bi-2212 high temperature superconductor)")
-            elif self.data.pf_coil.i_cs_superconductor == 3:
-                op.ocmmnt(self.outfile, "  (NbTi)")
-            elif self.data.pf_coil.i_cs_superconductor == 4:
-                op.ocmmnt(
-                    self.outfile,
-                    "  (ITER Nb3Sn critical surface model, user-defined parameters)",
-                )
-            elif self.data.pf_coil.i_cs_superconductor == 5:
-                op.ocmmnt(self.outfile, " (WST Nb3Sn critical surface model)")
-            elif self.data.pf_coil.i_cs_superconductor == 6:
-                op.ocmmnt(self.outfile, " (REBCO HTS)")
-            elif self.data.pf_coil.i_cs_superconductor == 7:
-                op.ocmmnt(
-                    self.outfile,
-                    " (Durham Ginzburg-Landau critical surface model for Nb-Ti)",
-                )
-            elif self.data.pf_coil.i_cs_superconductor == 8:
-                op.ocmmnt(
-                    self.outfile,
-                    " (Durham Ginzburg-Landau critical surface model for REBCO)",
-                )
-            elif self.data.pf_coil.i_cs_superconductor == 9:
-                op.ocmmnt(
-                    self.outfile,
-                    " (Hazelton experimental data + Zhai conceptual model for REBCO)",
-                )
+            op.ocmmnt(
+                self.outfile,
+                f"Superconductor used: "
+                f"{SuperconductorModel(self.data.pf_coil.i_cs_superconductor).full_name}",
+            )
+            op.oblnkl(self.outfile)
+            op.ovarre(
+                self.outfile,
+                "CS superconductor operating temperature (K)",
+                "(temp_cs_superconductor_operating)",
+                self.data.pf_coil.temp_cs_superconductor_operating,
+            )
+            op.ovarre(
+                self.outfile,
+                "CS temperature margin (K)",
+                "(temp_cs_superconductor_margin)",
+                self.data.pf_coil.temp_cs_superconductor_margin,
+                "OP ",
+            )
+            op.ovarre(
+                self.outfile,
+                "Minimum permitted temperature margin (K)",
+                "(temp_cs_superconductor_margin_min)",
+                tfv.temp_cs_superconductor_margin_min,
+            )
+
+            op.oblnkl(self.outfile)
+            op.ocmmnt(self.outfile, "----------------------------")
 
             op.osubhd(self.outfile, "Central Solenoid Current Density Limits :")
             op.ovarre(
@@ -2041,28 +2055,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 ",
@@ -2070,35 +2084,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,
             )
@@ -2117,72 +2131,9 @@ def outpf(self):
                     self.data.pf_coil.z_pf_cs_current_filaments[i],
                 )
             op.oblnkl(self.outfile)
-            # 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)",
-                "(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,
-                "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)",
-                "(a_cs_poloidal)",
-                self.data.pf_coil.a_cs_poloidal,
-                "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 conductor+void cross-sectional area (m2)",
-                "(awpoh)",
-                self.data.pf_coil.awpoh,
-                "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),
-                "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,
-                "OP ",
-            )
-            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,
-                "OP ",
-            )
-            op.ovarre(
-                self.outfile,
-                "CS steel area fraction",
-                "(f_a_cs_turn_steel)",
-                self.data.pf_coil.f_a_cs_turn_steel,
-            )
-            if self.data.pf_coil.i_cs_stress == 1:
-                op.ocmmnt(self.outfile, "Hoop + axial stress considered")
-            else:
-                op.ocmmnt(self.outfile, "Only hoop stress considered")
+
+            op.ocmmnt(self.outfile, "----------------------------")
+            op.osubhd(self.outfile, "CS Stresses:")
 
             op.ovarin(
                 self.outfile,
@@ -2190,6 +2141,10 @@ def outpf(self):
                 "(i_cs_stress)",
                 self.data.pf_coil.i_cs_stress,
             )
+            if self.data.pf_coil.i_cs_stress == 1:
+                op.ocmmnt(self.outfile, "Hoop + axial stress considered")
+            else:
+                op.ocmmnt(self.outfile, "Only hoop stress considered")
             op.ovarre(
                 self.outfile,
                 "Allowable stress in CS steel (Pa)",
@@ -2199,8 +2154,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_inner)",
+                self.data.pf_coil.stress_hoop_cs_inner,
                 "OP ",
             )
             op.ovarre(
@@ -2230,52 +2185,8 @@ def outpf(self):
                 "(str_cs_con_res)",
                 tfv.str_cs_con_res,
             )
-            op.ovarre(
-                self.outfile,
-                "Copper fraction in strand",
-                "(fcuohsu)",
-                self.data.pf_coil.fcuohsu,
-            )
-            # If REBCO material is used, print copperaoh_m2
-            if self.data.pf_coil.i_cs_superconductor in {6, 8, 9}:
-                op.ovarre(
-                    self.outfile,
-                    "CS current/copper area (A/m2)",
-                    "(copperaoh_m2)",
-                    rcv.copperaoh_m2,
-                )
-                op.ovarre(
-                    self.outfile,
-                    "Max CS current/copper area (A/m2)",
-                    "(copperaoh_m2_max)",
-                    rcv.copperaoh_m2_max,
-                )
 
-            op.ovarre(
-                self.outfile,
-                "Void (coolant) fraction in conductor",
-                "(f_a_cs_void)",
-                self.data.pf_coil.f_a_cs_void,
-            )
-            op.ovarre(
-                self.outfile,
-                "Helium coolant temperature (K)",
-                "(tftmp)",
-                tfv.tftmp,
-            )
-            op.ovarre(
-                self.outfile,
-                "CS temperature margin (K)",
-                "(temp_cs_superconductor_margin)",
-                self.data.pf_coil.temp_cs_superconductor_margin,
-                "OP ",
-            )
-            op.ovarre(
-                self.outfile,
-                "Minimum permitted temperature margin (K)",
-                "(temp_cs_superconductor_margin_min)",
-                tfv.temp_cs_superconductor_margin_min,
-            )
+            op.oblnkl(self.outfile)
             # only output CS fatigue model for pulsed reactor design
             if pv.f_c_plasma_inductive > 0.0e-4:
                 op.ovarre(
@@ -2302,42 +2213,7 @@ def outpf(self):
                     "(t_crack_radial)",
                     self.data.cs_fatigue.t_crack_radial,
                 )
-                op.ovarre(
-                    self.outfile,
-                    "CS turn area (m)",
-                    "(a_cs_turn)",
-                    self.data.pf_coil.a_cs_turn,
-                )
-                op.ovarre(
-                    self.outfile,
-                    "CS turn length (m)",
-                    "(dr_cs_turn)",
-                    self.data.pf_coil.dr_cs_turn,
-                )
-                op.ovarre(
-                    self.outfile,
-                    "CS turn internal cable space radius (m)",
-                    "(radius_cs_turn_cable_space)",
-                    self.data.pf_coil.radius_cs_turn_cable_space,
-                )
-                op.ovarre(
-                    self.outfile,
-                    "CS turn width (m)",
-                    "(dz_cs_turn)",
-                    self.data.pf_coil.dz_cs_turn,
-                )
-                op.ovarre(
-                    self.outfile,
-                    "CS structural vertical thickness (m)",
-                    "(dz_cs_turn_conduit)",
-                    self.data.cs_fatigue.dz_cs_turn_conduit,
-                )
-                op.ovarre(
-                    self.outfile,
-                    "CS structural radial thickness (m)",
-                    "(dr_cs_turn_conduit)",
-                    self.data.cs_fatigue.dr_cs_turn_conduit,
-                )
+
                 op.ovarre(
                     self.outfile,
                     "Allowable number of cycles till CS fracture",
@@ -2421,6 +2297,10 @@ def outpf(self):
         else:
             op.ocmmnt(self.outfile, "Resistive central solenoid")
 
+        op.oblnkl(self.outfile)
+        op.ocmmnt(self.outfile, "----------------------------")
+        op.oblnkl(self.outfile)
+
         if self.data.pf_coil.i_pf_conductor == 0:
             op.oblnkl(self.outfile)
             op.ocmmnt(self.outfile, "Superconducting PF coils")
@@ -2434,7 +2314,7 @@ def outpf(self):
 
             op.ocmmnt(
                 self.outfile,
-                f"  -> {SuperconductorModel(self.data.pf_coil.i_pf_superconductor).full_name}",
+                f"Superconductor used: {SuperconductorModel(self.data.pf_coil.i_pf_superconductor).full_name}",
             )
 
             op.ovarre(
@@ -2487,7 +2367,7 @@ def outpf(self):
         op.osubhd(self.outfile, "Geometry of PF coils, central solenoid and plasma:")
         op.write(
             self.outfile,
-            "coil\t\t\tR(m)\t\tZ(m)\t\tdR(m)\t\tdZ(m)\t\tturns",
+            "Coil\t\t\tR(m)\t\tZ(m)\t\tdR(m)\t\tdZ(m)\t\tturns",
         )
         op.oblnkl(self.outfile)
 
@@ -2622,15 +2502,19 @@ def outpf(self):
             f"Plasma\t\t\t{pv.rmajor:.2e}\t0.0e0\t\t{2.0e0 * pv.rminor:.2e}\t{2.0e0 * pv.rminor * pv.kappa:.2e}\t1.0e0",
         )
 
+        op.oblnkl(self.outfile)
+        op.ocmmnt(self.outfile, "----------------------------")
+        op.oblnkl(self.outfile)
+
         op.osubhd(self.outfile, "PF Coil Information at Peak Current:")
 
         op.write(
             self.outfile,
-            "coil\tcurrent\t\tallowed J\tactual J\tJ\t\tcond. mass\tsteel mass\tfield",
+            r"Coil   Peak Current   Critical J  Peak J        Critical J Ratio  Cond. Mass     Steel Mass     Field",
         )
         op.write(
             self.outfile,
-            "\t(MA)\t\t(A/m2)\t\t(A/m2)\t\f_temp_plasma_ion_electron\t\t(kg)\t\t(kg)\t\t(T)",
+            "       (MA)          (A/m²)      (A/m²)          (-)            (kg)            (kg)           (T)",
         )
 
         op.oblnkl(self.outfile)
@@ -2640,12 +2524,18 @@ def outpf(self):
             if self.data.pf_coil.i_pf_conductor == 0:
                 op.write(
                     self.outfile,
-                    f"PF {k}\t{self.data.pf_coil.c_pf_cs_coils_peak_ma[k]:.2e}\t{self.data.pf_coil.j_pf_wp_critical[k]:.2e}\t{self.data.pf_coil.j_pf_coil_wp_peak[k]:.2e}\t{self.data.pf_coil.j_pf_coil_wp_peak[k] / self.data.pf_coil.j_pf_wp_critical[k]:.2e}\t{self.data.pf_coil.m_pf_coil_conductor[k]:.2e}\t{self.data.pf_coil.m_pf_coil_structure[k]:.2e}\t{self.data.pf_coil.b_pf_coil_peak[k]:.2e}",
+                    f"PF {k}   {self.data.pf_coil.c_pf_cs_coils_peak_ma[k]:{'.3e' if self.data.pf_coil.c_pf_cs_coils_peak_ma[k] >= 0 else '.2e'}}    "
+                    f"{self.data.pf_coil.j_pf_wp_critical[k]:{'.3e' if self.data.pf_coil.j_pf_wp_critical[k] >= 0 else '.2e'}}    "
+                    f"{self.data.pf_coil.j_pf_coil_wp_peak[k]:{'.3e' if self.data.pf_coil.j_pf_coil_wp_peak[k] >= 0 else '.2e'}}      "
+                    f"{self.data.pf_coil.j_pf_coil_wp_peak[k] / self.data.pf_coil.j_pf_wp_critical[k]:{'.3e' if (self.data.pf_coil.j_pf_coil_wp_peak[k] / self.data.pf_coil.j_pf_wp_critical[k]) >= 0 else '.2e'}}          "
+                    f"{self.data.pf_coil.m_pf_coil_conductor[k]:{'.3e' if self.data.pf_coil.m_pf_coil_conductor[k] >= 0 else '.2e'}}          "
+                    f"{self.data.pf_coil.m_pf_coil_structure[k]:{'.3e' if self.data.pf_coil.m_pf_coil_structure[k] >= 0 else '.2e'}}          "
+                    f"{self.data.pf_coil.b_pf_coil_peak[k]:{'.3e' if self.data.pf_coil.b_pf_coil_peak[k] >= 0 else '.2e'}}",
                 )
             else:
                 op.write(
                     self.outfile,
-                    f"PF {k}\t{self.data.pf_coil.c_pf_cs_coils_peak_ma[k]:.2e}\t-1.0e0\t{self.data.pf_coil.j_pf_coil_wp_peak[k]:.2e}\t1.0e0\t{self.data.pf_coil.m_pf_coil_conductor[k]:.2e}\t{self.data.pf_coil.m_pf_coil_structure[k]:.2e}\t{self.data.pf_coil.b_pf_coil_peak[k]:.2e}\t",
+                    f"PF {k}   {self.data.pf_coil.c_pf_cs_coils_peak_ma[k]:.2e}\t-1.0e0\t{self.data.pf_coil.j_pf_coil_wp_peak[k]:.2e}\t1.0e0\t{self.data.pf_coil.m_pf_coil_conductor[k]:.2e}\t{self.data.pf_coil.m_pf_coil_structure[k]:.2e}\t{self.data.pf_coil.b_pf_coil_peak[k]:.2e}\t",
                 )
 
         # Central Solenoid, if present
@@ -2654,12 +2544,18 @@ def outpf(self):
                 # Issue #328
                 op.write(
                     self.outfile,
-                    f"CS\t\t{self.data.pf_coil.c_pf_cs_coils_peak_ma[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}\t{self.data.pf_coil.j_pf_wp_critical[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}\t{max(abs(self.data.pf_coil.j_cs_pulse_start), abs(self.data.pf_coil.j_cs_flat_top_end)):.2e}\t{max(abs(self.data.pf_coil.j_cs_pulse_start), abs(self.data.pf_coil.j_cs_flat_top_end)) / self.data.pf_coil.j_pf_wp_critical[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}\t{self.data.pf_coil.m_pf_coil_conductor[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}\t{self.data.pf_coil.m_pf_coil_structure[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}\t{self.data.pf_coil.b_pf_coil_peak[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}",
+                    f"CS    {self.data.pf_coil.c_pf_cs_coils_peak_ma[self.data.pf_coil.n_cs_pf_coils - 1]:{'.3e' if self.data.pf_coil.c_pf_cs_coils_peak_ma[self.data.pf_coil.n_cs_pf_coils - 1] >= 0 else '.2e'}}    "
+                    f"{self.data.pf_coil.j_pf_wp_critical[self.data.pf_coil.n_cs_pf_coils - 1]:{'.3e' if self.data.pf_coil.j_pf_wp_critical[self.data.pf_coil.n_cs_pf_coils - 1] >= 0 else '.2e'}}    "
+                    f"{max(abs(self.data.pf_coil.j_cs_pulse_start), abs(self.data.pf_coil.j_cs_flat_top_end)):{'.3e' if max(abs(self.data.pf_coil.j_cs_pulse_start), abs(self.data.pf_coil.j_cs_flat_top_end)) >= 0 else '.2e'}}      "
+                    f"{max(abs(self.data.pf_coil.j_cs_pulse_start), abs(self.data.pf_coil.j_cs_flat_top_end)) / self.data.pf_coil.j_pf_wp_critical[self.data.pf_coil.n_cs_pf_coils - 1]:{'.3e' if (max(abs(self.data.pf_coil.j_cs_pulse_start), abs(self.data.pf_coil.j_cs_flat_top_end)) / self.data.pf_coil.j_pf_wp_critical[self.data.pf_coil.n_cs_pf_coils - 1]) >= 0 else '.2e'}}          "
+                    f"{self.data.pf_coil.m_pf_coil_conductor[self.data.pf_coil.n_cs_pf_coils - 1]:{'.3e' if self.data.pf_coil.m_pf_coil_conductor[self.data.pf_coil.n_cs_pf_coils - 1] >= 0 else '.2e'}}          "
+                    f"{self.data.pf_coil.m_pf_coil_structure[self.data.pf_coil.n_cs_pf_coils - 1]:{'.3e' if self.data.pf_coil.m_pf_coil_structure[self.data.pf_coil.n_cs_pf_coils - 1] >= 0 else '.2e'}}          "
+                    f"{self.data.pf_coil.b_pf_coil_peak[self.data.pf_coil.n_cs_pf_coils - 1]:{'.3e' if self.data.pf_coil.b_pf_coil_peak[self.data.pf_coil.n_cs_pf_coils - 1] >= 0 else '.2e'}}",
                 )
             else:
                 op.write(
                     self.outfile,
-                    f"CS\t\t{self.data.pf_coil.c_pf_cs_coils_peak_ma[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}\t-1.0e0\t{max(abs(self.data.pf_coil.j_cs_pulse_start)):.2e}\t{abs(self.data.pf_coil.j_cs_flat_top_end):.2e}\t1.0e0\t{self.data.pf_coil.m_pf_coil_conductor[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}\t{self.data.pf_coil.m_pf_coil_structure[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}\t{self.data.pf_coil.b_pf_coil_peak[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}",
+                    f"CS     {self.data.pf_coil.c_pf_cs_coils_peak_ma[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}\t-1.0e0\t{max(abs(self.data.pf_coil.j_cs_pulse_start), abs(self.data.pf_coil.j_cs_flat_top_end)):.2e}\t1.0e0\t{self.data.pf_coil.m_pf_coil_conductor[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}\t{self.data.pf_coil.m_pf_coil_structure[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}\t{self.data.pf_coil.b_pf_coil_peak[self.data.pf_coil.n_cs_pf_coils - 1]:.2e}",
                 )
 
         # Miscellaneous totals
@@ -2676,6 +2572,10 @@ def outpf(self):
             + f"{self.data.pf_coil.m_pf_coil_conductor_total:.2e}\t{self.data.pf_coil.m_pf_coil_structure_total:.2e}",
         )
 
+        op.oblnkl(self.outfile)
+        op.ocmmnt(self.outfile, "----------------------------")
+        op.oblnkl(self.outfile)
+
         op.osubhd(self.outfile, "PF coil current scaling information :")
         op.ovarre(
             self.outfile,
@@ -2696,22 +2596,23 @@ def outvolt(self):
         """
         op.oheadr(self.outfile, "Volt Second Consumption")
 
-        op.write(self.outfile, "\t" * 3 + "volt-sec\t\t\tvolt-sec\t\tvolt-sec")
-        op.write(self.outfile, "\t" * 3 + "start-up\t\t\t_burn\t\t\ttotal")
+        pf = self.data.pf_coil
+        op.write(self.outfile, "\t" * 6 + "volt-sec\t\t\tvolt-sec\t\tvolt-sec")
+        op.write(self.outfile, "\t" * 6 + "start-up\t\t\tburn\t\t\t\ttotal")
         op.write(
             self.outfile,
-            f"PF coils:\t\t{self.data.pf_coil.vs_pf_coils_total_ramp:.2f}\t\t\t\t{self.data.pf_coil.vs_pf_coils_total_burn:.2f}\t\t\t{self.data.pf_coil.vs_pf_coils_total_pulse:.2f}",
+            f"PF coils:\t\t{pf.vs_pf_coils_total_ramp:.2f}\t\t\t\t{pf.vs_pf_coils_total_burn:.2f}\t\t\t{pf.vs_pf_coils_total_pulse:.2f}",
         )
         op.write(
             self.outfile,
-            f"CS coil:\t\t{self.data.pf_coil.vs_cs_ramp:.2f}\t\t\t\t{self.data.pf_coil.vs_cs_burn:.2f}\t\t\t{self.data.pf_coil.vs_cs_total_pulse:.2f}",
+            f"CS coil:\t\t{pf.vs_cs_ramp:.2f}\t\t\t\t{pf.vs_cs_burn:.2f}\t\t\t{pf.vs_cs_total_pulse:.2f}",
         )
         op.write(
-            self.outfile, "\t" * 3 + "-" * 7 + "\t" * 4 + "-" * 7 + "\t" * 3 + "-" * 7
+            self.outfile, "\t" * 6 + "-" * 7 + "\t" * 4 + "-" * 7 + "\t" * 3 + "-" * 7
         )
         op.write(
             self.outfile,
-            f"Total:\t\t\t{self.data.pf_coil.vs_cs_pf_total_ramp:.2f}\t\t\t\t{self.data.pf_coil.vs_cs_pf_total_burn:.2f}\t\t\t{self.data.pf_coil.vs_cs_pf_total_pulse:.2f}",
+            f"Total:\t\t\t{pf.vs_cs_pf_total_ramp:.2f}\t\t\t\t{pf.vs_cs_pf_total_burn:.2f}\t\t\t{pf.vs_cs_pf_total_pulse:.2f}",
         )
 
         op.oblnkl(self.outfile)
@@ -2731,8 +2632,7 @@ def outvolt(self):
         )
 
         op.osubhd(self.outfile, "Summary of volt-second consumption by circuit (Wb):")
-
-        op.write(self.outfile, "circuit\t\t\tBOP\t\t\tBOF\t\tEOF")
+        op.write(self.outfile, "Circuit\t\t\tBOP\t\t\tBOF\t\tEOF")
         op.oblnkl(self.outfile)
 
         for k in range(self.data.pf_coil.nef):
@@ -2741,9 +2641,10 @@ def outvolt(self):
                 f"\t{k}\t\t\t{self.data.pf_coil.vsdum[k, 0]:.3f}\t\t\t{self.data.pf_coil.vsdum[k, 1]:.3f}\t\t{self.data.pf_coil.vsdum[k, 2]:.3f}",
             )
 
+        n_cs = self.data.pf_coil.n_cs_pf_coils - 1
         op.write(
             self.outfile,
-            f"\tCS coil\t\t\t{self.data.pf_coil.vsdum[self.data.pf_coil.n_cs_pf_coils - 1, 0]:.3f}\t\t\t{self.data.pf_coil.vsdum[self.data.pf_coil.n_cs_pf_coils - 1, 1]:.3f}\t\t{self.data.pf_coil.vsdum[self.data.pf_coil.n_cs_pf_coils - 1, 2]:.3f}",
+            f"\tCS coil\t\t\t{self.data.pf_coil.vsdum[n_cs, 0]:.3f}\t\t\t{self.data.pf_coil.vsdum[n_cs, 1]:.3f}\t\t{self.data.pf_coil.vsdum[n_cs, 2]:.3f}",
         )
 
         op.oshead(self.outfile, "Waveforms")
@@ -2949,6 +2850,52 @@ 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²)"""
+    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
+    """Full radial thickness of CS coil (m)"""
+
+
+@dataclass
+class CSEUDEMOTurnGeometry:
+    """Data class to hold the geometry parameters of a CS turn using the
+    EU DEMO stadium-shaped model.
+    """
+
+    dz_cs_turn: float
+    """Vertical thickness of CS turn (m)"""
+    dr_cs_turn: float
+    """Length/ radial width of CS turn (m)"""
+    radius_cs_turn_cable_space: float
+    """Radius of CS turn cable space corners (m)"""
+    dr_cs_turn_conduit: float
+    """Radial thickness of CS turn conduit (m)"""
+    dz_cs_turn_conduit: float
+    """Vertical thickness of CS turn conduit (m)"""
+
+
 class CSCoil(Model):
     """Calculate central solenoid coil system parameters."""
 
@@ -2970,7 +2917,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
@@ -2986,27 +2933,30 @@ def calculate_cs_geometry(
 
         Returns
         -------
-        tuple[float, float, float, float, float, float, float, float]
-            Tuple containing:
-            - 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)
-            - 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)
+        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²]
+            - 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)
 
         # Scale the CS height relative to TF bore height
-        z_cs_inside_half = z_tf_inside_half * f_z_cs_tf_internal
+        z_cs_half = z_tf_inside_half * f_z_cs_tf_internal
 
-        dz_cs_full = 2.0e0 * z_cs_inside_half  # Full height of CS coil
+        dz_cs_full = 2.0e0 * z_cs_half  # Full height of CS coil
 
         # Z coordinates of CS coil edges
-        z_cs_coil_upper = z_cs_inside_half
+        z_cs_coil_upper = z_cs_half
         z_cs_coil_lower = -z_cs_coil_upper
 
         # (R,Z) coordinates of coil centre
@@ -3022,20 +2972,24 @@ def calculate_cs_geometry(
         # Full radial thickness of CS coil
         dr_cs_full = 2 * r_cs_coil_outer
 
-        # 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,
+        # Total poloidal cross-sectional area [m²]
+        a_cs_poloidal = dz_cs_full * 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,
+            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,
+            a_cs_toroidal=a_cs_toroidal,
+            dz_cs_full=dz_cs_full,
+            dr_cs_full=dr_cs_full,
         )
 
     def calculate_cs_turn_geometry_eu_demo(
@@ -3044,13 +2998,14 @@ def calculate_cs_turn_geometry_eu_demo(
         f_dr_dz_cs_turn: float,
         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.
+    ) -> CSEUDEMOTurnGeometry:
+        """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
@@ -3060,8 +3015,7 @@ def calculate_cs_turn_geometry_eu_demo(
 
         Returns
         -------
-        :
-            Tuple containing:
+        CSEUDEMOTurnGeometry
             - dz_cs_turn: Depth/width of CS turn conduit (m)
             - dr_cs_turn: Length of CS turn conduit (m)
             - radius_cs_turn_cable_space: Radius of CS turn cable space (m)
@@ -3070,14 +3024,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
@@ -3108,12 +3063,12 @@ def calculate_cs_turn_geometry_eu_demo(
             dr_cs_turn_conduit = 1.0e-3
             logger.error("CS turn conduit radial thickness < 1 mm, kludged to 1 mm")
 
-        return (
-            dz_cs_turn,
-            dr_cs_turn,
-            radius_cs_turn_cable_space,
-            dr_cs_turn_conduit,
-            dz_cs_turn_conduit,
+        return CSEUDEMOTurnGeometry(
+            dz_cs_turn=dz_cs_turn,
+            dr_cs_turn=dr_cs_turn,
+            radius_cs_turn_cable_space=radius_cs_turn_cable_space,
+            dr_cs_turn_conduit=dr_cs_turn_conduit,
+            dz_cs_turn_conduit=dz_cs_turn_conduit,
         )
 
     def place_cs_filaments(
@@ -3189,24 +3144,37 @@ 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.z_cs_upper = 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_cs_lower = 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_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.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
+
         # 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
@@ -3248,17 +3216,25 @@ def ohcalc(self):
             / self.data.pf_coil.n_pf_coil_turns[self.data.pf_coil.n_cs_pf_coils - 1]
         )
 
-        (
-            self.data.pf_coil.dz_cs_turn,
-            self.data.pf_coil.dr_cs_turn,
-            self.data.pf_coil.radius_cs_turn_cable_space,
-            self.data.cs_fatigue.dr_cs_turn_conduit,
-            self.data.cs_fatigue.dz_cs_turn_conduit,
-        ) = self.calculate_cs_turn_geometry_eu_demo(
-            a_cs_turn=self.data.pf_coil.a_cs_turn,
-            f_dr_dz_cs_turn=self.data.pf_coil.f_dr_dz_cs_turn,
-            radius_cs_turn_corners=self.data.pf_coil.radius_cs_turn_corners,
-            f_a_cs_turn_steel=self.data.pf_coil.f_a_cs_turn_steel,
+        eu_demo_turn_geometry: CSEUDEMOTurnGeometry = (
+            self.calculate_cs_turn_geometry_eu_demo(
+                a_cs_turn=self.data.pf_coil.a_cs_turn,
+                f_dr_dz_cs_turn=self.data.pf_coil.f_dr_dz_cs_turn,
+                radius_cs_turn_corners=self.data.pf_coil.radius_cs_turn_corners,
+                f_a_cs_turn_steel=self.data.pf_coil.f_a_cs_turn_steel,
+            )
+        )
+
+        self.data.pf_coil.dz_cs_turn = eu_demo_turn_geometry.dz_cs_turn
+        self.data.pf_coil.dr_cs_turn = eu_demo_turn_geometry.dr_cs_turn
+        self.data.pf_coil.radius_cs_turn_cable_space = (
+            eu_demo_turn_geometry.radius_cs_turn_cable_space
+        )
+        self.data.cs_fatigue.dr_cs_turn_conduit = (
+            eu_demo_turn_geometry.dr_cs_turn_conduit
+        )
+        self.data.cs_fatigue.dz_cs_turn_conduit = (
+            eu_demo_turn_geometry.dz_cs_turn_conduit
         )
 
         # Non-steel area void fraction for coolant
@@ -3267,10 +3243,10 @@ 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(
+        b_cs_self_peak_flat_top_end = 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
@@ -3285,17 +3261,20 @@ def ohcalc(self):
 
         # Peak field due to other PF coils plus plasma
         timepoint = 5
-        _bri, _bro, bzi, bzo = peak_b_field_at_pf_coil(
+        _, _, bzi, bzo = peak_b_field_at_pf_coil(
             n_coil=self.data.pf_coil.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)
+        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)
+
+        self.data.pf_coil.b_cs_self_outer_midplane = 0.0
+
         bohco = abs(bzo)
 
         # Peak field at the Beginning-Of-Pulse (BOP)
@@ -3315,7 +3294,7 @@ def ohcalc(self):
             )
         )
         timepoint = 2
-        _bri, _bro, bzi, bzo = peak_b_field_at_pf_coil(
+        _, _, bzi, bzo = peak_b_field_at_pf_coil(
             n_coil=self.data.pf_coil.n_cs_pf_coils,
             n_coil_group=99,
             t_b_field_peak=timepoint,
@@ -3340,11 +3319,22 @@ 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.r_pf_coil_inner[self.data.pf_coil.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=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
+                ],
+                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
             (
                 self.data.pf_coil.stress_z_cs_self_peak_midplane,
                 self.data.pf_coil.forc_z_cs_self_peak_midplane,
@@ -3352,14 +3342,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
@@ -3373,7 +3361,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,
+                    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,
@@ -3384,24 +3372,24 @@ def ohcalc(self):
             # equation is used for Central Solenoid stress
 
             # Area of steel in Central Solenoid
-            areaspf = (
+            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(
-                        self.data.pf_coil.sig_hoop
+                        self.data.pf_coil.stress_hoop_cs_inner
                         - self.data.pf_coil.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(0.0e0 - self.data.pf_coil.stress_hoop_cs_inner),
                 )
             else:
                 self.data.pf_coil.s_shear_cs_peak = max(
-                    abs(self.data.pf_coil.sig_hoop - 0.0e0),
+                    abs(self.data.pf_coil.stress_hoop_cs_inner - 0.0e0),
                     abs(0.0e0 - 0.0e0),
-                    abs(0.0e0 - self.data.pf_coil.sig_hoop),
+                    abs(0.0e0 - self.data.pf_coil.stress_hoop_cs_inner),
                 )
 
             # Thickness of hypothetical steel cylinders assumed to encase the CS along
@@ -3409,17 +3397,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.a_cs_steel_poloidal
                 / self.data.pf_coil.z_pf_coil_upper[self.data.pf_coil.n_cs_pf_coils - 1]
             )
 
         else:
-            areaspf = 0.0e0  # Resistive Central Solenoid - no steel needed
+            self.data.pf_coil.a_cs_steel_poloidal = (
+                0.0e0  # Resistive Central Solenoid - no steel needed
+            )
             self.data.pf_coil.pfcaseth[self.data.pf_coil.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
+            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]
@@ -3427,21 +3417,25 @@ def ohcalc(self):
         )
 
         # Non-steel cross-sectional area
-        self.data.pf_coil.awpoh = self.data.pf_coil.a_cs_poloidal - areaspf
+        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 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 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 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
+                self.data.pf_coil.a_cs_cable_space
                 * (1.0e0 - self.data.pf_coil.f_a_cs_void)
                 * 2.0e0
                 * np.pi
@@ -3452,7 +3446,7 @@ def ohcalc(self):
             self.data.pf_coil.m_pf_coil_conductor[
                 self.data.pf_coil.n_cs_pf_coils - 1
             ] = (
-                self.data.pf_coil.awpoh
+                self.data.pf_coil.a_cs_cable_space
                 * (1.0e0 - self.data.pf_coil.f_a_cs_void)
                 * 2.0e0
                 * np.pi
@@ -3479,13 +3473,13 @@ def ohcalc(self):
                             self.data.pf_coil.n_cs_pf_coils - 1
                         ]
                     )
-                    / self.data.pf_coil.awpoh
+                    / self.data.pf_coil.a_cs_cable_space
                 )
                 * 1.0e6,
                 self.data.pf_coil.i_cs_superconductor,
                 tfv.fhts,
                 tfv.str_cs_con_res,
-                tfv.tftmp,
+                self.data.pf_coil.temp_cs_superconductor_operating,
                 tfv.bcritsc,
                 tfv.tcritsc,
             )
@@ -3502,7 +3496,9 @@ def ohcalc(self):
                 )
 
             self.data.pf_coil.j_cs_critical_flat_top_end = (
-                jcritwp * self.data.pf_coil.awpoh / self.data.pf_coil.a_cs_poloidal
+                jcritwp
+                * self.data.pf_coil.a_cs_cable_space
+                / self.data.pf_coil.a_cs_poloidal
             )
 
             # Allowable coil overall current density at BOP
@@ -3522,19 +3518,21 @@ def ohcalc(self):
                             self.data.pf_coil.n_cs_pf_coils - 1
                         ]
                     )
-                    / self.data.pf_coil.awpoh
+                    / self.data.pf_coil.a_cs_cable_space
                 )
                 * 1.0e6,
                 self.data.pf_coil.i_cs_superconductor,
                 tfv.fhts,
                 tfv.str_cs_con_res,
-                tfv.tftmp,
+                self.data.pf_coil.temp_cs_superconductor_operating,
                 tfv.bcritsc,
                 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
+                jcritwp
+                * 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]
@@ -3567,14 +3565,14 @@ def ohcalc(self):
             )
         self.calculate_cs_self_midplane_axial_stress_time_profile()
 
-    def calculate_cs_self_peak_magnetic_field(
-        self,
+    @staticmethod
+    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
         ----------
@@ -3590,19 +3588,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,
-            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
 
-        # 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
@@ -3612,6 +3608,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
+        ----------
+        [1] 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)
@@ -3644,22 +3682,123 @@ def calculate_cs_self_peak_magnetic_field(
 
         return b_cs_peak
 
-    def output_cs_structure(self):
+    def output_cs_structure(self) -> None:
+        """Outputs the central solenoid structure parameters to the output file."""
         op.oheadr(self.outfile, "Central Solenoid Structure")
 
-        op.ocmmnt(self.outfile, "CS turn structure")
+        op.osubhd(self.outfile, "CS coil geometry:")
+
+        op.ovarre(
+            self.outfile,
+            "Inner radius of the CS coil [m]",
+            "(r_cs_inner)",
+            self.data.pf_coil.r_cs_inner,
+            "OP ",
+        )
+        op.ovarre(
+            self.outfile,
+            "Middle radius of the CS coil [m]",
+            "(r_cs_middle)",
+            self.data.pf_coil.r_cs_middle,
+            "OP ",
+        )
+        op.ovarre(
+            self.outfile,
+            "Outer radius of the CS coil [m]",
+            "(r_cs_outer)",
+            self.data.pf_coil.r_cs_outer,
+            "OP ",
+        )
+        op.oblnkl(self.outfile)
+        op.ovarre(
+            self.outfile,
+            "Full radial extent of the CS coil [m]",
+            "(dr_cs_full)",
+            self.data.pf_coil.dr_cs_full,
+            "OP ",
+        )
+        op.ovarre(
+            self.outfile,
+            "Radial thickness of the CS coil [m]",
+            "(dr_cs)",
+            self.data.build.dr_cs,
+            "OP ",
+        )
+        op.ovarre(
+            self.outfile,
+            "Radial thickness of the CS bore [m]",
+            "(dr_bore)",
+            self.data.build.dr_bore,
+            "OP ",
+        )
+        op.oblnkl(self.outfile)
+        op.ovarre(
+            self.outfile,
+            "Vertical top of the CS coil [m]",
+            "(z_cs_upper)",
+            self.data.pf_coil.z_cs_upper,
+            "OP ",
+        )
+        op.ovarre(
+            self.outfile,
+            "Vertical middle of the CS coil [m]",
+            "(z_cs_middle)",
+            self.data.pf_coil.z_cs_middle,
+            "OP ",
+        )
+        op.ovarre(
+            self.outfile,
+            "Vertical bottom of the CS coil [m]",
+            "(z_cs_lower)",
+            self.data.pf_coil.z_cs_lower,
+            "OP ",
+        )
         op.oblnkl(self.outfile)
 
         op.ovarre(
             self.outfile,
-            "Poloidal area of a CS turn [m^2]",
+            "Full vertical extent of the CS coil [m]",
+            "(dz_cs_full)",
+            self.data.pf_coil.dz_cs_full,
+            "OP ",
+        )
+        op.ovarre(
+            self.outfile,
+            "Central solenoid to TF coil internal edge height [m]",
+            "(f_z_cs_tf_internal)",
+            self.data.pf_coil.f_z_cs_tf_internal,
+            "OP ",
+        )
+        op.oblnkl(self.outfile)
+        op.ovarre(
+            self.outfile,
+            "CS poloidal cross-sectional area [m²]",
+            "(a_cs_poloidal)",
+            self.data.pf_coil.a_cs_poloidal,
+            "OP ",
+        )
+        op.ovarre(
+            self.outfile,
+            "CS total top-down toroidal cross-sectional area [m²]",
+            "(a_cs_toroidal)",
+            self.data.pf_coil.a_cs_toroidal,
+            "OP ",
+        )
+
+        op.oblnkl(self.outfile)
+        op.ocmmnt(self.outfile, "----------------------------")
+        op.osubhd(self.outfile, "CS turn structure:")
+
+        op.ovarre(
+            self.outfile,
+            "Poloidal area of a CS turn [m²]",
             "(a_cs_turn)",
             self.data.pf_coil.a_cs_turn,
             "OP ",
         )
         op.ovarre(
             self.outfile,
-            "Radial width a CS turn [m^2]",
+            "Radial width a CS turn [m²]",
             "(dz_cs_turn)",
             self.data.pf_coil.dz_cs_turn,
             "OP ",
@@ -3706,26 +3845,87 @@ def output_cs_structure(self):
             self.data.pf_coil.radius_cs_turn_corners,
             "OP ",
         )
+        op.ovarre(
+            self.outfile,
+            "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 [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),
+            "OP ",
+        )
+        op.ovarre(
+            self.outfile,
+            "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 [m²]",
+            "(a_cs_steel_poloidal)",
+            self.data.pf_coil.a_cs_steel_poloidal,
+            "OP ",
+        )
+        op.ovarre(
+            self.outfile,
+            "CS steel area fraction",
+            "(f_a_cs_turn_steel)",
+            self.data.pf_coil.f_a_cs_turn_steel,
+        )
+        op.ovarre(
+            self.outfile,
+            "Copper fraction in strand",
+            "(fcuohsu)",
+            self.data.pf_coil.fcuohsu,
+        )
+        # If REBCO material is used, print copperaoh_m2
+        if self.data.pf_coil.i_cs_superconductor in {6, 8, 9}:
+            op.ovarre(
+                self.outfile,
+                "CS current/copper area (A/m²)",
+                "(copperaoh_m2)",
+                rcv.copperaoh_m2,
+            )
+            op.ovarre(
+                self.outfile,
+                "Max CS current/copper area (A/m²)",
+                "(copperaoh_m2_max)",
+                rcv.copperaoh_m2_max,
+            )
+
+            op.ovarre(
+                self.outfile,
+                "Void (coolant) fraction in conductor",
+                "(f_a_cs_void)",
+                self.data.pf_coil.f_a_cs_void,
+            )
 
     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
         -------
@@ -3781,12 +3981,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
 
@@ -3798,9 +3995,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[
@@ -3810,10 +4004,20 @@ 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
 
-    def hoop_stress(self, r):
+    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,
+        f_a_cs_turn_steel: float,
+    ) -> float:
         """Calculation of hoop stress of central solenoid.
 
         This routine calculates the hoop stress of the central solenoid
@@ -3821,43 +4025,52 @@ def hoop_stress(self, r):
 
         Parameters
         ----------
-        r : float
-            radial position a < r < b
+        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)
+        f_a_cs_turn_steel : float
+            Steel area fraction of the central solenoid turn cross-section
+            (dimensionless)
 
         Returns
         -------
         float
-            hoop stress (MPa)
-        """
-        a = self.data.pf_coil.r_pf_coil_inner[self.data.pf_coil.n_cs_pf_coils - 1]
+            hoop stress at the specified radial location (MPa)
 
-        # Outer radius of central Solenoid [m]
-        b = self.data.pf_coil.r_pf_coil_outer[self.data.pf_coil.n_cs_pf_coils - 1]
+        References
+        ----------
+        - M. N. Wilson, Superconducting Magnets. Oxford University Press, USA, 1983.
+        ‌
+        """
+        alpha = r_cs_outer / r_cs_inner
 
         # alpha
-        alpha = b / a
+        alpha = r_cs_outer / r_cs_inner
 
         # epsilon
-        epsilon = r / a
-
-        # Field at inner radius of coil [T]
-        b_a = self.data.pf_coil.b_cs_peak_pulse_start
+        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 = self.data.pf_coil.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)))
+        hp_term_1 = k * ((2.0e0 + f_poisson_cs_structure) / (3.0e0 * (alpha + 1.0e0)))
 
         hp_term_2 = (
             alpha**2
@@ -3866,24 +4079,170 @@ def hoop_stress(self, r):
             + 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 / self.data.pf_coil.f_a_cs_turn_steel
+        return s_hoop_nom / f_a_cs_turn_steel
+
+    def calculate_cs_radial_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 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
+        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 same as for an infinite solenoid
+        b_cs_outer = 0.0e0
+
+        # K term
+        k = ((alpha * b_cs_inner - b_cs_outer) * j_cs * r_cs_inner) / (alpha - 1.0e0)
+
+        # M term
+        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)))
+
+        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
+
+    def plot_cs_radial_hoop_stress_profile(
+        self,
+        axis: plt.Axes,
+        mfile: 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,
+                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
+        ])
+
+        axis.plot(
+            radii,
+            stress_values / 1e6,
+            linewidth=2,
+            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("CS Hoop Stress at BOP")
+        axis.grid(True, alpha=0.3)
+
+    def plot_cs_radial_stress_profile(
+        self,
+        axis: plt.Axes,
+        mfile: 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,
+            linewidth=2,
+            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")
 
 
 def peak_b_field_at_pf_coil(
diff --git a/tests/unit/models/test_costs_1990.py b/tests/unit/models/test_costs_1990.py
index 164c7c4145..cbbcd339ae 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(
+        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 5ca0ba00ed..f24fe5ecdf 100644
--- a/tests/unit/models/test_pfcoil.py
+++ b/tests/unit/models/test_pfcoil.py
@@ -21,6 +21,8 @@
 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 (
+    CSEUDEMOTurnGeometry,
+    CSGeometry,
     PFCoil,
     calculate_b_field_at_point,
     fixb,
@@ -1970,82 +1972,34 @@ 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.
 
-    :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.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
 
@@ -2093,17 +2047,59 @@ 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,
+                0.5654866776461627,
+                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, 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,
+                3.1101767270538945,
+                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,
@@ -2270,27 +2266,22 @@ def test_calculate_cs_turn_geometry_eu_demo_basic(cs_coil):
     radius_cs_turn_corners = 0.01  # m
     f_a_cs_turn_steel = 0.4
 
-    (
-        dz_cs_turn,
-        dr_cs_turn,
-        radius_cs_turn_cable_space,
-        dr_cs_turn_conduit,
-        dz_cs_turn_conduit,
-    ) = cs_coil.calculate_cs_turn_geometry_eu_demo(
+    result = cs_coil.calculate_cs_turn_geometry_eu_demo(
         a_cs_turn, f_dr_dz_cs_turn, radius_cs_turn_corners, f_a_cs_turn_steel
     )
 
     # Check types and positivity
-    assert isinstance(dz_cs_turn, float)
-    assert isinstance(dr_cs_turn, float)
-    assert isinstance(radius_cs_turn_cable_space, float)
-    assert isinstance(dr_cs_turn_conduit, float)
-    assert isinstance(dz_cs_turn_conduit, float)
-    assert dz_cs_turn > 0
-    assert dr_cs_turn > 0
-    assert radius_cs_turn_cable_space > 0
-    assert dr_cs_turn_conduit > 0
-    assert dz_cs_turn_conduit > 0
+    assert isinstance(result, CSEUDEMOTurnGeometry)
+    assert isinstance(result.dz_cs_turn, float)
+    assert isinstance(result.dr_cs_turn, float)
+    assert isinstance(result.radius_cs_turn_cable_space, float)
+    assert isinstance(result.dr_cs_turn_conduit, float)
+    assert isinstance(result.dz_cs_turn_conduit, float)
+    assert result.dz_cs_turn > 0
+    assert result.dr_cs_turn > 0
+    assert result.radius_cs_turn_cable_space > 0
+    assert result.dr_cs_turn_conduit > 0
+    assert result.dz_cs_turn_conduit > 0
 
 
 @pytest.mark.parametrize(
@@ -2303,18 +2294,12 @@ def test_calculate_cs_turn_geometry_eu_demo_zero_corner(
     cs_coil, a_cs_turn, f_dr_dz_cs_turn, radius_cs_turn_corners, f_a_cs_turn_steel
 ):
     # Test with zero corner radius
-    (
-        dz_cs_turn,
-        dr_cs_turn,
-        radius_cs_turn_cable_space,
-        _dr_cs_turn_conduit,
-        _dz_cs_turn_conduit,
-    ) = cs_coil.calculate_cs_turn_geometry_eu_demo(
+    result = cs_coil.calculate_cs_turn_geometry_eu_demo(
         a_cs_turn, f_dr_dz_cs_turn, radius_cs_turn_corners, f_a_cs_turn_steel
     )
-    assert dz_cs_turn > 0
-    assert dr_cs_turn > 0
-    assert radius_cs_turn_cable_space > 0
+    assert result.dz_cs_turn > 0
+    assert result.dr_cs_turn > 0
+    assert result.radius_cs_turn_cable_space > 0
 
 
 @pytest.mark.parametrize(
@@ -2327,17 +2312,11 @@ def test_calculate_cs_turn_geometry_eu_demo_high_aspect(
     cs_coil, a_cs_turn, f_dr_dz_cs_turn, radius_cs_turn_corners, f_a_cs_turn_steel
 ):
     # High aspect ratio
-    (
-        dz_cs_turn,
-        dr_cs_turn,
-        radius_cs_turn_cable_space,
-        _dr_cs_turn_conduit,
-        _dz_cs_turn_conduit,
-    ) = cs_coil.calculate_cs_turn_geometry_eu_demo(
+    result = cs_coil.calculate_cs_turn_geometry_eu_demo(
         a_cs_turn, f_dr_dz_cs_turn, radius_cs_turn_corners, f_a_cs_turn_steel
     )
-    assert dr_cs_turn > dz_cs_turn
-    assert radius_cs_turn_cable_space > 0
+    assert result.dr_cs_turn > result.dz_cs_turn
+    assert result.radius_cs_turn_cable_space > 0
 
 
 @pytest.mark.parametrize(
@@ -2350,28 +2329,25 @@ def test_calculate_cs_turn_geometry_eu_demo_output_consistency(
     cs_coil, a_cs_turn, f_dr_dz_cs_turn, radius_cs_turn_corners, f_a_cs_turn_steel
 ):
     # Check that the sum of cable space and conduit thicknesses doesn't exceed half-width
-    (
-        dz_cs_turn,
-        _dr_cs_turn,
-        radius_cs_turn_cable_space,
-        _dr_cs_turn_conduit,
-        dz_cs_turn_conduit,
-    ) = cs_coil.calculate_cs_turn_geometry_eu_demo(
+    result = cs_coil.calculate_cs_turn_geometry_eu_demo(
         a_cs_turn, f_dr_dz_cs_turn, radius_cs_turn_corners, f_a_cs_turn_steel
     )
     # The cable space plus conduit thickness should not exceed half the turn width
-    assert radius_cs_turn_cable_space + dz_cs_turn_conduit <= dz_cs_turn / 2 + 1e-8
+    assert (
+        result.radius_cs_turn_cable_space + result.dz_cs_turn_conduit
+        <= result.dz_cs_turn / 2 + 1e-8
+    )
 
 
 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
@@ -2386,12 +2362,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
@@ -2401,12 +2377,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
@@ -2420,12 +2396,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
@@ -2611,7 +2587,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)