diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py
index cae753ed08..d69f96724d 100644
--- a/process/plasma_profiles.py
+++ b/process/plasma_profiles.py
@@ -6,7 +6,6 @@
 from process.fortran import (
     constants,
     divertor_variables,
-    maths_library,
     physics_variables,
 )
 
@@ -96,14 +95,13 @@ def parabolic_paramterisation(self):
         )
 
         #  Line averaged electron density (IPDG89)
-        #  0.5*gamfun(0.5) = 0.5*sqrt(pi) = 0.886227
 
         physics_variables.dnla = (
             physics_variables.dene
             * (1.0 + physics_variables.alphan)
             * 0.886227
-            * maths_library.gamfun(physics_variables.alphan + 1.0)
-            / maths_library.gamfun(physics_variables.alphan + 1.5e0)
+            * sp.special.gamma(physics_variables.alphan + 1.0)
+            / sp.special.gamma(physics_variables.alphan + 1.5)
         )
 
         #  Density-weighted temperatures
diff --git a/process/profiles.py b/process/profiles.py
index 58d2efbd63..37355877b2 100644
--- a/process/profiles.py
+++ b/process/profiles.py
@@ -1,9 +1,9 @@
 import numpy as np
 import logging
-from scipy import integrate
+import scipy as sp
 from abc import ABC, abstractmethod
 
-from process.fortran import maths_library, physics_variables, error_handling
+from process.fortran import physics_variables, error_handling
 
 logger = logging.getLogger(__name__)
 # Logging handler for console output
@@ -49,7 +49,7 @@ def integrate_profile_y(self):
         """
         Integrate profile_y values using scipy.integrate.simpson() function.
         """
-        self.profile_integ = integrate.simpson(
+        self.profile_integ = sp.integrate.simpson(
             self.profile_y, x=self.profile_x, dx=self.profile_dx
         )
 
@@ -261,20 +261,16 @@ def tcore(rhopedt, teped, tesep, tav, alphat, tbeta):
         :return: core temperature
         :rtype: numpy.array
         """
-        #  For integer values of alphat, the limit of
-        #  gamfun(-alphat)*sin(pi*alphat) needs to be calculated directly
 
         gamfac = (
-            maths_library.gamfun(1 + alphat + 2 / tbeta)
-            / maths_library.gamfun((2 + tbeta) / tbeta)
+            sp.special.gamma(1 + alphat + 2 / tbeta)
+            / sp.special.gamma((2 + tbeta) / tbeta)
             / rhopedt**2
         )
         if abs(alphat - np.around(alphat)) <= 1e-7:
-            gamfac = -gamfac / maths_library.gamfun(1 + alphat)
+            gamfac = -gamfac / sp.special.gamma(1 + alphat)
         else:
-            gamfac = (
-                gamfac * maths_library.gamfun(-alphat) * np.sin(np.pi * alphat) / np.pi
-            )
+            gamfac = gamfac * sp.special.gamma(-alphat) * np.sin(np.pi * alphat) / np.pi
 
         #  Calculate core temperature
 
diff --git a/tests/unit/test_plasma_profiles.py b/tests/unit/test_plasma_profiles.py
index c839b8fabd..e428d63562 100644
--- a/tests/unit/test_plasma_profiles.py
+++ b/tests/unit/test_plasma_profiles.py
@@ -86,15 +86,15 @@ class TProfileParam(NamedTuple):
             alphat=1.45,
             teped=5.5,
             expected_tprofile=[
-                18.84880504,
-                18.73828506,
-                18.40251429,
-                17.82801477,
-                16.98907408,
-                15.84091405,
-                14.3037138,
-                12.21867652,
-                9.17694539,
+                18.85,
+                18.739472621498333,
+                18.403679368327712,
+                17.829141392239833,
+                16.990144534125456,
+                15.841907634203302,
+                14.30460446612375,
+                12.219427594826554,
+                9.177492824141694,
                 1.0,
             ],
         ),
@@ -118,7 +118,7 @@ def test_tcore():
     tbeta = 2.0
 
     assert tprofile.tcore(rhopedt, tped, tsep, tav, alphat, tbeta) == pytest.approx(
-        28.089723663920328
+        28.09093632260765
     )
 
 
@@ -246,15 +246,15 @@ class PlasmaProfilesParam(NamedTuple):
             gradient_length_te=0.0,
             rminor=2.9264516129032256,
             expected_prn1=0.45623618297262686,
-            expected_ten=14.521871327399182,
-            expected_tin=14.521871327399182,
+            expected_ten=14.52233022043558,
+            expected_tin=14.52233022043558,
             expected_alphap=2.4500000000000002,
-            expected_te0=27.369013322953624,
-            expected_p0=868071.46874220832,
-            expected_pcoef=1.1110842637642833,
+            expected_te0=27.370104119511087,
+            expected_p0=868106.0658743214,
+            expected_pcoef=1.111119374172577,
             expected_ni0=9.210720071916929e19,
             expected_ne0=1.0585658890823703e20,
-            expected_ti0=27.369013322953624,
+            expected_ti0=27.370104119511087,
             expected_dnla=8.8687354645836431e19,
             expected_ti=13.07,
         ),
@@ -291,15 +291,15 @@ class PlasmaProfilesParam(NamedTuple):
             gradient_length_te=0.0,
             rminor=2.9264516129032256,
             expected_prn1=0.45623618297262686,
-            expected_ten=14.521871327399182,
-            expected_tin=14.521871327399182,
+            expected_ten=14.52233022043558,
+            expected_tin=14.52233022043558,
             expected_alphap=2.4500000000000002,
-            expected_te0=27.369013322953624,
-            expected_p0=868071.46874220832,
-            expected_pcoef=1.1110842637642833,
+            expected_te0=27.370104119511087,
+            expected_p0=868106.0658743214,
+            expected_pcoef=1.111119374172577,
             expected_ni0=9.210720071916929e19,
             expected_ne0=1.0585658890823703e20,
-            expected_ti0=27.369013322953624,
+            expected_ti0=27.370104119511087,
             expected_dnla=8.8687354645836431e19,
             expected_ti=13.07,
         ),