Hui-Hong HW1

Problem1.3/cal_Schwarzchild_radius.py

@@ -0,0 +1,15 @@
+"Problem1.3(a) Write a python script, using pint, which finds the Schwarzchild radius of the Sun, in m."
+import numpy as np
+import pint
+from astropy import constants as const
+
+units=pint.UnitRegistry()
+def Cal_R_Sch(m_in_solarms):
+    const_G=const.G.value*units.meter**3/(units.kilogram*units.second**2)
+    const_solarmass=const.M_sun.value*units.kilogram
+    const_sol=const.c.value*units.meter/units.second
+
+    return 2*const_G*const_solarmass*m_in_solarms/const_sol**2
+
+radius=Cal_R_Sch(1).magnitude
+print("The Schwarzchild radius of the Sun is "+"%.2f"%radius+"m.")

Problem1.3/matrix_multiply.py

@@ -0,0 +1,66 @@
+"Write a simple routine to multiply two matrices together in python"
+import numpy as np
+import time
+
+"routine using nested for loops"
+def matrix_mul_1(Matrix_A,Matrix_B):
+    rowA=len(Matrix_A);colA=len(Matrix_A[0])
+    rowB=len(Matrix_B);colB=len(Matrix_B[0])
+
+    Matrix_C=[]
+    for i in range(rowA):
+        row_C=[]
+        for j in range(colB):
+            sum=0
+            for m in range(colA):
+                sum+=Matrix_A[i][m]*Matrix_B[m][j]
+            row_C.append(sum)
+        Matrix_C.append(row_C)
+
+    return np.array(Matrix_C)
+
+"routine using a list comprehension"
+def matrix_mul_2(Matrix_A,Matrix_B):
+   return np.array([[np.sum(np.array([Matrix_A[j][m]*Matrix_B[m][i] for m in range(len(Matrix_A[0]))])) for i in range(len(Matrix_B[0]))] for j in range(len(Matrix_A))])
+
+"routine using the built-in numpy matrix multiplication"
+def matrix_mul_3(Matrix_A,Matrix_B):
+    return np.dot(Matrix_A,Matrix_B)
+
+
+np.set_printoptions(threshold=10)
+"test"
+"generate dimension=dim diagonal matrix so that it is inversible."
+dim=100
+mA=np.diag([i+1 for i in range(dim)])
+mB=np.linalg.inv(mA)
+print("Matrix A is")
+print(mA)
+print("Matrix B is")
+print(mB)
+
+"test1"
+print("test1 for routine using nested for loops")
+print("Multiplication of A and B is")
+stime=time.time()
+print(matrix_mul_1(mA,mB))
+etime=time.time()
+print("time="+"%.5f"%(etime-stime)+" s")
+print('\n')
+
+"test2"
+print("test2 using a list comprehension")
+print("Multiplication of A and B is")
+stime=time.time()
+print(matrix_mul_2(mA,mB))
+etime=time.time()
+print("time="+"%.5f"%(etime-stime)+" s")
+print('\n')
+
+"test3"
+print("test3 using the built-in numpy matrix multiplication")
+print("Multiplication of A and B is")
+stime=time.time()
+print(matrix_mul_3(mA,mB))
+etime=time.time()
+print("time="+"%.5f"%(etime-stime)+" s")

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)
         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)/rate
 
     def bias(self,mass):
         """The formula for halo bias in EPS theory (Mo & White 1996), eq. 13"""
@@ -319,7 +319,7 @@ def plot_concentration_vs_mass(redshift):
     plt.xticks(np.logspace(3,15,5))
     plt.xlabel(r"$M_\mathrm{vir}$ ($M_\odot/h$)")
     plt.ylabel(r"Concentration")
-    plt.ylim(1e-8, 1)
+    plt.ylim(1,1e2)
     plt.legend(loc=0)
     plt.savefig("concentration.pdf")
     plt.clf()