Homec4science
Diffusion Gear (master) Browse

Gear/PyCool/scripts/pycool_cooling_with_metals_plot_CoolingTime_vs_TemperatureDensity
5ec004b6ec8dmaster

Last Change
Raw File

pycool_cooling_with_metals_plot_CoolingTime_vs_TemperatureDensity
View Options

#!/usr/bin/env python
from PyCool import cooling
from numpy import *
import sys
from pNbody import units
from pNbody import ctes
from pNbody import thermodyn
from optparse import OptionParser
import Ptools as pt
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_postscript_options(parser)
parser = pt.add_limits_options(parser)
parser = pt.add_log_options(parser)
parser.add_option("--x",
action="store",
dest="x",
type="string",
default = 'r',
help="x value to plot",
metavar=" STRING")
parser.add_option("--y",
action="store",
dest="y",
type="string",
default = 'T',
help="y value to plot",
metavar=" STRING")
parser.add_option("--z",
action="store",
dest="z",
type="string",
default = None,
help="z value to plot",
metavar=" STRING")
parser.add_option("--legend",
action="store_true",
dest="legend",
default = False,
help="add a legend")
parser.add_option("--coolingdir",
action="store",
dest="coolingdir",
default = "/home/epfl/revaz/code/gear/PyCool/tables_wiersma/coolingtables/",
help="cooling directory")
parser.add_option("--HeMfaction",
action="store",
dest="HeMfaction",
type="float",
default = 0.2806,
help="Helium mass fraction",
metavar=" FLOAT")
parser.add_option("--nH",
action="store",
dest="nH",
type="float",
default = 1e-3,
help="Hydrogen density",
metavar=" FLOAT")
parser.add_option("--Z",
action="store",
dest="Z",
type="float",
default = 0.02,
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("--coolingfile",
action="store",
dest="coolingfile",
default = "cooling.dat",
help="cooling file")
parser.add_option("--Xi",
action="store",
dest="Xi",
type="float",
default = 0.76,
help="Hydrogen mass fraction",
metavar=" FLOAT")
parser.add_option("--gamma",
action="store",
dest="gamma",
type="float",
default = 5/3.,
help="Gas adiabatic index",
metavar=" FLOAT")
(options, args) = parser.parse_args()
#pt.check_files_number(args)
files = args
return files,options
########################################################################
# MAIN
########################################################################
files,opt = parse_options()
FeHSolar=0.00176604
Zsol =0.02
param = {}
param["CoolingType"]=2
param["CutofCoolingTemperature"]=1e1
param["InitGasMetallicity"]=-10
param["CoolingParameters_FeHSolar"]=FeHSolar
param["HubbleParam"]=1.0
files,opt = parse_options()
# init cooling using cooling file
cooling.init_cooling(opt.coolingfile,param)
# get parameters
#params = cooling.GetParameters()
# init cooling
cooling.init_cooling_with_metals()
#check tables
#cooling.check_cooling_table_with_metals()
HeMfaction = opt.HeMfaction
nH = opt.nH
Z = opt.Z;
Z_fe= opt.Z; # bof
Redshift=opt.Redshift
Metal = Z
Ts = 10**arange(1,9,0.02)
nHs = 10**arange(-7,0,0.02)
Ls = zeros( (len(Ts),len(nHs)) )
for i,T in enumerate(Ts):
for j,nH in enumerate(nHs):
Ls[i,j] = cooling.cooling_function_with_metals(T,Metal)
#if log10(T)<4:
# Ls[i,j] = 1e-22
########################################################################
# plot
########################################################################
#data = log10(fabs(Ls))
# rho : atom/cc
# T : K
# Lambda : normalised [erg cm3 / s]
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 = opt.Xi
gamma = opt.gamma
pars = {}
pars['k'] = k
pars['mh'] = mh
pars['gamma'] = gamma
pars['mu'] = 1 # unused here
rho = nHs * density_units.convertionFactorTo(cgs_units.UnitDensity) # to gr/cm3
U = thermodyn.Urt(rho,Ts,pars)
#Ls = where(Ls<0,1e-15,Ls)
Tcool = zeros(Ls.shape)
for i in range(Ls.shape[0]):
Tcool[i,:] = mh**2 / (fabs(Ls[i,:]*Xi**2)) * U[i]/rho
Tcool = Tcool * cgs_units.convertionFactorTo(time_units.UnitTime)
data = log10(Tcool)
# set ranges
xmin=min(log10(nHs))
xmax=max(log10(nHs))
ymin=min(log10(Ts))
ymax=max(log10(Ts))
vmin = -2
vmax = 10
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]}$")
zlabel = r'$\rm{log\,\,Cooling\,Time}\,\left[ \rm{Myr} \right]$'
cb = pt.colorbar()
cb.set_label(zlabel)
pt.show()
c4science · Help