Homec4science
Diffusion Gear (master) Browse

Gear/PyCool/scripts/pycool_NvsT
d89b118970f3master

Last Change
Raw File

pycool_NvsT
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=2
T=2e7
X=0.76
n_H=1.0
# 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,7,0.1)
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)
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)
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
############################
# plot
############################
import Ptools as pt
pt.plot(logTs,n_HIs)
pt.plot(logTs,n_HIIs)
pt.plot(logTs,n_HIIs+n_HIs)
pt.plot(logTs,n_HEIs)
pt.plot(logTs,n_HEIIs)
pt.plot(logTs,n_HEIIIs)
pt.plot(logTs,n_HEIs+n_HEIIs+n_HEIIIs)
pt.plot(logTs,n_Es,'.')
pt.plot(logTs,mus,'--')
pt.show()
c4science · Help