Updated internal energy profile plotting script

examples/CoolingHaloWithSpin/internal_energy_profile.py

@@ -21,16 +21,17 @@ def do_binning(x,y,x_bin_edges):
     return(count,y_totals)
 
 
-n_snaps = 100
-
 #for the plotting
-n_radial_bins = int(sys.argv[1])
+max_r = float(sys.argv[1])
+n_radial_bins = int(sys.argv[2])
+n_snaps = int(sys.argv[3])
 
 #some constants
 OMEGA = 0.3 # Cosmological matter fraction at z = 0
 PARSEC_IN_CGS = 3.0856776e18
 KM_PER_SEC_IN_CGS = 1.0e5
 CONST_G_CGS = 6.672e-8
+CONST_m_H_CGS = 1.67e-24
 h = 0.67777 # hubble parameter
 gamma = 5./3.
 eta = 1.2349
@@ -38,15 +39,17 @@ H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
 
 #read some header/parameter information from the first snapshot
 
-filename = "Hydrostatic_000.hdf5"
+filename = "CoolingHalo_000.hdf5"
 f = h5.File(filename,'r')
 params = f["Parameters"]
 unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["IsothermalPotential:vrot"])
+v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
+lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
+X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
 header = f["Header"]
 N = header.attrs["NumPart_Total"][0]
 box_centre = np.array(header.attrs["BoxSize"])
@@ -57,7 +60,7 @@ M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
 
 for i in range(n_snaps):
 
-    filename = "Hydrostatic_%03d.hdf5" %i
+    filename = "CoolingHalo_%03d.hdf5" %i
     f = h5.File(filename,'r')
     coords_dset = f["PartType0/Coordinates"]
     coords = np.array(coords_dset)
@@ -80,21 +83,25 @@ for i in range(n_snaps):
     u /= v_c**2/(2. * (gamma - 1.))
     r = radius_over_virial_radius
 
-    bin_edges = np.linspace(0,1,n_radial_bins + 1)
+    bin_edges = np.linspace(0,max_r,n_radial_bins + 1)
     (hist,u_totals) = do_binning(r,u,bin_edges)
     
-    bin_widths = 1. / n_radial_bins
-    radial_bin_mids = np.linspace(bin_widths / 2. , 1. - bin_widths / 2. , n_radial_bins) 
+    bin_widths = bin_edges[1] - bin_edges[0]
+    radial_bin_mids = np.linspace(bin_widths / 2. , max_r - bin_widths / 2. , n_radial_bins) 
     binned_u = u_totals / hist
 
+    #calculate cooling radius
+
+    r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs)
 
     plt.plot(radial_bin_mids,binned_u,'ko',label = "Numerical solution")
-    plt.plot((0,1),(1,1),label = "Analytic Solution")
+    #plt.plot((0,1),(1,1),label = "Analytic Solution")
+    plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,2),'r',label = "Cooling radius")
     plt.legend(loc = "lower right")
     plt.xlabel(r"$r / r_{vir}$")
     plt.ylabel(r"$u / (v_c^2 / (2(\gamma - 1)) $")
     plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
-    plt.ylim((0,1))
+    plt.ylim((0,2))
     plot_filename = "internal_energy_profile_%03d.png" %i
     plt.savefig(plot_filename,format = "png")
     plt.close()

examples/CoolingHaloWithSpin/radial_profile.py

@@ -3,16 +3,17 @@ import h5py as h5
 import matplotlib.pyplot as plt
 import sys
 
-n_snaps = 11
-
 #for the plotting
-#n_radial_bins = int(sys.argv[1])
+max_r = float(sys.argv[1]) #in units of the virial radius 
+n_radial_bins = int(sys.argv[2])
+n_snaps = int(sys.argv[3])
 
 #some constants
 OMEGA = 0.3 # Cosmological matter fraction at z = 0
 PARSEC_IN_CGS = 3.0856776e18
 KM_PER_SEC_IN_CGS = 1.0e5
 CONST_G_CGS = 6.672e-8
+CONST_m_H_CGS = 1.67e-24
 h = 0.67777 # hubble parameter
 gamma = 5./3.
 eta = 1.2349
@@ -20,15 +21,17 @@ H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
 
 #read some header/parameter information from the first snapshot
 
