diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index ef794eada1..547377a5a4 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -244,6 +244,9 @@
 beta_toroidal_vol_avg: float = None
 """Plasma volume averaged toroidal beta"""
 
+beta_thermal_toroidal_profile: list[float] = None
+"""toroidal beta profile"""
+
 
 beta_thermal_vol_avg: float = None
 """Plasma volume averaged thermal beta"""
@@ -288,6 +291,15 @@
 b_plasma_toroidal_on_axis: float = None
 """Plasma toroidal field on axis (T) (`iteration variable 2`)"""
 
+b_plasma_inboard_toroidal: float = None
+"""Plasma inboard toroidal field (T)"""
+
+b_plasma_outboard_toroidal: float = None
+"""Plasma outboard toroidal field (T)"""
+
+b_plasma_toroidal_profile: list[float] = None
+"""toroidal field profile in plasma (T)"""
+
 
 b_plasma_total: float = None
 """Sum of plasma total toroidal + poloidal field (T)"""
@@ -809,7 +821,7 @@
 pres_plasma_thermal_on_axis: float = None
 """Plasma central thermal pressure (no fast ions or beam pressure) (Pa)"""
 
-pres_plasma_total_profile: list[float] = None
+pres_plasma_thermal_total_profile: list[float] = None
 """Profile of total pressure in plasma (Pa)"""
 
 pres_plasma_electron_profile: list[float] = None
@@ -1383,6 +1395,7 @@ def init_physics_variables():
     global beta_poloidal_vol_avg
     global beta_poloidal_eps
     global beta_toroidal_vol_avg
+    global beta_thermal_toroidal_profile
     global beta_thermal_vol_avg
     global beta_thermal_poloidal_vol_avg
     global beta_thermal_toroidal_vol_avg
@@ -1394,6 +1407,9 @@ def init_physics_variables():
     global betbm0
     global b_plasma_poloidal_average
     global b_plasma_toroidal_on_axis
+    global b_plasma_toroidal_inboard
+    global b_plasma_toroidal_outboard
+    global b_plasma_toroidal_profile
     global b_plasma_total
     global e_plasma_magnetic_stored
     global burnup
@@ -1493,7 +1509,7 @@ def init_physics_variables():
     global nd_plasma_ions_on_axis
     global m_s_limit
     global pres_plasma_thermal_on_axis
-    global pres_plasma_total_profile
+    global pres_plasma_thermal_total_profile
     global pres_plasma_electron_profile
     global pres_plasma_ion_total_profile
     global pres_plasma_fuel_profile
@@ -1636,6 +1652,7 @@ def init_physics_variables():
     beta_poloidal_vol_avg = 0.0
     beta_poloidal_eps = 0.0
     beta_toroidal_vol_avg = 0.0
+    beta_thermal_toroidal_profile = []
     beta_thermal_vol_avg = 0.0
     beta_thermal_poloidal_vol_avg = 0.0
     beta_thermal_toroidal_vol_avg = 0.0
@@ -1647,6 +1664,9 @@ def init_physics_variables():
     betbm0 = 1.5
     b_plasma_poloidal_average = 0.0
     b_plasma_toroidal_on_axis = 5.68
+    b_plasma_toroidal_inboard = 0.0
+    b_plasma_toroidal_outboard = 0.0
+    b_plasma_toroidal_profile = []
     b_plasma_total = 0.0
     e_plasma_magnetic_stored = 0.0
     burnup = 0.0
@@ -1746,7 +1766,7 @@ def init_physics_variables():
     nd_plasma_ions_on_axis = 0.0
     m_s_limit = 0.3
     pres_plasma_thermal_on_axis = 0.0
-    pres_plasma_total_profile = []
+    pres_plasma_thermal_total_profile = []
     pres_plasma_electron_profile = []
     pres_plasma_ion_total_profile = []
     pres_plasma_fuel_profile = []
diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index 54a18014df..f31113121e 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -10840,8 +10840,8 @@ def plot_plasma_pressure_profiles(axis, mfile_data, scan):
         mfile_data.data[f"pres_plasma_ion_total_profile{i}"].get_scan(scan)
         for i in range(n_plasma_profile_elements)
     ]
