diff --git a/examples/HydrostaticHalo/radial_profile.py b/examples/HydrostaticHalo/radial_profile.py new file mode 100644 index 0000000000000000000000000000000000000000..efb4a155c07d85b64fe709cefe3d788c13b34efa --- /dev/null +++ b/examples/HydrostaticHalo/radial_profile.py @@ -0,0 +1,92 @@ +import numpy as np +import h5py as h5 +import matplotlib.pyplot as plt +import sys + +n_snaps = 11 + +#for the plotting +n_radial_bins = int(sys.argv[1]) + +#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 +h = 0.67777 # hubble parameter +gamma = 5./3. +eta = 1.2349 +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" +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_cgs = v_c * unit_velocity_cgs +header = f["Header"] +N = header.attrs["NumPart_Total"][0] +box_centre = np.array(header.attrs["BoxSize"]) + +#calculate r_vir and M_vir from v_c +r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA)) +M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS + +for i in range(n_snaps): + + filename = "Hydrostatic_%03d.hdf5" %i + f = h5.File(filename,'r') + coords_dset = f["PartType0/Coordinates"] + coords = np.array(coords_dset) +#translate coords by centre of box + header = f["Header"] + snap_time = header.attrs["Time"] + snap_time_cgs = snap_time * unit_time_cgs + coords[:,0] -= box_centre[0]/2. + coords[:,1] -= box_centre[1]/2. + coords[:,2] -= box_centre[2]/2. + radius = np.sqrt(coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2) + radius_cgs = radius*unit_length_cgs + radius_over_virial_radius = radius_cgs / r_vir_cgs + + 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 + +#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_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 + +#now divide hist by the volume so we have a density in each bin + + density = mass_hist / volume + + t = np.linspace(0.01,1.0,1000) + rho_analytic = t**(-2)/(4.*np.pi) + + plt.plot(radial_bin_mids,density,'ko',label = "Numerical solution") + plt.plot(t,rho_analytic,label = "Analytic Solution") + plt.legend(loc = "upper right") + plt.xlabel(r"$r / r_{vir}$") + plt.ylabel(r"$\rho / (M_{vir} / r_{vir}^3)$") + 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)) + plot_filename = "density_profile_%03d.png" %i + plt.savefig(plot_filename,format = "png") + plt.close() +