Page MenuHomec4science
Files F100639977

gplot_madau
No OneTemporary
Actions

File Metadata

Created
Sat, Feb 1, 10:02

gplot_madau
View Options

#!/home/revaz/local/bin/python
'''
Plot star formation rate of a cosmological model as a function of redshift
Yves Revaz
ven mai 19 10:55:04 CEST 2006
'''
from numarray import *
from pNbody import *
from pNbody import cosmo
import SM
import string
import sys
import os
from Nbody.libutil import histogram
from optparse import OptionParser
from Gtools import *
from Gtools import io
from Gtools import cgtools
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = add_postscript_options(parser)
parser = add_color_options(parser)
parser = add_limits_options(parser)
parser = add_units_options(parser)
parser = add_log_options(parser)
parser.add_option("-t",
action="store",
dest="ftype",
type="string",
default = None,
help="type of the file",
metavar=" TYPE")
parser.add_option("--cgs",
action="store_true",
dest="cgs",
default = 0,
help="invert into cgs units")
parser.add_option("-r",
action="store",
dest="reduc",
type="int",
default = 1,
help="reduction number",
metavar=" INT")
parser.add_option("--boxsize",
action="store",
dest="boxsize",
type="float",
default = 1,
help="box size in Mpc",
metavar=" FLOAT")
parser.add_option("--HubbleParam",
action="store",
dest="HubbleParam",
type="float",
default = 0.7,
help="HubbleParam",
metavar=" FLOAT")
parser.add_option("--Omega0",
action="store",
dest="Omega0",
type="float",
default = 0.3,
help="Omega0",
metavar=" FLOAT")
parser.add_option("--OmegaLambda",
action="store",
dest="OmegaLambda",
type="float",
default = 0.7,
help="OmegaLambda",
metavar=" FLOAT")
parser.add_option("--Hubble",
action="store",
dest="Hubble",
type="float",
default = 0.1,
help="Hubble",
metavar=" FLOAT")
(options, args) = parser.parse_args()
if options.colors!=None:
exec("options.colors = array([%s])"%(options.colors))
if len(args) == 0:
print "you must specify a filename"
sys.exit(0)
files = args
return files,options
#############################
# graph
#############################
# get options
files,options = parse_options()
ps = options.ps
col = options.colors
xmin = options.xmin
xmax = options.xmax
ymin = options.ymin
ymax = options.ymax
log = options.log
cgs = options.cgs
ftype = options.ftype
reduc = options.reduc
boxsize = options.boxsize
HubbleParam = options.HubbleParam
Omega0 = options.Omega0
OmegaLambda = options.OmegaLambda
Hubble = options.Hubble
localsystem = Set_SystemUnits_From_Options(options)
#######################################
# define output system of unit
#######################################
outputunits = units.UnitSystem('mks',[units.Unit_kpc, units.Unit_Ms, units.Unit_yr , units.Unit_K, units.Unit_mol, units.Unit_C])
#######################################
# open sm
#######################################
g = Graph_Init(ps)
Graph_SetDefaultsGraphSettings(g)
colors = Graph_SetColorsForFiles(files,col)
#######################################
# LOOP
#######################################
# read files
for file in files:
try:
vals = io.read_new_energy(file)
except "NotNewEgyFileError":
vals = io.read_energy(file,iobs=None)
a = vals['Time']
m = vals['MassComp5']
# reduc
n = reduc
c = fmod(arange(len(a)),n)
a = compress((c==0),a)
m = compress((c==0),m)
# other implementation to get dt
# convert a into t
t = cgtools.Age_a(a,Omega0,OmegaLambda,Hubble)
dt = -(t[1:] - t[0:-1])
# other implementation to get dt
da = (a[1:] - a[0:-1])
dt = da/cosmo.Adot_a(0.5*(a[1:] + a[0:-1]))
sfr = (m[1:] - m[0:-1])/dt
x = cosmo.Z_a(a[1:])
y = sfr
# convert into Msol/yr / Mpc^3
# (1) : 10**10/h Msol -> Msol
y = y*10**10/HubbleParam
# (2) : Gyr -> yr
y = y/ (0.978e9/HubbleParam)
# (3) : -> Mpc^3
y = y/ (boxsize**3)
# use the log
y = log10(y)
# use log
if log != None:
x,y = Graph_UseLog(x,y,log)
g.ticksize(-1,0,-1,-1)
if file == files[0]:
xmin,xmax,ymin,ymax = Graph_SetLimits(g,xmin,xmax,ymin,ymax,x,y)
Graph_DrawBox(g,xmin,xmax,ymin,ymax,log)
# plot points
g.ctype(colors[file])
g.connect(x,y)
g.ctype(0)
if log == 'xy' or log == 'yx':
g.xlabel('log Redshift')
g.ylabel('log SFR')
elif log == 'x':
g.xlabel('log Redshift')
g.ylabel('SFR')
elif log == 'y':
g.xlabel('Redshift')
g.ylabel('log SFR')
else:
g.xlabel('Redshift')
g.ylabel('SFR')
# -- end ---
Graph_End(g,ps)

Event Timeline

c4science · Help