diff --git a/src/openpile/construct.py b/src/openpile/construct.py
index add93a4..7d2e8aa 100644
--- a/src/openpile/construct.py
+++ b/src/openpile/construct.py
@@ -646,11 +646,11 @@ def get_soil_properties(self) -> pd.DataFrame:
 
     def _postinit(self):
         def check_springs(arr):
-                check_nan = np.isnan(arr).any()
-                check_negative = (arr < 0).any()
-                 
-                return check_nan or check_negative
-        
+            check_nan = np.isnan(arr).any()
+            check_negative = (arr < 0).any()
+
+            return check_nan or check_negative
+
         def get_coordinates() -> pd.DataFrame:
             # Primary discretisation over x-axis
             x = np.array([], dtype=np.float16)
@@ -766,14 +766,14 @@ def create_springs() -> np.ndarray:
                 # py curve
                 if layer.lateral_model is None:
                     pass
-                else: 
+                else:
                     # Set local layer parameters for each element of the layer
                     for i in elements_for_layer:
                         # vertical effective stress
                         sig_v = self.soil_properties[
                             ["sigma_v top [kPa]", "sigma_v bottom [kPa]"]
                         ].iloc[i]
-                        # elevation 
+                        # elevation
                         elevation = self.soil_properties[["x_top [m]", "x_bottom [m]"]].iloc[i]
                         # depth from ground
                         depth_from_ground = (
@@ -783,7 +783,9 @@ def create_springs() -> np.ndarray:
                         pile_width = self.element_properties["Diameter [m]"].iloc[i]
 
                         # p-y curves
-                        if layer.lateral_model.spring_signature[0] and self.distributed_lateral:  # True if py spring function exist    
+                        if (
+                            layer.lateral_model.spring_signature[0] and self.distributed_lateral
+                        ):  # True if py spring function exist
 
                             # calculate springs (top and) for each element
                             for j in [0, 1]:
@@ -798,7 +800,9 @@ def create_springs() -> np.ndarray:
                                     output_length=spring_dim,
                                 )
 
-                        if layer.lateral_model.spring_signature[2] and self.distributed_moment:  # True if mt spring function exist    
+                        if (
+                            layer.lateral_model.spring_signature[2] and self.distributed_moment
+                        ):  # True if mt spring function exist
 
                             # calculate springs (top and) for each element
                             for j in [0, 1]:
@@ -813,60 +817,77 @@ def create_springs() -> np.ndarray:
                                     output_length=spring_dim,
                                 )
 
-                        #check if pile tip is within layer
-                        if layer.top >= self.pile.bottom_elevation and layer.bottom <= self.pile.bottom_elevation:
+                        # check if pile tip is within layer
+                        if (
+                            layer.top >= self.pile.bottom_elevation
+                            and layer.bottom <= self.pile.bottom_elevation
+                        ):
 
-                            #Hb curve
+                            # Hb curve
                             sig_v_tip = self.soil_properties["sigma_v bottom [kPa]"].iloc[-1]
 
-                            if layer.lateral_model.spring_signature[1] and self.base_shear:    
+                            if layer.lateral_model.spring_signature[1] and self.base_shear:
 
                                 # calculate Hb spring
                                 (Hb[0, 0, 0], Hb[0, 0, 1]) = layer.lateral_model.Hb_spring_fct(
                                     sig=sig_v_tip,
                                     X=self.pile.bottom_elevation,
                                     layer_height=(layer.top - layer.bottom),
-                                    depth_from_top_of_layer=(layer.top - self.pile.bottom_elevation),
+                                    depth_from_top_of_layer=(
+                                        layer.top - self.pile.bottom_elevation
+                                    ),
                                     D=pile_width,
                                     L=(self.soil.top_elevation - self.pile.bottom_elevation),
-                                    below_water_table=self.pile.bottom_elevation <= self.soil.water_elevation,
+                                    below_water_table=self.pile.bottom_elevation
+                                    <= self.soil.water_elevation,
                                     output_length=spring_dim,
                                 )
-                            
-                            #Mb curve
-                            if layer.lateral_model.spring_signature[3] and self.base_moment:    
+
+                            # Mb curve
+                            if layer.lateral_model.spring_signature[3] and self.base_moment:
 
                                 (Mb[0, 0, 0], Mb[0, 0, 1]) = layer.lateral_model.Mb_spring_fct(
                                     sig=sig_v_tip,
                                     X=self.pile.bottom_elevation,
                                     layer_height=(layer.top - layer.bottom),
-                                    depth_from_top_of_layer=(layer.top - self.pile.bottom_elevation),
+                                    depth_from_top_of_layer=(
+                                        layer.top - self.pile.bottom_elevation
+                                    ),
                                     D=pile_width,
                                     L=(self.soil.top_elevation - self.pile.bottom_elevation),
-                                    below_water_table=self.pile.bottom_elevation <= self.soil.water_elevation,
+                                    below_water_table=self.pile.bottom_elevation
+                                    <= self.soil.water_elevation,
                                     output_length=spring_dim,
                                 )
 
             if check_springs(py):
-                print('py springs have negative or NaN values.')
-                print("""if using PISA type springs, this can be due to parameters behind out of the parameter space.
+                print("py springs have negative or NaN values.")
+                print(
+                    """if using PISA type springs, this can be due to parameters behind out of the parameter space.
                 Please check that: 2 < L/D < 6.
-                """)
+                """
+                )
             if check_springs(mt):
