Page MenuHomec4science
Files F85249193

test.py
No OneTemporary
Actions

File Metadata

Created
Fri, Sep 27, 19:19

test.py
View Options

#!/usr/bin/env python
from PyGrackle import grackle
from numpy import *
import sys
params = {}
# grackle parameters
params['GrackleCloudyTable'] = "CloudyData_UVB=HM2012.h5"
#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
grackle.init_grackle(params)
#params = grackle.GetParameters()
# some units stuffs
mass_h = 1.67262171e-24
kboltz = 1.3806504e-16
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;
kboltz_LU = kboltz/units_energy ;
mass_h_LU = mass_h/units_mass ;
units_temperature = mass_h_LU / kboltz_LU;
# compute the cooling function
Zsol = 0.02041;
energy = 2.199e-08 *1e4; # 1K
density = 0.000110712362196 *100; # 1 atom/cc
dtime = 0.01; # here in Myrs */
ne = 0.0; # mass fraction of eletron */ /* useless for GRACKLE_CHEMISTRY = 0 */
Z = Zsol; # mass fraction of metals */
redshift = 0.0;
gamma = 5/3.;
xi = 0.76;
mu = 4./(8.-5.*(1.-xi));
nH = xi*density/mass_h_LU;
a_now = 1. / (1. + redshift);
n = 10001
Ts = zeros(n,float)
Ls = zeros(n,float)
for i in xrange(n):
T = pow(10,i/1000.);
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([1e1,1e9,1e-29,1e-21])
pyplot.show()

Event Timeline

c4science · Help