diff --git a/process/physics.py b/process/physics.py
index e616abff60..0aa4345617 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -24,7 +24,6 @@
     numerics,
     stellarator_variables,
     process_output as po,
-    profiles_module,
 )
 
 
@@ -987,7 +986,8 @@ def physics(self):
         current_drive_variables.pscf_scene = ps_fraction_scene(physics_variables.beta)
 
         current_drive_variables.bscf_sauter = (
-            current_drive_variables.cboot * self.bootstrap_fraction_sauter()
+            current_drive_variables.cboot
+            * self.bootstrap_fraction_sauter(self.plasma_profile)
         )
 
         current_drive_variables.bscf_sakai = (
@@ -4553,7 +4553,7 @@ def bootstrap_fraction_nevins(
         return 1.0e6 * aibs / plascur
 
     @staticmethod
-    def bootstrap_fraction_sauter():
+    def bootstrap_fraction_sauter(plasma_profile):
         """Bootstrap current fraction from Sauter et al scaling
         author: P J Knight, CCFE, Culham Science Centre
         None
@@ -4566,38 +4566,20 @@ def bootstrap_fraction_sauter():
         O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA)
         Physics of Plasmas <B>9</B> (2002) 5140
         """
-        NR = 200
-
-        roa = np.arange(1, NR + 1, step=1) / NR
+        NR = plasma_profile.profile_size
 
+        roa = plasma_profile.neprofile.profile_x
         rho = np.sqrt(physics_variables.xarea / np.pi) * roa
+
         sqeps = np.sqrt(roa * (physics_variables.rminor / physics_variables.rmajor))
 
-        ne = 1e-19 * np.vectorize(
-            lambda r: profiles_module.nprofile(
-                r,
-                physics_variables.rhopedn,
-                physics_variables.ne0,
-                physics_variables.neped,
-                physics_variables.nesep,
-                physics_variables.alphan,
-            )
-        )(roa)
+        ne = plasma_profile.neprofile.profile_y * 1e-19
         ni = (physics_variables.dnitot / physics_variables.dene) * ne
-        tempe = np.vectorize(
-            lambda r: profiles_module.tprofile(
-                r,
-                physics_variables.rhopedt,
-                physics_variables.te0,
-                physics_variables.teped,
-                physics_variables.tesep,
-                physics_variables.alphat,
-                physics_variables.tbeta,
-            )
-        )(roa)
-        tempi = (physics_variables.ti / physics_variables.te) * tempe
 
-        zef = np.full_like(tempi, physics_variables.zeff)  # Flat Zeff profile assumed
+        tempe = plasma_profile.teprofile.profile_y
+        tempi = (physics_variables.ti / physics_variables.te) * tempe
+        # Flat Zeff profile assumed
+        zef = np.full_like(tempi, physics_variables.zeff)
 
         # mu = 1/safety factor
         # Parabolic q profile assumed
@@ -4620,13 +4602,11 @@ def bootstrap_fraction_sauter():
         # Looping from 2 because dcsa etc should return 0 @ j == 1
         nr_rng = np.arange(2, NR)
         nr_rng_1 = nr_rng - 1
-
         drho = rho[nr_rng] - rho[nr_rng_1]
         da = 2 * np.pi * rho[nr_rng_1] * drho  # area of annulus
         dlogte_drho = (np.log(tempe[nr_rng]) - np.log(tempe[nr_rng_1])) / drho
         dlogti_drho = (np.log(tempi[nr_rng]) - np.log(tempi[nr_rng_1])) / drho
         dlogne_drho = (np.log(ne[nr_rng]) - np.log(ne[nr_rng_1])) / drho
-
         jboot = (
             0.5
             * (
diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py
index 32bbac8ce6..dd95d87dae 100644
--- a/tests/unit/test_physics.py
+++ b/tests/unit/test_physics.py
@@ -350,33 +350,33 @@ class BootstrapFractionSauterParam(NamedTuple):
     "bootstrapfractionsauterparam",
     (
         BootstrapFractionSauterParam(
-            dnitot=6.6125550702454276e19,
+            dnitot=7.1297522422781575e19,
             rminor=2.6666666666666665,
             tesep=0.10000000000000001,
-            ti=12,
+            ti=12.570861186498382,
             triang=0.5,
             q0=1,
             afuel=2.5,
-            zeff=2.0909945616489103,
-            rhopedn=0.94000000000000006,
-            bt=5.7000000000000002,
-            plascur=18398455.678867526,
+            zeff=2.5211399464385624,
+            rhopedn=0.9400000000000001,
+            bt=5.326133750416047,
+            plascur=16528278.760008096,
             xarea=38.39822223637151,
             fhe3=0,
             teped=5.5,
-            dene=7.5e19,
-            te=12,
+            dene=8.016748468651018e19,
+            te=12.570861186498382,
             rmajor=8,
             q=3.5,
-            nesep=4.1177885154594193e19,
-            te0=24.402321098330372,
-            neped=7.000240476281013e19,
+            nesep=3.6992211545476006e19,
+            te0=25.986118047669795,
+            neped=6.2886759627309195e19,
             tbeta=2,
-            ne0=8.515060981068918e19,
+            ne0=1.054474759840606e20,
             alphan=1,
-            rhopedt=0.94000000000000006,
+            rhopedt=0.9400000000000001,
             alphat=1.45,
-            expected_bfs=0.27635918746616817,
+            expected_bfs=0.4110838247346975,
         ),
     ),
 )
@@ -462,8 +462,8 @@ def test_bootstrap_fraction_sauter(bootstrapfractionsauterparam, monkeypatch, ph
     monkeypatch.setattr(
         physics_variables, "alphat", bootstrapfractionsauterparam.alphat
     )
-
-    bfs = physics.bootstrap_fraction_sauter()
+    physics.plasma_profile.run()
+    bfs = physics.bootstrap_fraction_sauter(physics.plasma_profile)
 
     assert bfs == pytest.approx(bootstrapfractionsauterparam.expected_bfs)