-filename = "Hydrostatic_000.hdf5"
+filename = "CoolingHalo_000.hdf5"
 f = h5.File(filename,'r')
 params = f["Parameters"]
 unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["IsothermalPotential:vrot"])
+v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
+lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
+X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
 header = f["Header"]
 N = header.attrs["NumPart_Total"][0]
 box_centre = np.array(header.attrs["BoxSize"])
@@ -39,7 +42,7 @@ M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
 
 for i in range(n_snaps):
 
-    filename = "Hydrostatic_%03d.hdf5" %i
+    filename = "CoolingHalo_%03d.hdf5" %i
     f = h5.File(filename,'r')
     coords_dset = f["PartType0/Coordinates"]
     coords = np.array(coords_dset)
@@ -56,45 +59,61 @@ for i in range(n_snaps):
 
     r = radius_over_virial_radius
 
-    # bin_width = 1./n_radial_bins
-#     hist = np.histogram(r,bins = n_radial_bins)[0] # number of particles in each bin
+    bin_edges = np.linspace(0.,max_r,n_radial_bins+1)
+    bin_width = bin_edges[1] - bin_edges[0]
+    hist = np.histogram(r,bins = bin_edges)[0] # number of particles in each bin
 
-# #find the mass in each radial bin
+#find the mass in each radial bin
 
-#     mass_dset = f["PartType0/Masses"]
-# #mass of each particles should be equal
-#     part_mass = np.array(mass_dset)[0]
-#     part_mass_cgs = part_mass * unit_mass_cgs
-#     part_mass_over_virial_mass = part_mass_cgs / M_vir_cgs 
+    mass_dset = f["PartType0/Masses"]
+#mass of each particles should be equal
+    part_mass = np.array(mass_dset)[0]
+    part_mass_cgs = part_mass * unit_mass_cgs
+    part_mass_over_virial_mass = part_mass_cgs / M_vir_cgs 
 
-#     mass_hist = hist * part_mass_over_virial_mass
-#     radial_bin_mids = np.linspace(bin_width/2.,1 - bin_width/2.,n_radial_bins)
-# #volume in each radial bin
-#     volume = 4.*np.pi * radial_bin_mids**2 * bin_width
+    mass_hist = hist * part_mass_over_virial_mass
+    radial_bin_mids = np.linspace(bin_width/2.,max_r - bin_width/2.,n_radial_bins)
+#volume in each radial bin
+    volume = 4.*np.pi * radial_bin_mids**2 * bin_width
 
-# #now divide hist by the volume so we have a density in each bin
+#now divide hist by the volume so we have a density in each bin
 
-#     density = mass_hist / volume
+    density = mass_hist / volume
 
-    # read the densities
+    ##read the densities
 
-    density_dset = f["PartType0/Density"]
-    density = np.array(density_dset)
-    density_cgs = density * unit_mass_cgs / unit_length_cgs**3
-    rho = density_cgs * r_vir_cgs**3 / M_vir_cgs
+    # density_dset = f["PartType0/Density"]
+    # density = np.array(density_dset)
+    # density_cgs = density * unit_mass_cgs / unit_length_cgs**3
+    # rho = density_cgs * r_vir_cgs**3 / M_vir_cgs
 
-    t = np.linspace(0.01,2.0,1000)
+    t = np.linspace(10./n_radial_bins,10.0,1000)
     rho_analytic = t**(-2)/(4.*np.pi)
 
-    plt.plot(r,rho,'x',label = "Numerical solution")
-    plt.plot(t,rho_analytic,label = "Analytic Solution")
-    plt.legend(loc = "upper right")
+    #calculate cooling radius
+
+    r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs)
+
+    #initial analytic density profile
+    
+    if (i == 0):
+        r_0 = radial_bin_mids[0]
+        rho_0 = density[0]
+
+        rho_analytic_init = rho_0 * (radial_bin_mids/r_0)**(-2)
+    plt.plot(radial_bin_mids,density/rho_analytic_init,'ko',label = "Average density of shell")
+    #plt.plot(t,rho_analytic,label = "Initial analytic density profile"
     plt.xlabel(r"$r / r_{vir}$")
-    plt.ylabel(r"$\rho / (M_{vir} / r_{vir}^3)$")
+    plt.ylabel(r"$\rho / \rho_{init})$")
     plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
