diff --git a/examples/CoolingHaloWithSpin/internal_energy_profile.py b/examples/CoolingHaloWithSpin/internal_energy_profile.py index 854bdf223cfae75203a1924b4af6136b4b7aa6cd..a3e470cc24a939c9bc915371e927d9bd39196bff 100644 --- a/examples/CoolingHaloWithSpin/internal_energy_profile.py +++ b/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() diff --git a/examples/CoolingHaloWithSpin/radial_profile.py b/examples/CoolingHaloWithSpin/radial_profile.py index 335f7089b6835b65cf37e1bcd312a17966c295a7..ea282328e5b75530a128eab2dec5f065e46cf819 100644 --- a/examples/CoolingHaloWithSpin/radial_profile.py +++ b/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() diff --git a/examples/HydrostaticHalo/internal_energy_profile.py b/examples/HydrostaticHalo/internal_energy_profile.py index 854bdf223cfae75203a1924b4af6136b4b7aa6cd..65ba231ed7732bd386f663a364c053b5ec9a124e 100644 --- a/examples/HydrostaticHalo/internal_energy_profile.py +++ b/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() diff --git a/examples/HydrostaticHalo/radial_profile.py b/examples/HydrostaticHalo/radial_profile.py index 335f7089b6835b65cf37e1bcd312a17966c295a7..fd37df88048fbab997d6f3893254371c1c9ae022 100644 --- a/examples/HydrostaticHalo/radial_profile.py +++ b/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()