Page MenuHomec4science
Files F106603541

pycool_LHvsT
No OneTemporary
Actions

File Metadata

Created
Fri, Mar 28, 04:18

pycool_LHvsT
View Options

#!/usr/bin/env python
from optparse import OptionParser
import Ptools as pt
from pNbody import *
from pNbody import units
import string
from scipy import optimize
from PyCool import cooling
UnitLength_in_cm = 3.085e+21
UnitMass_in_g = 1.989e+43
UnitVelocity_in_cm_per_s = 20725573.785998672
UnitTime_in_s = 148849920000000.0
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_postscript_options(parser)
parser = pt.add_ftype_options(parser)
parser = pt.add_reduc_options(parser)
parser = pt.add_center_options(parser)
parser = pt.add_select_options(parser)
parser = pt.add_cmd_options(parser)
parser = pt.add_display_options(parser)
parser = pt.add_info_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")
(options, args) = parser.parse_args()
#pt.check_files_number(args)
files = args
return files,options
########################################################################
# MAIN
########################################################################
files,opt = parse_options()
Redshift=10
X=0.76
n_H=2.8e-3
# do
cooling.init_cooling()
# get parameters
params = cooling.GetParameters()
# set the redshift and compute flux, a.s.o
cooling.init_from_new_redshift(Redshift)
# compute elements densities
logTs = arange(1,8,0.01)
n = len(logTs)
n_Hs = zeros(n)
n_HIs = zeros(n)
n_HIIs = zeros(n)
n_HEIs = zeros(n)
n_HEIIs = zeros(n)
n_HEIIIs = zeros(n)
n_Es= zeros(n)
mus = zeros(n)
c1s = zeros(n)
c2s = zeros(n)
c3s = zeros(n)
c4s = zeros(n)
c5s = zeros(n)
c6s = zeros(n)
c7s = zeros(n)
c8s = zeros(n)
c9s = zeros(n)
c10s = zeros(n)
c11s = zeros(n)
c12s = zeros(n)
c13s = zeros(n)
h1s = zeros(n)
h2s = zeros(n)
h3s = zeros(n)
h4s = zeros(n)
for i,logT in enumerate(logTs):
T=10**logT
n_H,n_HI,n_HII,n_HEI,n_HEII,n_HEIII,n_E,mu = cooling.compute_densities(T,X,n_H)
c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,h1,h2,h3,h4=cooling.compute_cooling(T,X,n_H,n_HI,n_HII,n_HEI,n_HEII,n_HEIII,n_E,mu)
n_Hs[i]= n_H
n_HIs[i]= n_HI
n_HIIs[i]= n_HII
n_HEIs[i]= n_HEI
n_HEIIs[i]= n_HEII
n_HEIIIs[i]= n_HEIII
n_Es[i]= n_E
mus[i]= mu
c1s [i]= c1
c2s[i]= c2
c3s[i]= c3
c4s[i]= c4
c5s[i]= c5
c6s[i]= c6
c7s[i]= c7
c8s[i]= c8
c9s[i]= c9
c10s[i]= c10
c11s[i]= c11
c12s[i]= c12
c13s[i]= c13
h1s[i]= h1
h2s[i]= h2
h3s[i]= h3
h4s[i]= h4
############################
# plot
############################
import Ptools as pt
vz=-30
chs = c1s+c2s+c3s+c4s+c5s+c6s+c7s+c8s+c9s+c10s+c11s+c12s+c13s - (h1s+h2s+h3s+h4s)
chs = fabs(chs)
logc1s = where(c1s >0,log10(c1s),vz)
logc2s = where(c2s >0,log10(c2s),vz)
logc3s = where(c3s >0,log10(c3s),vz)
logc4s = where(c4s >0,log10(c4s),vz)
logc5s = where(c5s >0,log10(c5s),vz)
logc6s = where(c6s >0,log10(c6s),vz)
logc7s = where(c7s >0,log10(c7s),vz)
logc8s = where(c8s >0,log10(c8s),vz)
logc9s = where(c9s >0,log10(c9s),vz)
logc10s = where(c10s>0,log10(c10s),vz)
logc11s = where(c11s>0,log10(c11s),vz)
logc12s = where(c12s>0,log10(c12s),vz)
logc13s = where(c13s>0,log10(c13s),vz)
logh1s = where(h1s >0,log10(h1s),vz)
logh2s = where(h2s >0,log10(h2s),vz)
logh3s = where(h3s >0,log10(h3s),vz)
logh4s = where(h4s >0,log10(h4s),vz)
logchs = where(chs >0,log10(chs),vz)
pt.plot(logTs,logc1s)
pt.plot(logTs,logc2s)
pt.plot(logTs,logc3s)
pt.plot(logTs,logc4s)
pt.plot(logTs,logc5s)
pt.plot(logTs,logc6s)
pt.plot(logTs,logc7s)
pt.plot(logTs,logc8s)
pt.plot(logTs,logc9s)
pt.plot(logTs,logc10s)
pt.plot(logTs,logc11s)
pt.plot(logTs,logc12s)
pt.plot(logTs,logc13s)
pt.plot(logTs,logh1s,'--')
pt.plot(logTs,logh2s,'--')
pt.plot(logTs,logh3s,'--')
pt.plot(logTs,logh4s,'--')
pt.plot(logTs,logchs,'k-')
pt.axis([3.5,8,-25,-21])
pt.axis([0,8,-30,-20])
pt.show()

Event Timeline

c4science · Help