-    pres_plasma_total_profile = [
-        mfile_data.data[f"pres_plasma_total_profile{i}"].get_scan(scan)
+    pres_plasma_thermal_total_profile = [
+        mfile_data.data[f"pres_plasma_thermal_total_profile{i}"].get_scan(scan)
         for i in range(n_plasma_profile_elements)
     ]
     pres_plasma_profile_fuel = [
@@ -10851,7 +10851,9 @@ def plot_plasma_pressure_profiles(axis, mfile_data, scan):
     pres_plasma_profile_kpa = [p / 1000.0 for p in pres_plasma_profile]
     pres_plasma_profile_ion_kpa = [p / 1000.0 for p in pres_plasma_profile_ion]
     pres_plasma_profile_fuel_kpa = [p / 1000.0 for p in pres_plasma_profile_fuel]
-    pres_plasma_profile_total_kpa = [p / 1000.0 for p in pres_plasma_total_profile]
+    pres_plasma_profile_total_kpa = [
+        p / 1000.0 for p in pres_plasma_thermal_total_profile
+    ]
 
     axis.plot(
         np.linspace(0, 1, len(pres_plasma_profile_kpa)),
@@ -10903,7 +10905,7 @@ def plot_plasma_pressure_gradient_profiles(axis, mfile_data, scan):
         for i in range(n_plasma_profile_elements)
     ]
     pres_plasma_profile_total = [
-        mfile_data.data[f"pres_plasma_total_profile{i}"].get_scan(scan)
+        mfile_data.data[f"pres_plasma_thermal_total_profile{i}"].get_scan(scan)
         for i in range(n_plasma_profile_elements)
     ]
     pres_plasma_profile_fuel = [
@@ -11500,6 +11502,152 @@ def plot_fusion_rate_contours(
     # ===============================================
 
 
+def plot_magnetic_fields_in_plasma(axis, mfile_data, scan):
+    # Plot magnetic field profiles inside the plasma boundary
+
+    n_plasma_profile_elements = int(
+        mfile_data.data["n_plasma_profile_elements"].get_scan(scan)
+    )
+
+    # Get toroidal magnetic field profile (in Tesla)
+    b_plasma_toroidal_profile = [
+        mfile_data.data[f"b_plasma_toroidal_profile{i}"].get_scan(scan)
+        for i in range(2 * n_plasma_profile_elements)
+    ]
+
+    # Get major and minor radius for x-axis in metres
+    rmajor = mfile_data.data["rmajor"].get_scan(scan)
+    rminor = mfile_data.data["rminor"].get_scan(scan)
+
+    # Plot magnetic field first (background)
+    axis.plot(
+        np.linspace(rmajor - rminor, rmajor + rminor, len(b_plasma_toroidal_profile)),
+        b_plasma_toroidal_profile,
+        color="blue",
+        label="Toroidal B-field [T]",
+        linewidth=2,
+    )
+
+    # Plot plasma on top of magnetic field, displaced vertically by bt
+    plot_plasma(axis, mfile_data, scan, colour_scheme=1)
+
+    # Plot plasma centre dot
+    axis.plot(rmajor, 0, marker="o", color="red", markersize=8, label="Plasma Centre")
+
+    # Plot vertical lines at plasma edge
+    axis.axvline(
+        rmajor - rminor,
+        color="green",
+        linestyle="--",
+        linewidth=1,
+    )
+    axis.axvline(rmajor + rminor, color="green", linestyle="--", linewidth=1)
+
+    # Plot horizontal line for toroidal magnetic field at plasma inboard
+    axis.axhline(
+        mfile_data.data[f"b_plasma_toroidal_profile{0}"].get_scan(scan),
+        color="blue",
+        linestyle="--",
+        linewidth=1.0,
+    )
+
+    # Plot horizontal line for toroidal magnetic field at plasma centre
+    axis.axhline(
+        mfile_data.data["b_plasma_toroidal_on_axis"].get_scan(scan),
+        color="blue",
+        linestyle="--",
+        linewidth=1.0,
+    )
+
+    # Plot horizontal line for toroidal magnetic field at plasma outboard
+    axis.axhline(
+        b_plasma_toroidal_profile[-1],
+        color="blue",
+        linestyle="--",
+        linewidth=1.0,
+    )
+
+    # Text box for inboard toroidal field
+    axis.text(
+        0.1,
+        0.025,
+        f"$B_{{\\text{{T,inboard}}}}={mfile_data.data['b_plasma_inboard_toroidal'].get_scan(scan):.2f}$ T",
+        verticalalignment="center",
+        horizontalalignment="center",
+        transform=axis.transAxes,
+        bbox={
+            "boxstyle": "round",
+            "facecolor": "wheat",
+            "alpha": 1.0,
+            "linewidth": 2,
+        },
+    )
+
+    # Text box for outboard toroidal field
+    axis.text(
+        0.9,
+        0.1,
+        f"$B_{{\\text{{T,outboard}}}}={mfile_data.data['b_plasma_outboard_toroidal'].get_scan(scan):.2f}$ T",
+        verticalalignment="center",
+        horizontalalignment="center",
+        transform=axis.transAxes,
+        bbox={
+            "boxstyle": "round",
+            "facecolor": "wheat",
+            "alpha": 1.0,
+            "linewidth": 2,
+        },
+    )
+
+    axis.set_xlabel("Radial Position [m]")
+    axis.set_ylabel("Toroidal Magnetic Field [T]")
+    axis.set_title("Toroidal Magnetic Field Profile in Plasma")
+    axis.minorticks_on()
+    # Enable grid for both major and minor ticks
+    axis.grid(which="both", linestyle="--", alpha=0.5)
+    axis.grid(which="minor", linestyle=":", alpha=0.3)
+    axis.legend(loc="lower right")
+    axis.set_xlim(rmajor - 1.25 * rminor, rmajor + 1.25 * rminor)
+
+
+def plot_beta_profiles(axis, mfile_data, scan):
+    # Plot the beta profiles on the given axis
+
+    n_plasma_profile_elements = int(
+        mfile_data.data["n_plasma_profile_elements"].get_scan(scan)
+    )
+
+    beta_plasma_toroidal_profile = [
+        mfile_data.data[f"beta_thermal_toroidal_profile{i}"].get_scan(scan)
+        for i in range(2 * n_plasma_profile_elements)
+    ]
+
+    axis.plot(
+        np.linspace(-1, 1, 2 * n_plasma_profile_elements),
+        beta_plasma_toroidal_profile,
+        color="blue",
+        label="$\\beta_t$",
+    )
+
+    axis.axhline(
+        mfile_data.data["beta_thermal_toroidal_vol_avg"].get_scan(scan),
+        color="blue",
+        linestyle="--",
+        linewidth=1.0,
+        label="$\\langle \\beta_t \\rangle_{\\text{V}}$",
+    )
+
+    axis.set_xlabel("$\\rho$ [r/a]")
+    axis.set_ylabel("$\\beta$")
+    axis.minorticks_on()
+    axis.grid(which="minor", linestyle=":", linewidth=0.5, alpha=0.5)
+    axis.set_title("Thermal Beta Profiles")
+    axis.legend()
+    axis.axvline(x=0, color="black", linestyle="--", linewidth=1)
+    axis.grid(True, linestyle="--", alpha=0.5)
+    axis.set_ylim(bottom=0.0)
+
+
 def main_plot(
     fig0,
     fig1,
@@ -11524,6 +11672,7 @@ def main_plot(
     fig20,
     fig21,
     fig22,
+    fig23,
     m_file_data,
     scan,
     imp="../data/lz_non_corona_14_elements/",
@@ -11660,9 +11809,12 @@ def main_plot(
         )
         ax.axis("off")
 
+    plot_magnetic_fields_in_plasma(fig9.add_subplot(122), m_file_data, scan)
+    plot_beta_profiles(fig9.add_subplot(221), m_file_data, scan)
+
     # Plot poloidal cross-section
     poloidal_cross_section(
-        fig9.add_subplot(121, aspect="equal"),
+        fig10.add_subplot(121, aspect="equal"),
         m_file_data,
         scan,
         demo_ranges,
@@ -11671,7 +11823,7 @@ def main_plot(
 
     # Plot toroidal cross-section
     toroidal_cross_section(
-        fig9.add_subplot(122, aspect="equal"),
+        fig10.add_subplot(122, aspect="equal"),
         m_file_data,
         scan,
         demo_ranges,
@@ -11679,19 +11831,19 @@ def main_plot(
     )
 
     # Plot color key
-    ax17 = fig9.add_subplot(222)
+    ax17 = fig10.add_subplot(222)
     ax17.set_position([0.5, 0.5, 0.5, 0.5])
     color_key(ax17, m_file_data, scan, colour_scheme)
 
-    ax18 = fig10.add_subplot(211)
+    ax18 = fig11.add_subplot(211)
     ax18.set_position([0.1, 0.33, 0.8, 0.6])
     plot_radial_build(ax18, m_file_data, colour_scheme)
 
     # Make each axes smaller vertically to leave room for the legend
-    ax185 = fig11.add_subplot(211)
+    ax185 = fig12.add_subplot(211)
     ax185.set_position([0.1, 0.61, 0.8, 0.32])
 
-    ax18b = fig11.add_subplot(212)
+    ax18b = fig12.add_subplot(212)
     ax18b.set_position([0.1, 0.13, 0.8, 0.32])
     plot_upper_vertical_build(ax185, m_file_data, colour_scheme)
     plot_lower_vertical_build(ax18b, m_file_data, colour_scheme)
@@ -11699,57 +11851,57 @@ def main_plot(
     # Can only plot WP and turn structure if superconducting coil at the moment
     if m_file_data.data["i_tf_sup"].get_scan(scan) == 1:
         # TF coil with WP
-        ax19 = fig12.add_subplot(221, aspect="equal")
+        ax19 = fig13.add_subplot(221, aspect="equal")
         ax19.set_position([
             0.025,
             0.45,
             0.5,
             0.5,
         ])  # Half height, a bit wider, top left
-        plot_superconducting_tf_wp(ax19, m_file_data, scan, fig12)
+        plot_superconducting_tf_wp(ax19, m_file_data, scan, fig13)
 
         # TF coil turn structure
-        ax20 = fig13.add_subplot(325, aspect="equal")
+        ax20 = fig14.add_subplot(325, aspect="equal")
         ax20.set_position([0.025, 0.5, 0.4, 0.4])
-        plot_tf_cable_in_conduit_turn(ax20, fig13, m_file_data, scan)
-        plot_205 = fig13.add_subplot(223, aspect="equal")
+        plot_tf_cable_in_conduit_turn(ax20, fig14, m_file_data, scan)
+        plot_205 = fig14.add_subplot(223, aspect="equal")
         plot_205.set_position([0.075, 0.1, 0.3, 0.3])
-        plot_cable_in_conduit_cable(plot_205, fig16, m_file_data, scan)
+        plot_cable_in_conduit_cable(plot_205, fig14, m_file_data, scan)
     else:
-        ax19 = fig12.add_subplot(211, aspect="equal")
+        ax19 = fig13.add_subplot(211, aspect="equal")
         ax19.set_position([0.06, 0.55, 0.675, 0.4])
-        plot_resistive_tf_wp(ax19, m_file_data, scan, fig12)
+        plot_resistive_tf_wp(ax19, m_file_data, scan, fig13)
 
     plot_tf_coil_structure(
-        fig14.add_subplot(111, aspect="equal"), m_file_data, scan, colour_scheme
+        fig15.add_subplot(111, aspect="equal"), m_file_data, scan, colour_scheme
     )
 
-    axes = fig15.subplots(nrows=3, ncols=1, sharex=True).flatten()
+    axes = fig16.subplots(nrows=3, ncols=1, sharex=True).flatten()
     plot_tf_stress(axes)
 
-    plot_bootstrap_comparison(fig16.add_subplot(221), m_file_data, scan)
-    plot_h_threshold_comparison(fig16.add_subplot(224), m_file_data, scan)
-    plot_density_limit_comparison(fig17.add_subplot(221), m_file_data, scan)
-    plot_confinement_time_comparison(fig17.add_subplot(224), m_file_data, scan)
-    plot_current_profiles_over_time(fig18.add_subplot(111), m_file_data, scan)
+    plot_bootstrap_comparison(fig17.add_subplot(221), m_file_data, scan)
+    plot_h_threshold_comparison(fig17.add_subplot(224), m_file_data, scan)
+    plot_density_limit_comparison(fig18.add_subplot(221), m_file_data, scan)
+    plot_confinement_time_comparison(fig18.add_subplot(224), m_file_data, scan)
+    plot_current_profiles_over_time(fig19.add_subplot(111), m_file_data, scan)
 
     plot_cs_coil_structure(
-        fig19.add_subplot(121, aspect="equal"), fig19, m_file_data, scan
+        fig20.add_subplot(121, aspect="equal"), fig20, m_file_data, scan
     )
     plot_cs_turn_structure(
-        fig19.add_subplot(224, aspect="equal"), fig19, m_file_data, scan
+        fig20.add_subplot(224, aspect="equal"), fig20, m_file_data, scan
     )
 
     plot_first_wall_top_down_cross_section(
-        fig20.add_subplot(221, aspect="equal"), m_file_data, scan
+        fig21.add_subplot(221, aspect="equal"), m_file_data, scan
     )
-    plot_first_wall_poloidal_cross_section(fig20.add_subplot(122), m_file_data, scan)
-    plot_fw_90_deg_pipe_bend(fig20.add_subplot(337), m_file_data, scan)
+    plot_first_wall_poloidal_cross_section(fig21.add_subplot(122), m_file_data, scan)
+    plot_fw_90_deg_pipe_bend(fig21.add_subplot(337), m_file_data, scan)
 
-    plot_blkt_pipe_bends(fig21, m_file_data, scan)
+    plot_blkt_pipe_bends(fig22, m_file_data, scan)
 
     plot_main_power_flow(
-        fig22.add_subplot(111, aspect="equal"), m_file_data, scan, fig22
+        fig23.add_subplot(111, aspect="equal"), m_file_data, scan, fig23
     )
 
 
@@ -12066,6 +12218,7 @@ def main(args=None):
     page20 = plt.figure(figsize=(12, 9), dpi=80)
     page21 = plt.figure(figsize=(12, 9), dpi=80)
     page22 = plt.figure(figsize=(12, 9), dpi=80)
+    page23 = plt.figure(figsize=(12, 9), dpi=80)
 
     # run main_plot
     main_plot(
@@ -12092,6 +12245,7 @@ def main(args=None):
         page20,
         page21,
         page22,
+        page23,
         m_file,
         scan=scan,
         demo_ranges=demo_ranges,
@@ -12123,6 +12277,7 @@ def main(args=None):
         pdf.savefig(page20)
         pdf.savefig(page21)
         pdf.savefig(page22)
+        pdf.savefig(page23)
 
     # show fig if option used
     if args.show:
@@ -12151,6 +12306,7 @@ def main(args=None):
     plt.close(page20)
     plt.close(page21)
     plt.close(page22)
+    plt.close(page23)
 
 
 if __name__ == "__main__":
diff --git a/process/physics.py b/process/physics.py
index 736ae760e0..d2b2a8caa3 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -1735,6 +1735,34 @@ def physics(self):
             + physics_variables.b_plasma_poloidal_average**2
         )
 
+        # Calculate the toroidal field across the plasma
+        # Calculate the toroidal field profile across the plasma (1/R dependence)
+        # Double element size to include both sides of the plasma
+        rho = np.linspace(
+            physics_variables.rmajor - physics_variables.rminor,
+            physics_variables.rmajor + physics_variables.rminor,
+            2 * physics_variables.n_plasma_profile_elements,
+        )
+
+        # Calculate the inboard and outboard toroidal field
+        physics_variables.b_plasma_inboard_toroidal = (
+            physics_variables.rmajor
+            * physics_variables.b_plasma_toroidal_on_axis
+            / (physics_variables.rmajor - physics_variables.rminor)
+        )
+
+        physics_variables.b_plasma_outboard_toroidal = (
+            physics_variables.rmajor
+            * physics_variables.b_plasma_toroidal_on_axis
+            / (physics_variables.rmajor + physics_variables.rminor)
+        )
+
+        # Avoid division by zero at the magnetic axis
+        rho = np.where(rho == 0, 1e-10, rho)
+        physics_variables.b_plasma_toroidal_profile = (
+            physics_variables.rmajor * physics_variables.b_plasma_toroidal_on_axis / rho
+        )
+
         # ============================================
 
         # -----------------------------------------------------
@@ -1780,6 +1808,25 @@ def physics(self):
             )
             ** 2
         )
+
+        # =======================================================
+
+        # Mirror the pressure profiles to match the doubled toroidal field profile
+        pres_profile_total = np.concatenate([
+            physics_variables.pres_plasma_thermal_total_profile[::-1],
+            physics_variables.pres_plasma_thermal_total_profile,
+        ])
+
+        physics_variables.beta_thermal_toroidal_profile = np.array([
+            self.calculate_plasma_beta(
+                pres_plasma=pres_profile_total[i],
+                b_field=physics_variables.b_plasma_toroidal_profile[i],
+            )
+            for i in range(len(physics_variables.b_plasma_toroidal_profile))
+        ])
+
+        # =======================================================
+
         physics_variables.beta_norm_thermal = (
             1.0e8
             * physics_variables.beta_thermal_vol_avg
@@ -3853,6 +3900,22 @@ def calculate_plasma_current(
 
         return b_plasma_poloidal_average, qstar, plasma_current
 
+    def calculate_plasma_beta(self, pres_plasma: float, b_field: float) -> float:
+        """
+        Calculate the plasma beta for a given pressure and field.
+
+        :param pres_plasma: Plasma pressure (in Pascals).
+        :type pres_plasma: float
+        :param b_field: Magnetic field strength (in Tesla).
+        :type b_field: float
+        :returns: The plasma beta (dimensionless).
+        :rtype: float
+
+        Plasma beta is the ratio of plasma pressure to magnetic pressure.
+        """
+
+        return 2 * constants.RMU0 * pres_plasma / (b_field**2)
+
     def outtim(self):
         po.oheadr(self.outfile, "Times")
 
@@ -4238,6 +4301,26 @@ def outplas(self):
                 "(b_plasma_toroidal_on_axis)",
                 physics_variables.b_plasma_toroidal_on_axis,
             )
+            po.ovarrf(
+                self.outfile,
+                "Toroidal field at plasma inboard (T)",
+                "(b_plasma_inboard_toroidal)",
+                physics_variables.b_plasma_inboard_toroidal,
+            )
+            po.ovarrf(
+                self.outfile,
+                "Toroidal field at plasma outboard (T)",
+                "(b_plasma_outboard_toroidal)",
+                physics_variables.b_plasma_outboard_toroidal,
+            )
+
+            for i in range(len(physics_variables.b_plasma_toroidal_profile)):
+                po.ovarre(
+                    self.mfile,
+                    f"Toroidal field in plasma at point {i}",
+                    f"b_plasma_toroidal_profile{i}",
+                    physics_variables.b_plasma_toroidal_profile[i],
+                )
             po.ovarrf(
                 self.outfile,
                 "Average poloidal field (T)",
@@ -4401,6 +4484,14 @@ def outplas(self):
             physics_variables.beta_toroidal_vol_avg,
             "OP ",
         )
+        for i in range(len(physics_variables.beta_thermal_toroidal_profile)):
+            po.ovarre(
+                self.mfile,
+                f"Beta toroidal profile at point {i}",
+                f"beta_thermal_toroidal_profile{i}",
+                physics_variables.beta_thermal_toroidal_profile[i],
+            )
+
         po.ovarre(
             self.outfile,
             "Fast alpha beta",
@@ -4654,12 +4745,12 @@ def outplas(self):
         )
         # As array output is not currently supported, each element is output as a float instance
         # Output plasma pressure profiles to mfile
-        for i in range(len(physics_variables.pres_plasma_total_profile)):
+        for i in range(len(physics_variables.pres_plasma_thermal_total_profile)):
             po.ovarre(
                 self.mfile,
                 f"Total plasma pressure at point {i}",
-                f"(pres_plasma_total_profile{i})",
-                physics_variables.pres_plasma_total_profile[i],
+                f"(pres_plasma_thermal_total_profile{i})",
+                physics_variables.pres_plasma_thermal_total_profile[i],
             )
         for i in range(len(physics_variables.pres_plasma_electron_profile)):
             po.ovarre(
diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py
index 2953afd120..d9860fa8ce 100644
--- a/process/plasma_profiles.py
+++ b/process/plasma_profiles.py
@@ -299,7 +299,7 @@ def calculate_profile_factors(self) -> None:
         )
 
         # Total pressure profile (Pa)
-        physics_variables.pres_plasma_total_profile = (
+        physics_variables.pres_plasma_thermal_total_profile = (
             physics_variables.pres_plasma_electron_profile
             + physics_variables.pres_plasma_ion_total_profile
         )