-    #plt.ylim((0.1,40))
-    plt.xscale('log')
-    plt.yscale('log')
+    #plt.ylim((1.e-2,1.e1))
+    plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,20),'r',label = "Cooling radius")
+    plt.xlim((radial_bin_mids[0],max_r))
+    plt.ylim((0,20))
+    plt.plot((0,max_r),(1,1))
+    #plt.xscale('log')
+    #plt.yscale('log')
+    plt.legend(loc = "upper right")
     plot_filename = "density_profile_%03d.png" %i
     plt.savefig(plot_filename,format = "png")
     plt.close()

examples/HydrostaticHalo/internal_energy_profile.py

@@ -21,16 +21,17 @@ def do_binning(x,y,x_bin_edges):
     return(count,y_totals)
 
 
-n_snaps = 100
-
 #for the plotting
-n_radial_bins = int(sys.argv[1])
+max_r = float(sys.argv[1])
+n_radial_bins = int(sys.argv[2])
+n_snaps = int(sys.argv[3])
 
 #some constants
 OMEGA = 0.3 # Cosmological matter fraction at z = 0
 PARSEC_IN_CGS = 3.0856776e18
 KM_PER_SEC_IN_CGS = 1.0e5
 CONST_G_CGS = 6.672e-8
+CONST_m_H_CGS = 1.67e-24
 h = 0.67777 # hubble parameter
 gamma = 5./3.
 eta = 1.2349
@@ -47,6 +48,8 @@ unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"]
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
 v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
+#lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
+#X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
 header = f["Header"]
 N = header.attrs["NumPart_Total"][0]
 box_centre = np.array(header.attrs["BoxSize"])
@@ -80,21 +83,25 @@ for i in range(n_snaps):
     u /= v_c**2/(2. * (gamma - 1.))
     r = radius_over_virial_radius
 
-    bin_edges = np.linspace(0,1,n_radial_bins + 1)
+    bin_edges = np.linspace(0,max_r,n_radial_bins + 1)
     (hist,u_totals) = do_binning(r,u,bin_edges)
     
-    bin_widths = 1. / n_radial_bins
-    radial_bin_mids = np.linspace(bin_widths / 2. , 1. - bin_widths / 2. , n_radial_bins) 
+    bin_widths = bin_edges[1] - bin_edges[0]
+    radial_bin_mids = np.linspace(bin_widths / 2. , max_r - bin_widths / 2. , n_radial_bins) 
     binned_u = u_totals / hist
 
+    #calculate cooling radius
+
+    #r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs)
 
     plt.plot(radial_bin_mids,binned_u,'ko',label = "Numerical solution")
-    plt.plot((0,1),(1,1),label = "Analytic Solution")
+    #plt.plot((0,1),(1,1),label = "Analytic Solution")
+    #plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,2),'r',label = "Cooling radius")
     plt.legend(loc = "lower right")
     plt.xlabel(r"$r / r_{vir}$")
     plt.ylabel(r"$u / (v_c^2 / (2(\gamma - 1)) $")
     plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
-    plt.ylim((0,1))
+    plt.ylim((0,2))
     plot_filename = "internal_energy_profile_%03d.png" %i
     plt.savefig(plot_filename,format = "png")
     plt.close()

examples/HydrostaticHalo/radial_profile.py

@@ -3,16 +3,17 @@ import h5py as h5
 import matplotlib.pyplot as plt
 import sys
 
-n_snaps = 11
-
 #for the plotting
-#n_radial_bins = int(sys.argv[1])
+max_r = float(sys.argv[1]) #in units of the virial radius 
+n_radial_bins = int(sys.argv[2])
+n_snaps = int(sys.argv[3])
 
 #some constants
 OMEGA = 0.3 # Cosmological matter fraction at z = 0
 PARSEC_IN_CGS = 3.0856776e18
 KM_PER_SEC_IN_CGS = 1.0e5
 CONST_G_CGS = 6.672e-8
+CONST_m_H_CGS = 1.67e-24
 h = 0.67777 # hubble parameter
 gamma = 5./3.
 eta = 1.2349
@@ -29,6 +30,8 @@ unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"]
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
 v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
