diff --git a/halo_mass_function.py b/halo_mass_function.py
index 43b5edf..104c913 100644
--- a/halo_mass_function.py
+++ b/halo_mass_function.py
@@ -230,7 +230,7 @@ def sigma_Pk(self, R):
            parameter from 0 to infinity.
            Note that R is in h^-1 Mpc (comoving)
         """
-        result = integ.quad(self.sigma_squared_integrand,0,np.Inf, args=(R,),epsrel=1e-2)[0]
+        result = integ.quad(self.sigma_squared_integrand,0,np.inf, args=(R,),epsrel=1e-2)[0]
         return result / (2.0 * math.pi**2)
 
     def sigma_squared_of_R(self, R):
@@ -258,7 +258,7 @@ def _logsigma_of_R(self, R):
         """For the halo mass function we also need M * d log sigma/dM.
         This routine computes a table of that quantity.
         We use d log sigma /dM = 1/(2 sigma^2) dR /dM int(k^2 P(k) d/dR(W^2 (kR) dk"""
-        result = integ.quad(self.dlogsigma_integrand,0,np.Inf, args=(R,),epsrel=1e-2)[0]
+        result = integ.quad(self.dlogsigma_integrand,0,np.inf, args=(R,),epsrel=1e-2)[0]
         return result / (2.0 * math.pi**2) * R / (self.sigma_squared_of_R(R) * 3.) - self.use_pbh*self.sigma_square_poisson(R) / (2*self.sigma_squared_of_R(R))
 
     def dlogsigma_integrand(self,k,R):
diff --git a/matrices.py b/matrices.py
new file mode 100644
index 0000000..974a4f3
--- /dev/null
+++ b/matrices.py
@@ -0,0 +1,83 @@
+import time
+import numpy as np
+
+
+def multiply_matrices(A: list[list[float]], B: list[list[float]]) -> list[list[float]]:
+    """
+    Multiply two matrices A and B using nested for loops.
+    
+    :param A: Matrix of size l x m
+    :type A: list[list[float]]
+    :param B: Matrix of size m x n
+    :type B: list[list[float]]
+    :return: Matrix C of size l x n
+    :rtype: list[list[float]]
+    """
+
+    assert len(A[0]) == len(B), "The columns of A do not match the rows of B."
+
+    C = []
+    for row in range(len(A)):
+        C.append([0] * len(B[0]))
+
+    for row in range(len(A)):
+        for column in range(len(B[0])):
+            for k in range(len(B)):
+                C[row][column] += A[row][k] * B[k][column]
+
+    return C
+
+def multiply_matrices_list(A: list[list[float]], B: list[list[float]]) -> list[list[float]]:
+    """
+    Multiply two matrices A and B using nested list comprehension.
+    
+    :param A: Matrix of size l x m
+    :type A: list[list[float]]
+    :param B: Matrix of size m x n
+    :type B: list[list[float]]
+    :return: Matrix C of size l x n
+    :rtype: list[list[float]]
+    """
+
+    assert len(A[0]) == len(B), "The columns of A do not match the rows of B."
+
+    C = [[sum(a*b for a, b in zip(A_row, B_column)) for B_column in zip(*B)] for A_row in A]
+
+    return C
+
+def multiply_matrices_np(A: list[list[float]], B: list[list[float]]) -> list[list[float]]:
+    """
+    Multiply matrices using numpy.
+    
+    :param A: Matrix of size l x m
+    :type A: list[list[float]]
+    :param B: Matrix of size m x n
+    :type B: list[list[float]]
+    :return: Matrix C of size l x n
+    :rtype: list[list[float]]
+    """
+
+    assert np.shape(A)[1] == np.shape(B)[0], "The columns of A do not match the rows of B."
+
+    C = np.dot(A, B)
+
+    return C
+
+# compute time
+A = np.random.rand(50,50) + np.eye(50) # A 50x50 matrix containing random integers
+B = np.linalg.inv(A) # Matrix inverse
+
+start_for = time.time()
+print(np.allclose(multiply_matrices(A, B), np.identity(len(A))))
+end_for = time.time()
+print("For loop time:", end_for - start_for, "seconds")
+
+start_list = time.time()
+print(np.allclose(multiply_matrices_list(A, B), np.identity(len(A))))
+end_list = time.time()
+print("List time:", end_list - start_list, "seconds")
+
+start_np = time.time()
+print(np.allclose(multiply_matrices_np(A, B), np.identity(len(A))))
+end_np = time.time()
+print("np time:", end_np - start_np, "seconds")
\ No newline at end of file
diff --git a/pbhmergers.py b/pbhmergers.py
index e9617d7..c651adc 100644
--- a/pbhmergers.py
+++ b/pbhmergers.py
@@ -63,7 +63,7 @@ def rhocrit(self):
         hubz2 = (self.overden.omega_matter0/aa**3 + self.overden.omega_lambda0) * hubble**2
         #Critical density at redshift in units of kg m^-3
         rhocrit = 3 * hubz2 / (8*math.pi* self.ureg.newtonian_constant_of_gravitation)
-        print "rhocrit = ", rhocrit
+        print("rhocrit = ", rhocrit) # Bug #1
         return rhocrit.to_base_units()
 
     def R200(self, mass):
@@ -224,7 +224,7 @@ def mergerhalflife(self,mass,threefac=True, bhmass=None):
             threefac = self.threebodyratio(mass)
             threefac = np.max([threefac, np.ones_like(threefac)],axis=0)
             rate *= threefac
-        return 0.5*(mass/bhmass)/rat
+        return 0.5*(mass/bhmass) # Bug #2
 
     def bias(self,mass):
         """The formula for halo bias in EPS theory (Mo & White 1996), eq. 13"""
@@ -313,7 +313,7 @@ def plot_concentration_vs_mass(redshift):
     mass = np.logspace(2,16)
     hh = NFWHalo(redshift)
     plt.loglog(mass, hh.concentration(mass), ls='-', label="Ludlow concentration")
-    hh.conc_model = concentration.PradaConcentration(hh.overden.omega_matter0)
+    hh.conc_model = concentration.LudlowConcentration(hh.overden.Dofz) # Bug #3 (?)
     plt.loglog(mass, hh.concentration(mass), ls='--', label="Prada concentration")
     plt.xlim(500,1e16)
     plt.xticks(np.logspace(3,15,5))
diff --git a/power_specs.py b/power_specs.py
index 1e00346..bb7b3ba 100644
--- a/power_specs.py
+++ b/power_specs.py
@@ -638,7 +638,6 @@ def plot_z(self,Knot,redshift,title="",ylabel=""):
         """Get a power spectrum in the flat format we use"""
         """for inputting some cosmomc tables"""
         def GetFlat(self,dir):
-                Pk_sdss=np.empty([11, 12])
                 #For z=2.07 we need to average snap_011 and snap_010
                 z=2.07
                 (k,PF_a)=self.loadpk(dir+self.suf+"snapshot_011"+self.ext, self.box)
@@ -660,12 +659,6 @@ def GetFlat(self,dir):
         flux=flux_pow()
         matter=matter_pow()
         fpdf=flux_pdf()
-        A_knot=knot(("A0.54/","A0.74/","A0.84/","A1.04/", "A1.14/","A1.34/"), (0.54,0.74,0.84,1.04,1.14,1.34),0.94,"best-fit/", 60)
-        AA_knot=knot(("AA0.54/","AA0.74/","AA1.14/","AA1.34/"), (0.54,0.74,1.14,1.34),0.94,"boxcorr400/", 120)
-        B_knot=knot(("B0.33/","B0.53/","B0.73/", "B0.83/", "B1.03/","B1.13/", "B1.33/"), (0.33,0.53,0.73,0.83,1.03, 1.13,1.33),0.93,"best-fit/", 60)
-        C_knot=knot(("C0.11/", "C0.31/","C0.51/","C0.71/","C1.11/","C1.31/","C1.51/"),(0.11, 0.31,0.51,0.71, 1.11,1.31,1.51),0.91,"bf2/", 60)
-        D_knot=knot(("D0.50/","D0.70/","D1.10/","D1.20/", "D1.30/", "D1.50/", "D1.70/"),(0.50, 0.70,1.10,1.20,1.30, 1.50, 1.70),0.90,"bfD/", 48)
-        interp=flux_interp(flux, (AA_knot, B_knot, C_knot, D_knot))
 #Thermal parameters for models
 
 #Models where gamma is varied
@@ -692,6 +685,3 @@ def GetFlat(self,dir):
         bf2ag=np.array([1.0184905 ,1.0332849 ,1.0404092 ,1.0538235 ,1.0598905 ,1.0695964 ,1.0771414 ,1.0827464 ,1.0904784 ,1.0964669 ,1.1009428,1.1068357])
         bf2at=np.array([26914.856, 27497.555, 28407.030, 28357.785, 28919.376, 28610.038, 28478.749, 28484.856, 27841.160, 27335.328, 26933.453,26099.643])/1e3
 
-        G_knot=knot(("bf2z/","bf2b/","bf2c/","bf2d/", "bf2T15/","bf2T35/","bf2T45/","bf2a/"),(bf2zg,bf2bg,bf2cg,bf2dg,T15g,T35g,T45g,bf2ag),bf2g,"bf2/",60, (bf2zt,bf2bt,bf2ct,bf2dt,T15t,T35t,T45t,bf2at),bf2t)
-        g_int=flux_interp(flux,(G_knot))
-
diff --git a/schwarzchild_radius.py b/schwarzchild_radius.py
new file mode 100644
index 0000000..da22ab7
--- /dev/null
+++ b/schwarzchild_radius.py
@@ -0,0 +1,18 @@
+import pint
+
+
+ureg = pint.UnitRegistry()
+
+# Define constants
+# https://physics.nist.gov/cgi-bin/cuu/Value?bg
+G = 6.67430E-11 * ureg.meter**3 / ureg.kilogram / ureg.second**2 # Newton's constant
+# https://web.archive.org/web/20140811195806/http://www.bipm.org/en/si/new_si/explicit_constant.html
+c = 299792458 * ureg.meter / ureg.second # speed of light
+# https://www.mccc.edu/~dornemam/Planet_Walk/Sun/the_sun.htm#:~:text=The%20sun%20has%20a%20mass%20of%201.9891x1030,it%20could%20hold%20over%201%20million%20Earths.
+M_sun = 1.9891E30 * ureg.kilogram # solar mass
+
+
+if __name__ == "__main__":
+    # Kutner 2003, p. 148
+    # https://archive.org/details/astronomyphysica00kutn/mode/2up
+    print(2*G*M_sun/(c**2)) # Schwarzchild radius formula
\ No newline at end of file
diff --git a/test_matrices.py b/test_matrices.py
new file mode 100644
index 0000000..5022dad
--- /dev/null
+++ b/test_matrices.py
@@ -0,0 +1,21 @@
+import numpy as np
+from matrices import *
+
+
+def test_multiply_matrices():
+    A = np.random.rand(50,50) + np.eye(50) # A 50x50 matrix containing random integers
+    B = np.linalg.inv(A) # Matrix inverse
+
+    assert np.allclose(multiply_matrices(A,B), np.identity(len(A))), "Matrix multiplication fails."
+
+def test_multiply_matrices_list():
+    A = np.random.rand(50,50) + np.eye(50) # A 50x50 matrix containing random integers
+    B = np.linalg.inv(A) # Matrix inverse
+
+    assert np.allclose(multiply_matrices_list(A,B), np.identity(len(A))), "Matrix multiplication fails."
+
+def test_multiply_matrices_np():
+    A = np.random.rand(50,50) + np.eye(50) # A 50x50 matrix containing random integers
+    B = np.linalg.inv(A) # Matrix inverse
+
+    assert np.allclose(multiply_matrices_np(A,B), np.identity(len(A))), "Matrix multiplication fails."
\ No newline at end of file