diff --git a/.gitignore b/.gitignore index 894a44c..a9c339a 100644 --- a/.gitignore +++ b/.gitignore @@ -102,3 +102,7 @@ venv.bak/ # mypy .mypy_cache/ +.DS_Store +concentration.pdf +halomergerrate.pdf +volumemergerrate.pdf diff --git a/HW2/Schwarzchild_radius.py b/HW2/Schwarzchild_radius.py new file mode 100644 index 0000000..881b0c8 --- /dev/null +++ b/HW2/Schwarzchild_radius.py @@ -0,0 +1,12 @@ +import pint + +def calculate_schwarzchild_radius(mass): + ureg = pint.UnitRegistry() + ureg.define("Msun = 1.98855*10**30 * kg") + + schwarzchild_radius = 2 * ureg.newtonian_constant_of_gravitation * ureg.Msun * mass / ureg.c**2 + return schwarzchild_radius.to_base_units() + +if __name__ == "__main__": + mass = 1 # in solar mass unit + print(calculate_schwarzchild_radius(mass)) \ No newline at end of file diff --git a/HW2/matrix_multiplication.py b/HW2/matrix_multiplication.py new file mode 100644 index 0000000..53b21b9 --- /dev/null +++ b/HW2/matrix_multiplication.py @@ -0,0 +1,51 @@ +import time + +import numpy as np + +np.random.seed(2046) + + + +def loop_matrix_multiply(m1, m2): + result = np.zeros((m1.shape[0], m2.shape[1])) + for i in range(m1.shape[0]): + for j in range(m2.shape[1]): + for k in range(m1.shape[1]): + result[i, j] += m1[i, k] * m2[k, j] + return result + +def test_loop(m1, m2): + start = time.time() + assert np.allclose(loop_matrix_multiply(m1, m2), np.eye(m1.shape[0])) + end = time.time() + + print(f"Loop Method: {end - start:.6f} seconds") + +def list_matrix_multiply(m1, m2): + return np.array([[sum(m1[i, :] * m2[:, j]) for j in range(m2.shape[1])] for i in range(m1.shape[0])]) + +def test_list(m1, m2): + start = time.time() + assert np.allclose(list_matrix_multiply(m1, m2), np.eye(m1.shape[0])) + end = time.time() + + print(f"List Method: {end - start:.6f} seconds") + +def numpy_matrix_multiply(m1, m2): + return np.dot(m1, m2) + +def test_numpy(m1, m2): + start = time.time() + assert np.allclose(numpy_matrix_multiply(m1, m2), np.eye(m1.shape[0])) + end = time.time() + + print(f"Numpy Method: {end - start:.6f} seconds") + +if __name__ == "__main__": + m1 = np.random.rand(500, 500) + m2 = np.linalg.inv(m1) + + test_loop(m1, m2) + test_list(m1, m2) + test_numpy(m1, m2) + diff --git a/pbhmergers.py b/pbhmergers.py index 95b3d81..c72dbb6 100644 --- a/pbhmergers.py +++ b/pbhmergers.py @@ -25,9 +25,11 @@ def __init__(self,*args,conc_model="ludlow", conc_value=1.,hubble=0.67, **kwargs #Hubble constant and related objects! #The try except is because the method names change between pint 0.6 and pint 0.8 try: - self.ureg.define(pint.unit.UnitDefinition('hub', '', (),pint.unit.ScaleConverter(hubble))) - except AttributeError: + # self.ureg.define(pint.unit.UnitDefinition('hub', '', (),pint.unit.ScaleConverter(hubble))) + # except AttributeError: self.ureg.define(pint.definitions.UnitDefinition('hub', '', (),pint.converters.ScaleConverter(hubble))) + except AttributeError: + self.ureg.define('hub = %g' % hubble) self.ureg.define("Msolarh = Msolar / hub") self.ureg.define("Mpch = Mpc / hub") #Factor of R_s at which the maximum circular velocity of the halo is reached. @@ -63,7 +65,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): @@ -116,7 +118,7 @@ def cross_section(self, mass): gammaint = sigma**(10/7)*scipy.special.gammainc(5/7,vbysig**2)* scipy.special.gamma(5/7)/2 #We also need to normalise the probability function for v: #Integrate[4*Pi*v^2*P[v, sigma, Vvir], {v, 0, Vvir}] - probnorm = math.pi**(3/2)*sigma**3*scipy.special.erf(vbysig) - 2*math.pi/3*np.exp(-(vvir**2/sigma**2))*(3*sigma**2*vvir + 2*vvir**3) *0 + probnorm = math.pi**(3/2)*sigma**3*scipy.special.erf(vbysig) - 2*math.pi/3*np.exp(-(vvir**2/sigma**2))*(3*sigma**2*vvir + 2*vvir**3) assert np.all(probnorm.magnitude > 0) cross_section = prefac*(gammaint + cutoff)/probnorm assert self.ureg.get_dimensionality('[length]**3 [time]**(-1) [mass]**(-2)') == self.ureg.get_dimensionality(cross_section) @@ -224,7 +226,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""" diff --git a/power_specs.py b/power_specs.py index 1e00346..571aabd 100644 --- a/power_specs.py +++ b/power_specs.py @@ -638,7 +638,7 @@ 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]) + # 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 +660,12 @@ 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)) + # 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 +692,6 @@ 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)) + # 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))