+#lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
+#X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
 header = f["Header"]
 N = header.attrs["NumPart_Total"][0]
 box_centre = np.array(header.attrs["BoxSize"])
@@ -56,45 +59,61 @@ for i in range(n_snaps):
 
     r = radius_over_virial_radius
 
-    # bin_width = 1./n_radial_bins
-#     hist = np.histogram(r,bins = n_radial_bins)[0] # number of particles in each bin
+    bin_edges = np.linspace(0.,max_r,n_radial_bins+1)
+    bin_width = bin_edges[1] - bin_edges[0]
+    hist = np.histogram(r,bins = bin_edges)[0] # number of particles in each bin
 
-# #find the mass in each radial bin
+#find the mass in each radial bin
 
-#     mass_dset = f["PartType0/Masses"]
-# #mass of each particles should be equal
-#     part_mass = np.array(mass_dset)[0]
-#     part_mass_cgs = part_mass * unit_mass_cgs
-#     part_mass_over_virial_mass = part_mass_cgs / M_vir_cgs 
+    mass_dset = f["PartType0/Masses"]
+#mass of each particles should be equal
+    part_mass = np.array(mass_dset)[0]
+    part_mass_cgs = part_mass * unit_mass_cgs
+    part_mass_over_virial_mass = part_mass_cgs / M_vir_cgs 
 
-#     mass_hist = hist * part_mass_over_virial_mass
-#     radial_bin_mids = np.linspace(bin_width/2.,1 - bin_width/2.,n_radial_bins)
-# #volume in each radial bin
-#     volume = 4.*np.pi * radial_bin_mids**2 * bin_width
+    mass_hist = hist * part_mass_over_virial_mass
+    radial_bin_mids = np.linspace(bin_width/2.,max_r - bin_width/2.,n_radial_bins)
+#volume in each radial bin
+    volume = 4.*np.pi * radial_bin_mids**2 * bin_width
 
-# #now divide hist by the volume so we have a density in each bin
+#now divide hist by the volume so we have a density in each bin
 
-#     density = mass_hist / volume
+    density = mass_hist / volume
 
-    # read the densities
+    ##read the densities
 
-    density_dset = f["PartType0/Density"]
-    density = np.array(density_dset)
-    density_cgs = density * unit_mass_cgs / unit_length_cgs**3
-    rho = density_cgs * r_vir_cgs**3 / M_vir_cgs
+    # density_dset = f["PartType0/Density"]
+    # density = np.array(density_dset)
+    # density_cgs = density * unit_mass_cgs / unit_length_cgs**3
+    # rho = density_cgs * r_vir_cgs**3 / M_vir_cgs
 
-    t = np.linspace(0.01,2.0,1000)
+    t = np.linspace(10./n_radial_bins,10.0,1000)
     rho_analytic = t**(-2)/(4.*np.pi)
 
-    plt.plot(r,rho,'x',label = "Numerical solution")
-    plt.plot(t,rho_analytic,label = "Analytic Solution")
-    plt.legend(loc = "upper right")
+    #calculate cooling radius
+
+    #r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs)
+
+    #initial analytic density profile
+    
+    if (i == 0):
+        r_0 = radial_bin_mids[0]
+        rho_0 = density[0]
+
+        rho_analytic_init = rho_0 * (radial_bin_mids/r_0)**(-2)
+    plt.plot(radial_bin_mids,density/rho_analytic_init,'ko',label = "Average density of shell")
+    #plt.plot(t,rho_analytic,label = "Initial analytic density profile"
     plt.xlabel(r"$r / r_{vir}$")
-    plt.ylabel(r"$\rho / (M_{vir} / r_{vir}^3)$")
+    plt.ylabel(r"$\rho / \rho_{init})$")
     plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
-    #plt.ylim((0.1,40))
-    plt.xscale('log')
-    plt.yscale('log')
+    #plt.ylim((1.e-2,1.e1))
+    #plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,20),'r',label = "Cooling radius")
+    plt.xlim((radial_bin_mids[0],max_r))
+    plt.ylim((0,20))
+    plt.plot((0,max_r),(1,1))
+    #plt.xscale('log')
+    #plt.yscale('log')
+    plt.legend(loc = "upper right")
     plot_filename = "density_profile_%03d.png" %i
     plt.savefig(plot_filename,format = "png")
     plt.close()