-                print('mt springs have negative or NaN values.') 
-                print("""if using PISA type springs, this can be due to parameters behind out of the parameter space.
+                print("mt springs have negative or NaN values.")
+                print(
+                    """if using PISA type springs, this can be due to parameters behind out of the parameter space.
                 Please check that: 2 < L/D < 6.
-                """)            
+                """
+                )
             if check_springs(Hb):
-                print('Hb spring has negative or NaN values.') 
-                print("""if using PISA type springs, this can be due to parameters behind out of the parameter space.
+                print("Hb spring has negative or NaN values.")
+                print(
+                    """if using PISA type springs, this can be due to parameters behind out of the parameter space.
                 Please check that: 2 < L/D < 6.
-                """)
+                """
+                )
             if check_springs(Mb):
-                print('Mb spring has negative or NaN values.') 
-                print("""if using PISA type springs, this can be due to parameters behind out of the parameter space.
+                print("Mb spring has negative or NaN values.")
+                print(
+                    """if using PISA type springs, this can be due to parameters behind out of the parameter space.
                 Please check that: 2 < L/D < 6.
-                """)
+                """
+                )
             return py, mt, Hb, Mb, tz
 
         # creates mesh coordinates
@@ -1237,7 +1258,7 @@ def py_springs(self) -> pd.DataFrame:
         Returns
         -------
         pd.DataFrame
-            Table with p-y springs. 
+            Table with p-y springs.
         """
         return misc.get_springs(
             springs=self._py_springs,
@@ -1253,6 +1274,6 @@ def embedment(self) -> float:
         Returns
         -------
         float
-            Pile embedment 
+            Pile embedment
         """
-        return (self.soil.top_elevation - self.pile.bottom_elevation)
\ No newline at end of file
+        return self.soil.top_elevation - self.pile.bottom_elevation
diff --git a/src/openpile/core/kernel.py b/src/openpile/core/kernel.py
index 7781b89..f27db96 100644
--- a/src/openpile/core/kernel.py
+++ b/src/openpile/core/kernel.py
@@ -397,7 +397,6 @@ def elem_mt_stiffness_matrix(model, u, kind):
     d = model.element_properties["Diameter [m]"].to_numpy(dtype=float).reshape((-1, 1, 1))
     wt = model.element_properties["Wall thickness [m]"].to_numpy(dtype=float).reshape((-1, 1, 1))
 
-
     # calculate shear component in stiffness matrix (if Timorshenko)
     if model.element_type == "EulerBernoulli":
         N = 0 * L
@@ -423,14 +422,18 @@ def elem_mt_stiffness_matrix(model, u, kind):
 
         else:
             raise ValueError("Timoshenko beams cannot be used yet for non-circular pile types")
-        
-        phi = (12*omega+1)**2
+
+        phi = (12 * omega + 1) ** 2
 
         N = 0 * L
