Commit 969f0a24 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'Stefan_add_README_files' into 'stats_include_e_pot_ext'

Stefan add readme files



See merge request !274
parents 6578b2fa c03fb411
;
;To make the initial conditions we distribute gas particles randomly in a cube
;with a side length twice that of the virial radius. The density profile
;of the gas is proportional to r^(-2) where r is the distance from the centre
;of the cube.
;
;The parameter v_rot (in makeIC.py and cooling.yml) sets the circular velocity
;of the halo, and by extension, the viral radius, viral mass, and the
;internal energy of the gas such that hydrostatic equilibrium is achieved.
;
;While the system is initially in hydrostatic equilibrium, the cooling of the
;gas means that the halo will collapse.
;
;To run this example, make such that the code is compiled with either the
;isothermal potential or softened isothermal potential, and 'const_lambda'
;cooling, set in src/const.h. In
;the latter case, a (small) value of epsilon needs to be set in cooling.yml.
;0.1 kpc should work well.
;
;The plotting scripts produce a plot of the density, internal energy and radial
;velocity profile for each snapshot. test_energy_conservation.py shows the
;evolution of energy with time. These can be used to check if the example
;has run properly.
;
;
\ No newline at end of file
...@@ -6,8 +6,8 @@ python makeIC.py 10000 ...@@ -6,8 +6,8 @@ python makeIC.py 10000
../swift -g -s -C -t 16 cooling_halo.yml 2>&1 | tee output.log ../swift -g -s -C -t 16 cooling_halo.yml 2>&1 | tee output.log
# python radial_profile.py 10 python radial_profile.py 2. 200 100
# python internal_energy_profile.py 10 python internal_energy_profile.py 2. 200 100
# python test_energy_conservation.py python test_energy_conservation.py 2. 200 100
import numpy as np
import h5py as h5
import matplotlib.pyplot as plt
import sys
def do_binning(x,y,x_bin_edges):
#x and y are arrays, where y = f(x)
#returns number of elements of x in each bin, and the total of the y elements corresponding to those x values
n_bins = x_bin_edges.size - 1
count = np.zeros(n_bins)
y_totals = np.zeros(n_bins)
for i in range(n_bins):
ind = np.intersect1d(np.where(x > bin_edges[i])[0],np.where(x <= bin_edges[i+1])[0])
count[i] = ind.size
binned_y = y[ind]
y_totals[i] = np.sum(binned_y)
return(count,y_totals)
#for the plotting
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
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 = "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["SoftenedIsothermalPotential: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 = "CoolingHalo_%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
#get the internal energies
vel_dset = f["PartType0/Velocities"]
vel = np.array(vel_dset)
#make dimensionless
vel /= v_c
r = radius_over_virial_radius
#find radial component of velocity
v_r = np.zeros(r.size)
for j in range(r.size):
v_r[j] = -np.dot(coords[j,:],vel[j,:])/radius[j]
bin_edges = np.linspace(0,max_r,n_radial_bins + 1)
(hist,v_r_totals) = do_binning(r,v_r,bin_edges)
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_v_r = v_r_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_v_r,'ko',label = "Average radial velocity in shell")
#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 = "upper right")
plt.xlabel(r"$r / r_{vir}$")
plt.ylabel(r"$v_r / v_c$")
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,2))
plot_filename = "./plots/radial_velocity_profile/velocity_profile_%03d.png" %i
plt.savefig(plot_filename,format = "png")
plt.close()
;
;To make the initial conditions we distribute gas particles randomly in a cube
;with a side length twice that of the virial radius. The density profile
;of the gas is proportional to r^(-2) where r is the distance from the centre
;of the cube.
;
;The parameter v_rot (in makeIC.py and cooling.yml) sets the circular velocity
;of the halo, and by extension, the viral radius, viral mass, and the
;internal energy of the gas such that hydrostatic equilibrium is achieved.
;
;The gas is given some angular momentum about the z-axis. This is defined by
;the 'spin_lambda' variable in makeIC.py
;
;While the system is initially in hydrostatic equilibrium, the cooling of the
;gas and the non-zero angular momentum means that the halo will collapse into
;a spinning disc.
;
;To run this example, make such that the code is compiled with either the
;isothermal potential or softened isothermal potential, and 'const_lambda'
;cooling, set in src/const.h. In
;the latter case, a (small) value of epsilon needs to be set in cooling.yml.
;0.1 kpc should work well.
;
;The plotting scripts produce a plot of the density, internal energy and radial
;velocity profile for each snapshot. test_energy_conservation.py shows the
;evolution of energy with time. These can be used to check if the example
;has run properly.
;
;
\ No newline at end of file
import numpy as np
import h5py as h5
import matplotlib.pyplot as plt
import sys
def do_binning(x,y,x_bin_edges):
#x and y are arrays, where y = f(x)
#returns number of elements of x in each bin, and the total of the y elements corresponding to those x values
n_bins = x_bin_edges.size - 1
count = np.zeros(n_bins)
y_totals = np.zeros(n_bins)
for i in range(n_bins):
ind = np.intersect1d(np.where(x > bin_edges[i])[0],np.where(x <= bin_edges[i+1])[0])
count[i] = ind.size
binned_y = y[ind]
y_totals[i] = np.sum(binned_y)
return(count,y_totals)
#for the plotting
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
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 = "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["SoftenedIsothermalPotential: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 = "CoolingHalo_%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
#get the internal energies
vel_dset = f["PartType0/Velocities"]
vel = np.array(vel_dset)
#make dimensionless
vel /= v_c
r = radius_over_virial_radius
#find radial component of velocity
v_r = np.zeros(r.size)
for j in range(r.size):
v_r[j] = -np.dot(coords[j,:],vel[j,:])/radius[j]
bin_edges = np.linspace(0,max_r,n_radial_bins + 1)
(hist,v_r_totals) = do_binning(r,v_r,bin_edges)
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_v_r = v_r_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_v_r,'ko',label = "Average radial velocity in shell")
#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 = "upper right")
plt.xlabel(r"$r / r_{vir}$")
plt.ylabel(r"$v_r / v_c$")
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,2))
plot_filename = "./plots/radial_velocity_profile/velocity_profile_%03d.png" %i
plt.savefig(plot_filename,format = "png")
plt.close()
; ;
;In this example we distribute a set of gas particles in sphere with an r^-2 ;To make the initial conditions we distribute gas particles randomly in a cube
;density profile, with an external isothermal potential. The particles are ;with a side length twice that of the virial radius. The density profile
;stationary with a pressure gradient such that we have hydrostatic equilibrium. ;of the gas is proportional to r^(-2) where r is the distance from the centre
;of the cube.
; ;
;The mass and radius of the halo is set by 'v_rot' ;The parameter v_rot (in makeIC.py and cooling.yml) sets the circular velocity
\ No newline at end of file ;of the halo, and by extension, the viral radius, viral mass, and the
;internal energy of the gas such that hydrostatic equilibrium is achieved.
;
;To run this example, make such that the code is compiled with either the
;isothermal potential or softened isothermal potential set in src/const.h. In
;the latter case, a (small) value of epsilon needs to be set in cooling.yml.
;0.1 kpc should work well.
;
;The plotting scripts produce a plot of the density, internal energy and radial
;velocity profile for each snapshot. test_energy_conservation.py shows the
;evolution of energy with time. These can be used to check if the example
;has run properly.
;
\ No newline at end of file
...@@ -28,7 +28,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"]) ...@@ -28,7 +28,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"]) unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"]) unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
unit_time_cgs = unit_length_cgs / unit_velocity_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 v_c_cgs = v_c * unit_velocity_cgs
#lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"]) #lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
#X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"]) #X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
...@@ -114,7 +114,7 @@ for i in range(n_snaps): ...@@ -114,7 +114,7 @@ for i in range(n_snaps):
#plt.xscale('log') #plt.xscale('log')
#plt.yscale('log') #plt.yscale('log')
plt.legend(loc = "upper right") plt.legend(loc = "upper right")
plot_filename = "density_profile_%03d.png" %i plot_filename = "./plots/density_profile/density_profile_%03d.png" %i
plt.savefig(plot_filename,format = "png") plt.savefig(plot_filename,format = "png")
plt.close() plt.close()
...@@ -9,7 +9,7 @@ InternalUnitSystem: ...@@ -9,7 +9,7 @@ InternalUnitSystem:
# Parameters governing the time integration # Parameters governing the time integration
TimeIntegration: TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units). time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 10.0 # The end time of the simulation (in internal units). time_end: 1.0 # The end time of the simulation (in internal units).
dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units). dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units). dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
...@@ -38,10 +38,11 @@ InitialConditions: ...@@ -38,10 +38,11 @@ InitialConditions:
shift_z: 0. shift_z: 0.
# External potential parameters # External potential parameters
IsothermalPotential: SoftenedIsothermalPotential:
position_x: 0. # location of centre of isothermal potential in internal units position_x: 0. # location of centre of isothermal potential in internal units
position_y: 0. position_y: 0.
position_z: 0. position_z: 0.
vrot: 200. # rotation speed of isothermal potential in internal units vrot: 200. # rotation speed of isothermal potential in internal units
epsilon: 0.1
timestep_mult: 0.03 # controls time step timestep_mult: 0.03 # controls time step
...@@ -46,7 +46,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"]) ...@@ -46,7 +46,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"]) unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"]) unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
unit_time_cgs = unit_length_cgs / unit_velocity_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 v_c_cgs = v_c * unit_velocity_cgs
#lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"]) #lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
#X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"]) #X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
...@@ -102,7 +102,7 @@ for i in range(n_snaps): ...@@ -102,7 +102,7 @@ for i in range(n_snaps):
plt.ylabel(r"$u / (v_c^2 / (2(\gamma - 1)) $") 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.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
plt.ylim((0,2)) plt.ylim((0,2))
plot_filename = "internal_energy_profile_%03d.png" %i plot_filename = "./plots/internal_energy/internal_energy_profile_%03d.png" %i
plt.savefig(plot_filename,format = "png") plt.savefig(plot_filename,format = "png")
plt.close() plt.close()
......
...@@ -6,8 +6,10 @@ python makeIC.py 10000 ...@@ -6,8 +6,10 @@ python makeIC.py 10000
../swift -g -s -t 16 hydrostatic.yml 2>&1 | tee output.log ../swift -g -s -t 16 hydrostatic.yml 2>&1 | tee output.log
python radial_profile.py 10 python density_profile.py 2. 200 100
python internal_energy_profile.py 10 python internal_energy_profile.py 2. 200 100
python velocity_profile.py 2. 200 100
python test_energy_conservation.py python test_energy_conservation.py
...@@ -91,6 +91,5 @@ plt.xlabel(r"$t / t_{dyn}$") ...@@ -91,6 +91,5 @@ plt.xlabel(r"$t / t_{dyn}$")
plt.ylabel(r"$E / v_c^2$") plt.ylabel(r"$E / v_c^2$")
plt.title(r"$%d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(N,v_c)) plt.title(r"$%d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(N,v_c))
plt.ylim((-2,2)) plt.ylim((-2,2))
#plot_filename = "density_profile_%03d.png" %i plt.savefig("energy_conservation.png",format = 'png')
plt.show()
import numpy as np
import h5py as h5
import matplotlib.pyplot as plt
import sys
def do_binning(x,y,x_bin_edges):
#x and y are arrays, where y = f(x)
#returns number of elements of x in each bin, and the total of the y elements corresponding to those x values
n_bins = x_bin_edges.size - 1
count = np.zeros(n_bins)
y_totals = np.zeros(n_bins)
for i in range(n_bins):
ind = np.intersect1d(np.where(x > bin_edges[i])[0],np.where(x <= bin_edges[i+1])[0])
count[i] = ind.size
binned_y = y[ind]
y_totals[i] = np.sum(binned_y)
return(count,y_totals)
#for the plotting
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
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["SoftenedIsothermalPotential: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
#get the internal energies
vel_dset = f["PartType0/Velocities"]
vel = np.array(vel_dset)
#make dimensionless
vel /= v_c
r = radius_over_virial_radius
#find radial component of velocity
v_r = np.zeros(r.size)
for j in range(r.size):
v_r[j] = -np.dot(coords[j,:],vel[j,:])/radius[j]
bin_edges = np.linspace(0,max_r,n_radial_bins + 1)
(hist,v_r_totals) = do_binning(r,v_r,bin_edges)
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_v_r = v_r_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_v_r,'ko',label = "Average radial velocity in shell")
#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 = "upper right")
plt.xlabel(r"$r / r_{vir}$")
plt.ylabel(r"$v_r / v_c$")
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,2))
plot_filename = "./plots/radial_velocity_profile/velocity_profile_%03d.png" %i
plt.savefig(plot_filename,format = "png")
plt.close()
...@@ -164,7 +164,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( ...@@ -164,7 +164,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
__attribute__((always_inline)) INLINE static float __attribute__((always_inline)) INLINE static float
external_gravity_get_potential_energy( external_gravity_get_potential_energy(
const struct external_potential* potential, const struct external_potential* potential,
const struct phys_const* const phys_const, const struct part* p) { const struct phys_const* const phys_const, const struct gpart* p) {
return 0.f; return 0.f;
} }
......
...@@ -78,7 +78,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( ...@@ -78,7 +78,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
__attribute__((always_inline)) INLINE static float __attribute__((always_inline)) INLINE static float
external_gravity_get_potential_energy( external_gravity_get_potential_energy(
const struct external_potential* potential, const struct external_potential* potential,
const struct phys_const* const phys_const, const struct part* g) { const struct phys_const* const phys_const, const struct gpart* g) {
return 0.f; return 0.f;
} }
......
...@@ -127,7 +127,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( ...@@ -127,7 +127,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
__attribute__((always_inline)) INLINE static float __attribute__((always_inline)) INLINE static float
external_gravity_get_potential_energy( external_gravity_get_potential_energy(
const struct external_potential* potential, const struct external_potential* potential,
const struct phys_const* const phys_const, const struct part* g) { const struct phys_const* const phys_const, const struct gpart* g) {