Homec4science
Diffusion Gear (master) Browse

Gear/PyCool/scripts/pycool_cooling_wiersma_plot_CoolingTime_vs_TemperatureDensity
b87f34a0ce7emaster

Last Change
Raw File

pycool_cooling_wiersma_plot_CoolingTime_vs_TemperatureDensity
View Options

#!/usr/bin/env python
import PyCool.CoolingWiersma as CoolingTables
from numpy import *
import sys
from pNbody import ctes
from pNbody import units
from pNbody import thermodyn
from optparse import OptionParser
import Ptools as pt
ZHsol = 0.02 # (MZ/MH)_sun
FeHsol = 0.00176604 # (MFe/MH)_sun
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_postscript_options(parser)
parser.add_option("--xmin",
action="store",
dest="xmin",
type="float",
default=1e1,
help="min value in x")
parser.add_option("--xmax",
action="store",
dest="xmax",
type="float",
default=1e9,
help="max value in x")
parser.add_option("--ymin",
action="store",
dest="ymin",
type="float",
default=1e-5,
help="min value in y")
parser.add_option("--ymax",
action="store",
dest="ymax",
type="float",
default=1e5,
help="max value in y")
parser.add_option("--log",
action="store",
dest="log",
type="string",
default='xy',
help="Use log to plot values (x,y,xy)")
parser.add_option("--legend",
action="store_true",
dest="legend",
default = False,
help="add a legend")
parser.add_option("--FeHSolar",
action="store",
dest="FeHSolar",
type="float",
default = FeHsol,
help="Fe solar abundance",
metavar=" FLOAT")
parser.add_option("--ZHSolar",
action="store",
dest="ZHSolar",
type="float",
default = ZHsol,
help="Z solar abundance",
metavar=" FLOAT")
parser.add_option("--Z",
action="store",
dest="Z",
type="float",
default = ZHsol,
help="Metallicity",
metavar=" FLOAT")
parser.add_option("--Redshift",
action="store",
dest="Redshift",
type="float",
default = 0.0,
help="Redshift",
metavar=" FLOAT")
parser.add_option("--Temperature",
action="store",
dest="Temperature",
type="float",
default = 1e5,
help="Temperature",
metavar=" FLOAT")
parser.add_option("--Density",
action="store",
dest="Density",
type="float",
default = 100,
help="Density [a/cc]",
metavar=" FLOAT")
parser.add_option("--HeMfaction",
action="store",
dest="HeMfaction",
type="float",
default = 0.2806,
help="Helium mass fraction",
metavar=" FLOAT")
parser.add_option("--UV",
action="store",
dest="UV",
type="int",
default = 1,
help="UV on/off",
metavar=" FLOAT")
parser.add_option("--coolingdir",
action="store",
dest="coolingdir",
default = "/home/epfl/revaz/code/gear/PyCool/tables_wiersma/coolingtables/",
help="cooling directory")
(options, args) = parser.parse_args()
#pt.check_files_number(args)
files = args
return files,options
########################################################################
# MAIN
########################################################################
files,opt = parse_options()
files,opt = parse_options()
# init cooling
CoolingTables.Init(opt.coolingdir)
CoolingTables.setTablesFromRedshift(opt.Redshift)
# define MFe/MH from Z (simple scaling)
Metal = opt.Z/ZHsol * FeHsol
print "[Fe/H] = %g"%(log10((Metal/FeHsol)))
# initialize temperature and coolingfct
Ts = 10**arange(1,9,0.02)
nHs = 10**arange(-7,2,0.02)
Cs = zeros( (len(Ts),len(nHs)) )
# unit stuffs
cgs_units = units.UnitSystem('cgs',[units.Unit_cm,units.Unit_g,units.Unit_s,units.Unit_K])
Unit_atom = ctes.PROTONMASS.into(units.cgs)*units.Unit_g
density_units = units.UnitSystem('local',[units.Unit_cm,Unit_atom,units.Unit_Myr,units.Unit_K])
time_units = units.UnitSystem('local',[units.Unit_cm,Unit_atom,units.Unit_Myr,units.Unit_K])
# get k (boltzmann) in cgs_units
k = ctes.BOLTZMANN.into(cgs_units)
# get mh (hydrogen mass) in cgs_units
mh = ctes.PROTONMASS.into(cgs_units)
Xi = 0.76
gamma = 5/3.
pars = {}
pars['k'] = k
pars['mh'] = mh
pars['gamma'] = gamma
pars['mu'] = 1 # unused here
nH = opt.Density
for i,T in enumerate(Ts):
for j,nH in enumerate(nHs):
L = CoolingTables.compute_LambdaTotal(T,nH,opt.HeMfaction,opt.Z,ZHsol)
rho = nH * density_units.convertionFactorTo(cgs_units.UnitDensity) # to gr/cm3
U = thermodyn.Urt(rho,T,pars)
Tcool = mh**2 / (fabs(L*Xi**2)) * U/rho
Tcool = Tcool * cgs_units.convertionFactorTo(time_units.UnitTime)
Cs[i,j] = Tcool
########################################################################
# plot
########################################################################
xmin=min(log10(nHs))
xmax=max(log10(nHs))
ymin=min(log10(Ts))
ymax=max(log10(Ts))
vmin = -5
vmax = 5
data = log10(fabs(Cs))
cmap = pt.GetColormap('rainbow4',revesed=False)
im = pt.imshow(data, interpolation='bilinear', origin='lower',cmap=cmap, extent=(xmin,xmax,ymin,ymax),aspect='auto',vmin=vmin,vmax=vmax)
pt.xlabel(r"$\rm{log10}\,\,\rm{Density}\,\rm{[a/cm^3]}$")
pt.ylabel(r"$\rm{log10}\,\,\rm{Temperature}\,\rm{[K]}$")
cb = pt.colorbar()
cb.set_label('log10 Cooling\,Time [Myr]')
pt.show()
c4science · Help