-        A = 6 * (120*omega**2 + 20*omega + 1) / ((5 * L) * phi) + 12*I / ( A * L**3 * phi)
-        B = 1 / (10 * phi) + 6*I / ( A*L**2 * phi)
-        C = (2 * L * (90*omega**2 + 15*omega +1) ) / ( 15 * phi ) + 4*I*(36*omega**2 + 6*omega + 1)/(A*L*phi)
-        D = - L * (360*omega**2 + 60*omega + 1) / (30 * phi) - 2*I*(72*omega**2 + 12*omega - 1)/(A*L*phi)
+        A = 6 * (120 * omega**2 + 20 * omega + 1) / ((5 * L) * phi) + 12 * I / (A * L**3 * phi)
+        B = 1 / (10 * phi) + 6 * I / (A * L**2 * phi)
+        C = (2 * L * (90 * omega**2 + 15 * omega + 1)) / (15 * phi) + 4 * I * (
+            36 * omega**2 + 6 * omega + 1
+        ) / (A * L * phi)
+        D = -L * (360 * omega**2 + 60 * omega + 1) / (30 * phi) - 2 * I * (
+            72 * omega**2 + 12 * omega - 1
+        ) / (A * L * phi)
 
     else:
         raise ValueError(
@@ -569,9 +572,9 @@ def struct_internal_force(model, u) -> np.ndarray:
         elif model.distributed_moment:
             k += elem_mt_stiffness_matrix(model, u, kind)
         elif model.base_shear:
-            k[-1,-2,-2] += calculate_base_spring_stiffness(u[-2], model._Hb_spring, kind)
+            k[-1, -2, -2] += calculate_base_spring_stiffness(u[-2], model._Hb_spring, kind)
         elif model.base_moment:
-            k[-1,-1,-1] += calculate_base_spring_stiffness(u[-1], model._Mb_spring, kind)
+            k[-1, -1, -1] += calculate_base_spring_stiffness(u[-1], model._Mb_spring, kind)
 
     # create array u of shape [n_elem x 6 x 1]
     u = global_dof_vector_to_consistent_stacked_array(u, ndof_per_node * node_per_element)
@@ -599,7 +602,7 @@ def calculate_base_spring_stiffness(
     Returns
     -------
     k: float
-        secant or tangent stiffness for all elements. 
+        secant or tangent stiffness for all elements.
     """
 
     # displacemet with same dimension as spring
@@ -634,6 +637,7 @@ def calculate_base_spring_stiffness(
 
     return k
 
+
 @njit(cache=True)
 def calculate_py_springs_stiffness(
     u: np.ndarray, springs: np.ndarray, kind: Literal["initial", "secant", "tangent"]
@@ -776,11 +780,11 @@ def calculate_mt_springs_stiffness(
                     else:
                         m1 = np.interp(dt, t, m)
                 elif kind == "tangent":
-                    dt = min(0.01*t[1], d[i, j, 0, 0])
-                    if (d[i, j, 0, 0]-dt) > t[-1]:
+                    dt = min(0.01 * t[1], d[i, j, 0, 0])
+                    if (d[i, j, 0, 0] - dt) > t[-1]:
                         m0 = m[-1]
                     else:
-                        m0 = np.interp(d[i, j, 0, 0]-dt, t, m)
+                        m0 = np.interp(d[i, j, 0, 0] - dt, t, m)
                     if (d[i, j, 0, 0]) > t[-1]:
                         m1 = m[-1]
                     else:
diff --git a/src/openpile/core/misc.py b/src/openpile/core/misc.py
index ed04029..d4765f0 100644
--- a/src/openpile/core/misc.py
+++ b/src/openpile/core/misc.py
@@ -107,33 +107,31 @@ def springs_to_df(springs: np.ndarray, elevations: np.ndarray, flag) -> pd.DataF
     return springs_to_df(springs, elevations, flag=kind)
 
 
-
-
 @njit(cache=True)
 def conic(
-        x_u : float,
-        n : float,
-        k : float,
-        y_u : float,
-        output_length : int,
+    x_u: float,
+    n: float,
+    k: float,
+    y_u: float,
+    output_length: int,
 ):
     # Create x vector with 10% extension
-    x = np.linspace(0,x_u,output_length-1).astype(np.float32)
-    x = np.append(x,1.1*x_u)
+    x = np.linspace(0, x_u, output_length - 1).astype(np.float32)
+    x = np.append(x, 1.1 * x_u)
 
-    a = 1-2*n
+    a = 1 - 2 * n
 
-    y = np.zeros((len(x)),dtype=np.float32)
+    y = np.zeros((len(x)), dtype=np.float32)
 
     for i in range(len(x)):
-        if abs(x[i]-x_u) < 1e-4:
+        if abs(x[i] - x_u) < 1e-4:
             y[i] = y_u
         elif x[i] < x_u:
-            b = 2*n*x[i]/x_u - (1-n)*(1+x[i]*k/y_u)
-            c = x[i] * (k / y_u) * (1-n) - n * (x[i]**2 / x_u**2)
+            b = 2 * n * x[i] / x_u - (1 - n) * (1 + x[i] * k / y_u)
+            c = x[i] * (k / y_u) * (1 - n) - n * (x[i] ** 2 / x_u**2)
 
-            y[i] = y_u * 2*c / (-b + ( b**2 - 4*a*c)**0.5 )
+            y[i] = y_u * 2 * c / (-b + (b**2 - 4 * a * c) ** 0.5)
         else:
             y[i] = y_u
 
-    return x, y
\ No newline at end of file
+    return x, y
diff --git a/src/openpile/soilmodels.py b/src/openpile/soilmodels.py
index c1aa2f1..3e12a1d 100644
--- a/src/openpile/soilmodels.py
+++ b/src/openpile/soilmodels.py
@@ -46,10 +46,11 @@ class LateralModel(ConstitutiveModel):
 class AxialModel(ConstitutiveModel):
     pass
 
+
 @dataclass(config=PydanticConfigFrozen)
 class API_clay(AxialModel):
     #: undrained shear strength [kPa], if a variation in values, two values can be given.
-    Su: Union[PositiveFloat, conlist(PositiveFloat, min_items=1, max_items=2)]    
+    Su: Union[PositiveFloat, conlist(PositiveFloat, min_items=1, max_items=2)]
     #: t-multiplier
     t_multiplier: confloat(ge=0.0) = 1.0
     #: z-multiplier
@@ -63,7 +64,6 @@ def __str__(self):
         return f"\tAPI clay\n\tSu = {var_to_str(self.Su)} kPa"
 
 
-
 @dataclass(config=PydanticConfigFrozen)
 class Cowden_clay(LateralModel):
     """A class to establish the PISA Cowden clay model.
@@ -127,7 +127,7 @@ def py_spring_fct(
         G0_t, G0_b = from_list2x_parse_top_bottom(self.G0)
         Gmax = G0_t + (G0_b - G0_t) * depth_from_top_of_layer / layer_height
 
-        p , y = py_curves.cowden_clay(
+        p, y = py_curves.cowden_clay(
             X=X,
             Su=Su,
             G0=Gmax,
@@ -194,7 +194,7 @@ def mt_spring_fct(
         G0_t, G0_b = from_list2x_parse_top_bottom(self.G0)
         Gmax = G0_t + (G0_b - G0_t) * depth_from_top_of_layer / layer_height
 
-        p_array , _ = py_curves.cowden_clay(
+        p_array, _ = py_curves.cowden_clay(
             X=X,
             Su=Su,
             G0=Gmax,
@@ -202,11 +202,11 @@ def mt_spring_fct(
             output_length=output_length,
         )
 
-        m = np.zeros((output_length,output_length),dtype=np.float32)
-        t = np.zeros((output_length,output_length),dtype=np.float32)
+        m = np.zeros((output_length, output_length), dtype=np.float32)
+        t = np.zeros((output_length, output_length), dtype=np.float32)
 
-        for count, _ in enumerate(p_array):  
-            m[count,:], t[count,:] = mt_curves.cowden_clay(
+        for count, _ in enumerate(p_array):
+            m[count, :], t[count, :] = mt_curves.cowden_clay(
                 X=X,
                 Su=Su,
                 G0=Gmax,
@@ -251,7 +251,6 @@ def Mb_spring_fct(
         return Mb, y
 
 
-
 @dataclass(config=PydanticConfigFrozen)
 class Dunkirk_sand(LateralModel):
     """A class to establish the PISA Dunkirk sand model.
@@ -385,7 +384,7 @@ def mt_spring_fct(
         G0_t, G0_b = from_list2x_parse_top_bottom(self.G0)
         Gmax = G0_t + (G0_b - G0_t) * depth_from_top_of_layer / layer_height
 
-        p_array , _ = py_curves.dunkirk_sand(
+        p_array, _ = py_curves.dunkirk_sand(
             sig=sig,
             X=X,
             Dr=Dr,
@@ -395,16 +394,16 @@ def mt_spring_fct(
             output_length=output_length,
         )
 
-        m = np.zeros((output_length,output_length),dtype=np.float32)
-        t = np.zeros((output_length,output_length),dtype=np.float32)
+        m = np.zeros((output_length, output_length), dtype=np.float32)
+        t = np.zeros((output_length, output_length), dtype=np.float32)
 
-        for count, p_iter in enumerate(p_array):  
-            m[count], t[count,:] = mt_curves.dunkirk_sand(
+        for count, p_iter in enumerate(p_array):
+            m[count], t[count, :] = mt_curves.dunkirk_sand(
                 sig=sig,
                 X=X,
                 Dr=Dr,
                 G0=Gmax,
-                p = p_iter,
+                p=p_iter,
                 D=D,
                 L=L,
                 output_length=output_length,
@@ -447,6 +446,7 @@ def Mb_spring_fct(
 
         return Mb, y
 
+
 @dataclass(config=PydanticConfigFrozen)
 class API_sand(LateralModel):
     """A class to establish the API sand model.
diff --git a/src/openpile/utils/Hb_curves.py b/src/openpile/utils/Hb_curves.py
index c65429e..9447df1 100644
--- a/src/openpile/utils/Hb_curves.py
+++ b/src/openpile/utils/Hb_curves.py
@@ -19,12 +19,12 @@ def dunkirk_sand(
     X: float,
     Dr: float,
     G0: float,
-    D: float, 
-    L: float, 
-    output_length: int = 20,  
+    D: float,
+    L: float,
+    output_length: int = 20,
 ):
     """
-    Creates the base shear spring from the PISA sand formulation 
+    Creates the base shear spring from the PISA sand formulation
     published by Burd et al (2020).
     Also called the General Dunkirk Sand Model (GDSM).
 
@@ -47,17 +47,17 @@ def dunkirk_sand(
 
     Returns
     -------
-    1darray 
+    1darray
         Hb vector [unit: kN]
-    1darray 
+    1darray
         y vector [unit: m]
     """
-    #correct relative density for decimal value
-    Dr = Dr/100
+    # correct relative density for decimal value
+    Dr = Dr / 100
 
     # Generalised Dunkirk Sand Model parameters
-    v_hu1 = 0.5150 + 2.883 * Dr 
-    v_hu2 = 0.1695 - 0.7018 * Dr 
+    v_hu1 = 0.5150 + 2.883 * Dr
+    v_hu2 = 0.1695 - 0.7018 * Dr
     k_h1 = 6.505 - 2.985 * Dr
     k_h2 = -0.007969 - 0.4299 * Dr
     n_h1 = 0.09978 + 0.7974 * Dr
@@ -66,17 +66,16 @@ def dunkirk_sand(
     p_u2 = 0.03988 - 0.1606 * Dr
 
     # Depth variation parameters
-    v_max = v_hu1 + v_hu2 * L/D
-    k = k_h1 + k_h2 * L/D
-    n = n_h1 + n_h2 * L/D
-    p_max = p_u1 + p_u2 * L/D
+    v_max = v_hu1 + v_hu2 * L / D
+    k = k_h1 + k_h2 * L / D
+    n = n_h1 + n_h2 * L / D
+    p_max = p_u1 + p_u2 * L / D
 
     # calculate normsalised conic function
     y, p = conic(v_max, n, k, p_max, output_length)
 
     # return non-normalised curve
-    return p*(sig*D**2), y*(sig*D/G0)
-
+    return p * (sig * D**2), y * (sig * D / G0)
 
 
 @njit(cache=True)
@@ -84,13 +83,13 @@ def cowden_clay(
     X: float,
     Su: float,
     G0: float,
-    D: float, 
+    D: float,
     L: float,
-    output_length: int = 20,  
+    output_length: int = 20,
 ):
     """
-    Creates the base shear spring from the PISA clay formulation 
-    published by Byrne et al (2020) and calibrated based pile 
+    Creates the base shear spring from the PISA clay formulation
+    published by Byrne et al (2020) and calibrated based pile
     load tests at Cowden (north east coast of England).
 
     Parameters
@@ -110,9 +109,9 @@ def cowden_clay(
 
     Returns
     -------
-    1darray 
+    1darray
         Hb vector [unit: kN]
-    1darray 
+    1darray
         y vector [unit: m]
     """
 
@@ -127,13 +126,13 @@ def cowden_clay(
     p_u2 = 0.04812
 
     # Depth variation parameters
-    v_max = v_hu1 + v_hu2 * L/D
-    k = k_h1 + k_h2 * L/D
-    n = n_h1 + n_h2 * L/D
-    p_max = p_u1 + p_u2 * L/D
+    v_max = v_hu1 + v_hu2 * L / D
+    k = k_h1 + k_h2 * L / D
+    n = n_h1 + n_h2 * L / D
+    p_max = p_u1 + p_u2 * L / D
 
     # calculate normsalised conic function
     y, p = conic(v_max, n, k, p_max, output_length)
 
     # return non-normalised curve
-    return p*(Su*D**2), y*(Su*D/G0)
\ No newline at end of file
+    return p * (Su * D**2), y * (Su * D / G0)
diff --git a/src/openpile/utils/Mb_curves.py b/src/openpile/utils/Mb_curves.py
index ca4efe9..5c21516 100644
--- a/src/openpile/utils/Mb_curves.py
+++ b/src/openpile/utils/Mb_curves.py
@@ -18,13 +18,13 @@ def cowden_clay(
     X: float,
     Su: float,
     G0: float,
-    D: float, 
+    D: float,
     L: float,
-    output_length: int = 20,  
+    output_length: int = 20,
 ):
     """
-    Create the base moment springs from the PISA clay formulation 
-    published by Byrne et al (2020) and calibrated based pile load tests 
+    Create the base moment springs from the PISA clay formulation
+    published by Byrne et al (2020) and calibrated based pile load tests
     at Cowden (north east coast of England).
 
     Parameters
@@ -44,9 +44,9 @@ def cowden_clay(
 
     Returns
     -------
-    1darray 
+    1darray
         Mb vector [unit: kN]
-    1darray 
+    1darray
         t vector of length [unit: rad]
     """
 
@@ -60,16 +60,16 @@ def cowden_clay(
     psi_u = 173.1
 
     # Depth variation parameters
-    k = k_m1 + k_m2 * L/D
-    n = n_m1 + n_m2 * L/D
-    m_max = m_m1 + m_m2 * L/D
+    k = k_m1 + k_m2 * L / D
+    n = n_m1 + n_m2 * L / D
+    m_max = m_m1 + m_m2 * L / D
     psi_max = psi_u
 
     # calculate normsalised conic function
     t, m = conic(psi_max, n, k, m_max, output_length)
 
     # return non-normalised curve
-    return m*(Su*D**3), t*(Su/G0)
+    return m * (Su * D**3), t * (Su / G0)
 
 
 @njit(cache=True)
@@ -78,14 +78,14 @@ def dunkirk_sand(
     X: float,
     Dr: float,
     G0: float,
-    D: float, 
-    L: float, 
-    output_length: int = 20,  
+    D: float,
+    L: float,
+    output_length: int = 20,
 ):
     """
-    Create the base moment spring from the PISA sand 
+    Create the base moment spring from the PISA sand
     formulation published by Burd et al (2020).
-    Also called the General Dunkirk Sand Model (GDSM) 
+    Also called the General Dunkirk Sand Model (GDSM)
 
     Parameters
     ----------
@@ -106,13 +106,13 @@ def dunkirk_sand(
 
     Returns
     -------
-    1darray 
+    1darray
         Mb vector [unit: kN]
-    1darray 
+    1darray
         t vector of length [unit: rad]
     """
-    #correct relative density for decimal value
-    Dr = Dr/100
+    # correct relative density for decimal value
+    Dr = Dr / 100
 
     # Generalised Dunkirk Sand Model parameters
     k_m1 = 0.3515
@@ -123,13 +123,13 @@ def dunkirk_sand(
     psi_u = 44.89
 
     # Depth variation parameters
-    k = k_m1 + k_m2 * L/D
+    k = k_m1 + k_m2 * L / D
     n = n_m
-    m_max = m_u1 + m_u2 * L/D
+    m_max = m_u1 + m_u2 * L / D
     psi_max = psi_u
 
     # calculate normsalised conic function
     t, m = conic(psi_max, n, k, m_max, output_length)
 
     # return non-normalised curve
-    return m*(sig*D**3), t*(sig/G0)
+    return m * (sig * D**3), t * (sig / G0)
diff --git a/src/openpile/utils/mt_curves.py b/src/openpile/utils/mt_curves.py
index 94a5d8e..e61c805 100644
--- a/src/openpile/utils/mt_curves.py
+++ b/src/openpile/utils/mt_curves.py
@@ -18,12 +18,12 @@ def cowden_clay(
     X: float,
     Su: float,
     G0: float,
-    D: float, 
-    output_length: int = 20,  
+    D: float,
+    output_length: int = 20,
 ):
     """
-    Create the rotational springs from the PISA clay formulation 
-    published by Byrne et al (2020) and calibrated based pile load tests 
+    Create the rotational springs from the PISA clay formulation
+    published by Byrne et al (2020) and calibrated based pile load tests
     at Cowden (north east coast of England).
 
     Parameters
@@ -41,9 +41,9 @@ def cowden_clay(
 
     Returns
     -------
-    1darray 
+    1darray
         m vector [unit: kN]
-    1darray 
+    1darray
         t vector of length [unit: rad]
     """
 
@@ -56,16 +56,16 @@ def cowden_clay(
     m_m2 = -0.04775
 
     # Depth variation parameters
-    k = k_m1 + k_m2 * X/D
-    n = n_m1 + n_m2 * X/D
-    m_max = m_m1 + m_m2 * X/D
+    k = k_m1 + k_m2 * X / D
+    n = n_m1 + n_m2 * X / D
+    m_max = m_m1 + m_m2 * X / D
     psi_max = m_max / k
 
     # calculate normsalised conic function
     t, m = conic(psi_max, n, k, m_max, output_length)
 
     # return non-normalised curve
-    return m*(Su*D**2), t*(Su/G0)
+    return m * (Su * D**2), t * (Su / G0)
 
 
 @njit(cache=True)
@@ -75,14 +75,14 @@ def dunkirk_sand(
     Dr: float,
     G0: float,
     p: float,
-    D: float, 
-    L: float, 
-    output_length: int = 20,  
+    D: float,
+    L: float,
+    output_length: int = 20,
 ):
     """
-    Create the rotational springs from the PISA sand 
+    Create the rotational springs from the PISA sand
     formulation published by Burd et al (2020).
-    Also called the General Dunkirk Sand Model (GDSM) 
+    Also called the General Dunkirk Sand Model (GDSM)
 
     Parameters
     ----------
@@ -95,7 +95,7 @@ def dunkirk_sand(
     G0 : float
         Small-strain shear modulus [unit: kPa]
     p : float
-        radial stress computed via p-y curves [unit: kN/m] 
+        radial stress computed via p-y curves [unit: kN/m]
     D : float
         Pile diameter [unit: m]
     L : float
@@ -105,18 +105,17 @@ def dunkirk_sand(
 
     Returns
     -------
-    1darray 
+    1darray
         m vector [unit: kN]
-    1darray 
+    1darray
         t vector of length [unit: rad]
     """
-    #correct relative density for decimal value
-    Dr = Dr/100
+    # correct relative density for decimal value
+    Dr = Dr / 100
 
     # correct p value
     p = abs(p)
 
-
     # Generalised Dunkirk Sand Model parameters
     k_m1 = 17.00
     k_m2 = 0.00
@@ -125,13 +124,13 @@ def dunkirk_sand(
     m_u2 = -0.1989 + 0.2019 * Dr
 
     # Depth variation parameters
-    k = k_m1 + k_m2 * X/D
+    k = k_m1 + k_m2 * X / D
     n = n_m
-    m_max = m_u1 + m_u2 * X/L
+    m_max = m_u1 + m_u2 * X / L
     psi_max = m_max / k
 
     # calculate normsalised conic function
     t, m = conic(psi_max, n, k, m_max, output_length)
 
     # return non-normalised curve
-    return m*(p*D), t*(sig/G0)
\ No newline at end of file
+    return m * (p * D), t * (sig / G0)
diff --git a/src/openpile/utils/py_curves.py b/src/openpile/utils/py_curves.py
index 7e2b46c..e9afc3e 100644
--- a/src/openpile/utils/py_curves.py
+++ b/src/openpile/utils/py_curves.py
@@ -18,12 +18,12 @@ def cowden_clay(
     X: float,
     Su: float,
     G0: float,
-    D: float, 
-    output_length: int = 20,  
+    D: float,
+    output_length: int = 20,
 ):
     """
-    Creates the lateral springs from the PISA clay formulation 
-    published by Byrne et al (2020) and calibrated based pile 
+    Creates the lateral springs from the PISA clay formulation
+    published by Byrne et al (2020) and calibrated based pile
     load tests at Cowden (north east coast of England).
 
     Parameters
@@ -41,14 +41,14 @@ def cowden_clay(
 
     Returns
     -------
-    1darray 
+    1darray
         p vector [unit: kN/m]
-    1darray 
+    1darray
         y vector [unit: m]
     """
 
     # # Cowden clay parameters
-    v_pu = 241.4 
+    v_pu = 241.4
     k_p1 = 10.6
     k_p2 = -1.650
     n_p1 = 0.9390
@@ -58,15 +58,15 @@ def cowden_clay(
 
     # Depth variation parameters
     v_max = v_pu
-    k = k_p1 + k_p2 * X/D
-    n = n_p1 + n_p2 * X/D
-    p_max = p_u1 + p_u2 * m.exp(-0.3085*X/D)
+    k = k_p1 + k_p2 * X / D
+    n = n_p1 + n_p2 * X / D
+    p_max = p_u1 + p_u2 * m.exp(-0.3085 * X / D)
 
     # calculate normsalised conic function
     y, p = conic(v_max, n, k, p_max, output_length)
 
     # return non-normalised curve
-    return p*(Su*D), y*(Su*D/G0)
+    return p * (Su * D), y * (Su * D / G0)
 
 
 @njit(cache=True)
@@ -75,12 +75,12 @@ def dunkirk_sand(
     X: float,
     Dr: float,
     G0: float,
-    D: float, 
-    L: float, 
-    output_length: int = 20,  
+    D: float,
+    L: float,
+    output_length: int = 20,
 ):
     """
-    Creates the lateral spring from the PISA sand formulation 
+    Creates the lateral spring from the PISA sand formulation
     published by Burd et al (2020).
     Also called the General Dunkirk Sand Model (GDSM).
 
@@ -103,16 +103,16 @@ def dunkirk_sand(
 
     Returns
     -------
-    1darray 
+    1darray
         p vector [unit: kN/m]
-    1darray 
+    1darray
         y vector [unit: m]
     """
-    #correct relative density for decimal value
-    Dr = Dr/100
+    # correct relative density for decimal value
+    Dr = Dr / 100
 
     # Generalised Dunkirk Sand Model parameters
-    v_pu = 146.1 - 92.11 * Dr 
+    v_pu = 146.1 - 92.11 * Dr
     k_p1 = 8.731 - 0.6982 * Dr
     k_p2 = -0.9178
     n_p = 0.917 + 0.06193 * Dr
@@ -121,15 +121,16 @@ def dunkirk_sand(
 
     # Depth variation parameters
     v_max = v_pu
-    k = k_p1 + k_p2 * X/D
+    k = k_p1 + k_p2 * X / D
     n = n_p
-    p_max = p_u1 + p_u2 * X/L
+    p_max = p_u1 + p_u2 * X / L
 
     # calculate normsalised conic function
     y, p = conic(v_max, n, k, p_max, output_length)
 
     # return non-normalised curve
-    return p*(sig*D), y*(sig*D/G0)
+    return p * (sig * D), y * (sig * D / G0)
+
 
 # API sand function
 @njit(parallel=True, cache=True)
@@ -164,12 +165,12 @@ def api_sand(
         maximum value of y, default goes to 99.9% of ultimate resistance
     output_length: int, by default 20
         Number of discrete point along the springs
-    
+
     Returns
     -------
-    1darray 
+    1darray
         p vector [unit: kN/m]
-    1darray 
+    1darray
         y vector [unit: m]
     """
     # A value - only thing that changes between cyclic or static
@@ -268,12 +269,12 @@ def api_clay(
         maximum value of y, if null the maximum is calculated such that the whole curve is computed
     output_length: int, by default 20
         Number of discrete point along the springs
-    
+
     Returns
     -------
-    1darray 
+    1darray
         p vector [unit: kN/m]
-    1darray 
+    1darray
         y vector [unit: m]
     ---------
     """
diff --git a/src/openpile/utils/qz_curves.py b/src/openpile/utils/qz_curves.py
index 103522e..f6a94a1 100644
--- a/src/openpile/utils/qz_curves.py
+++ b/src/openpile/utils/qz_curves.py
@@ -41,7 +41,7 @@ def api_sand(
     numpy 1darray
         z vector [unit: m]
     """
-    #cannot have less than 7
+    # cannot have less than 7
     if output_length < 7:
         output_length = 7
 
@@ -50,15 +50,15 @@ def api_sand(
     Nq_table = np.array([8, 8, 12, 20, 40, 50, 50], dtype=np.float32)
     Qmax_table = np.array([1900, 1900, 2900, 4800, 9600, 12000, 12000], dtype=np.float32)
 
-    Nq = np.interp(delta, delta_table, Nq_table) 
-    Qmax = np.interp(delta, delta_table, Qmax_table) 
+    Nq = np.interp(delta, delta_table, Nq_table)
+    Qmax = np.interp(delta, delta_table, Qmax_table)
 
     # Unit end-bearing [kPa]
-    f = min(Qmax, sig*Nq)
+    f = min(Qmax, sig * Nq)
 
     # piecewise function
     zlist = [0.0, 0.002, 0.013, 0.042, 0.073, 0.100, 0.200]
-    Qlist = [0.0 , 0.25, 0.50, 0.75, 0.90, 1.00, 1.00]
+    Qlist = [0.0, 0.25, 0.50, 0.75, 0.90, 1.00, 1.00]
 
     # determine z vector
     z = np.array(zlist, dtype=np.float32) * D
@@ -71,11 +71,11 @@ def api_sand(
     add_Q_values = np.zeros((add_values), dtype=np.float32)
 
     for i in range(add_values):
-        add_z_values[i] = (0.1 + random()*0.1) * D
-        add_Q_values[i] =  f
+        add_z_values[i] = (0.1 + random() * 0.1) * D
+        add_Q_values[i] = f
 
-    z = np.append(z,add_z_values)
-    Q = np.append(Q,add_Q_values)
+    z = np.append(z, add_z_values)
+    Q = np.append(Q, add_Q_values)
 
     z = np.sort(z)
     z_id_sorted = np.argsort(z)
@@ -111,16 +111,16 @@ def api_clay(
     numpy 1darray
         z vector [unit: m]
     """
-    #cannot have less than 7
+    # cannot have less than 7
     if output_length < 7:
         output_length = 7
 
     # Unit end-bearing [kPa]
-    f = 9*Su
+    f = 9 * Su
 
     # piecewise function
     zlist = [0.0, 0.002, 0.013, 0.042, 0.073, 0.100, 0.200]
-    Qlist = [0.0 , 0.25, 0.50, 0.75, 0.90, 1.00, 1.00]
+    Qlist = [0.0, 0.25, 0.50, 0.75, 0.90, 1.00, 1.00]
 
     # determine z vector
     z = np.array(zlist, dtype=np.float32) * D
@@ -133,15 +133,15 @@ def api_clay(
     add_Q_values = np.zeros((add_values), dtype=np.float32)
 
     for i in range(add_values):
-        add_z_values[i] = (0.1 + random()*0.1) * D
-        add_Q_values[i] =  f
+        add_z_values[i] = (0.1 + random() * 0.1) * D
+        add_Q_values[i] = f
 
-    z = np.append(z,add_z_values)
-    Q = np.append(Q,add_Q_values)
+    z = np.append(z, add_z_values)
+    Q = np.append(Q, add_Q_values)
 
     z = np.sort(z)
     z_id_sorted = np.argsort(z)
 
     Q = Q[z_id_sorted]
 
-    return Q, z
\ No newline at end of file
+    return Q, z
diff --git a/src/openpile/utils/tz_curves.py b/src/openpile/utils/tz_curves.py
index db39669..24eafcd 100644
--- a/src/openpile/utils/tz_curves.py
+++ b/src/openpile/utils/tz_curves.py
@@ -33,8 +33,8 @@ def api_clay(
     D: float
         Pile diameter [unit: m]
     residual: float
-        residual strength after peak strength, according to API-RP-2A, 
-        this value is between 0.7 and 0.9, default to 0.9 
+        residual strength after peak strength, according to API-RP-2A,
+        this value is between 0.7 and 0.9, default to 0.9
     output_length : int, optional
         Number of discrete point along the springs, cannot be lower than 7, by default 7
 
@@ -45,7 +45,7 @@ def api_clay(
     numpy 1darray
         z vector [unit: m]
     """
-    #cannot have less than 8
+    # cannot have less than 8
     if output_length < 8:
         output_length = 8
 
@@ -56,16 +56,16 @@ def api_clay(
         psi = Su / sig
 
     if psi > 1.0:
-        alpha = min( 0.5 * psi**(-0.25), 1.0 )
+        alpha = min(0.5 * psi ** (-0.25), 1.0)
     else:
-        alpha = min( 0.5 * psi**(-0.5), 1.0 )
+        alpha = min(0.5 * psi ** (-0.5), 1.0)
 
     # Unit skin friction [kPa]
     f = alpha * Su
 
     # piecewise function
     zlist = [0.0, 0.0016, 0.0031, 0.0057, 0.0080, 0.0100, 0.0200, 0.0300]
-    tlist = [0.0 , 0.3, 0.5, 0.75, 0.90, 1.00, residual, residual]
+    tlist = [0.0, 0.3, 0.5, 0.75, 0.90, 1.00, residual, residual]
 
     # determine z vector
     z = np.array(zlist, dtype=np.float32) * D
@@ -78,11 +78,11 @@ def api_clay(
     add_t_values = np.zeros((add_values), dtype=np.float32)
 
     for i in range(add_values):
-        add_z_values[i] = (0.02 + random()*0.01) *D 
+        add_z_values[i] = (0.02 + random() * 0.01) * D
         add_t_values[i] = residual * f
 
-    z = np.append(z,add_z_values)
-    t = np.append(t,add_t_values)
+    z = np.append(z, add_z_values)
+    t = np.append(t, add_t_values)
 
     z = np.sort(z)
     z_id_sorted = np.argsort(z)
@@ -121,7 +121,7 @@ def api_sand(
     numpy 1darray
         z vector [unit: m]
     """
-    #cannot have less than 4
+    # cannot have less than 4
     if output_length < 4:
         output_length = 4
 
@@ -130,14 +130,14 @@ def api_sand(
     fs_max_table = np.array([47.8, 47.8, 67, 81.3, 95.7, 114.8, 114.8], dtype=np.float32)
 
     # limit unit skin friction according to API ref page 59
-    fs_max = np.interp(delta, delta_table, fs_max_table) 
+    fs_max = np.interp(delta, delta_table, fs_max_table)
 
     # Unit skin friction [kPa]
-    f = min(fs_max, K*sig*m.tan(delta*m.pi/180.0))
+    f = min(fs_max, K * sig * m.tan(delta * m.pi / 180.0))
 
     # piecewise function
     zlist = [0.0, 0.0254, 0.03, 0.04]
-    tlist = [0.0 , 1.0, 1.0, 1.0]
+    tlist = [0.0, 1.0, 1.0, 1.0]
 
     # determine z vector
     z = np.array(zlist, dtype=np.float32)
@@ -150,11 +150,11 @@ def api_sand(
     add_t_values = np.zeros((add_values), dtype=np.float32)
 
     for i in range(add_values):
-        add_z_values[i] = (0.03 + random()*0.01)
-        add_t_values[i] =  f
+        add_z_values[i] = 0.03 + random() * 0.01
+        add_t_values[i] = f
 
-    z = np.append(z,add_z_values)
-    t = np.append(t,add_t_values)
+    z = np.append(z, add_z_values)
+    t = np.append(t, add_t_values)
 
     z = np.sort(z)
     z_id_sorted = np.argsort(z)
@@ -162,5 +162,3 @@ def api_sand(
     t = t[z_id_sorted]
 
     return t, z
-
-