Homec4science
Diffusion Gear (master) Browse

Gear/PyGrackle/scripts/grackle_cooling_function.py
a136c2e88a2amaster

Last Change
Raw File

grackle_cooling_function.py
View Options

#!/usr/bin/env python
from PyGrackle import grackle
from numpy import *
import sys
from pNbody import ctes
from pNbody import units
from pNbody import thermodyn
from optparse import OptionParser
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser.add_option("--CoolingFile",
action="store",
dest="CoolingFile",
type="string",
default = "CloudyData_UVB=HM2012.h5",
help="cloudy cooling file",
metavar=" FLOAT")
parser.add_option("--Density",
action="store",
dest="Density",
type="float",
default = 100,
help="Density [a/cc]",
metavar=" FLOAT")
parser.add_option("--Z",
action="store",
dest="Z",
type="float",
default = 1.,
help="Metallicity in Solar units",
metavar=" FLOAT")
parser.add_option("--Redshift",
action="store",
dest="Redshift",
type="float",
default = 0.0,
help="Redshift",
metavar=" FLOAT")
parser.add_option("--UV",
action="store",
dest="UV",
type="int",
default = 1,
help="UV on/off",
metavar=" FLOAT")
(options, args) = parser.parse_args()
#pt.check_files_number(args)
files = args
return files,options
########################################################################
# MAIN
########################################################################
files,opt = parse_options()
params = {}
# grackle parameters
params['GrackleCloudyTable'] = opt.CoolingFile
params['GrackleUVbackground'] = opt.UV
#units parameters
params['UnitLength_in_cm'] = 3.085e+21
params['UnitMass_in_g'] = 1.989e+43
params['UnitVelocity_in_cm_per_s'] = 1e5
params['HubbleParam'] = 1.0
'''
params['UnitLength_in_cm'] = 3.085e+21;
params['UnitMass_in_g'] = 4.435693e+44;
params['UnitVelocity_in_cm_per_s'] = 97824708.2699;
params['HubbleParam'] = 1.0
'''
# dsph units
params['UnitLength_in_cm'] = 3.085e+21
params['UnitMass_in_g'] = 1.989e+43
params['UnitVelocity_in_cm_per_s'] = 20725573.785998672
params['HubbleParam'] = 1.0
# init grackle
grackle.init_grackle(params)
#params = grackle.GetParameters()
# some units stuffs
UnitMass_in_g = params['UnitMass_in_g']
UnitLength_in_cm = params['UnitLength_in_cm']
UnitVelocity_in_cm_per_s = params['UnitVelocity_in_cm_per_s']
HubbleParam = params['HubbleParam']
UnitDensity_in_cgs=UnitMass_in_g/pow(UnitLength_in_cm,3);
UnitTime_in_s=UnitLength_in_cm /UnitVelocity_in_cm_per_s;
units_density = UnitDensity_in_cgs * HubbleParam * HubbleParam;
units_length = UnitLength_in_cm / HubbleParam;
units_time = UnitTime_in_s / HubbleParam;
units_mass = units_density*pow(units_length,3);
units_velocity = units_length / units_time;
units_energy = units_velocity * units_velocity * units_mass;
code_units = units.Set_SystemUnits_From_Params(params)
Unit_atom = ctes.PROTONMASS.into(units.cgs)*units.Unit_g
Unit_atom.set_symbol('atom')
out_units = units.UnitSystem('local',[units.Unit_cm,Unit_atom,units.Unit_s,units.Unit_K])
atomCCtoDensityCU = out_units.convertionFactorTo(code_units.UnitDensity)
mass_h_LU = ctes.PROTONMASS.into(code_units)
kboltz_LU = ctes.BOLTZMANN.into(code_units)
units_temperature = mass_h_LU / kboltz_LU;
# compute the cooling function
Zsol = 0.02041;
density = atomCCtoDensityCU *opt.Density; # 1 atom/cc
ne = 0.0; # mass fraction of eletron */ /* useless for GRACKLE_CHEMISTRY = 0 */
Z = Zsol*opt.Z; # mass fraction of metals */
redshift = opt.Redshift;
gamma = 5/3.;
xi = 0.76;
mu = 4./(8.-5.*(1.-xi));
nH = xi*density/mass_h_LU;
a_now = 1. / (1. + redshift);
thermopars = {"k":kboltz_LU,"mh":mass_h_LU,"mu":mu,"gamma":gamma}
n = 10001
Ts = zeros(n,float)
Ls = zeros(n,float)
for i in xrange(n):
T = pow(10,i/1000.);
energy = thermodyn.Urt(None,T,thermopars) # here, mu is correctly computed
#energy = T/(gamma-1)/mu / units_temperature;
cooling_time = grackle.GetCoolingTime(energy, density, ne, Z, a_now);
dudt = energy/cooling_time;
Lambda_n = dudt * density /pow(nH,2) * units_energy*pow(units_length,3)/units_time ;
Ts[i] = T
Ls[i] = Lambda_n
#####################################################
# plot
#####################################################
from matplotlib import pyplot
pyplot.loglog(Ts,Ls)
pyplot.xlabel('T [K]')
pyplot.ylabel('$\\Lambda$ [erg s$^{-1}$ cm$^{3}$]')
pyplot.axis([1e2,1e9,1e-26,1e-21])
pyplot.show()
c4science · Help