Page MenuHomec4science
Files F100634943

pSFRvsTime
No OneTemporary
Actions

File Metadata

Created
Sat, Feb 1, 08:49

pSFRvsTime
View Options

#!/usr/bin/python
from optparse import OptionParser
import Ptools as pt
from pNbody import *
from pNbody import units
from pNbody import ctes
from pNbody import thermodyn
import string
from scipy import optimize
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)
parse = pt.add_units_options(parser)
parser.add_option("-o",
action="store",
dest="obs",
type="string",
default = 'SFR',
help="observable name",
metavar=" NAME")
parser.add_option("--x",
action="store",
dest="x",
type="string",
default = 'r',
help="x value to plot",
metavar=" STRING")
parser.add_option("--nx",
action="store",
type="int",
dest="nx",
default = 500,
help="number of bins")
parser.add_option("--legend",
action="store_true",
dest="legend",
default = False,
help="add a legend")
parser.add_option("--forceComovingIntegrationOn",
action="store_true",
dest="forceComovingIntegrationOn",
default = False,
help="force the model to be in in comoving integration")
parser.add_option("--curve",
action="store_true",
dest="curve",
default = False,
help="draw curve instead of histogram")
(options, args) = parser.parse_args()
pt.check_files_number(args)
files = args
return files,options
def get_value_and_label(nb,val):
if val == 'Age':
label = r'$\rm{Age}\,\left[ \rm{Gyr} \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Gyr,units.Unit_K])
val = nb.StellarAge(units=out_units.UnitTime)
return val,label
elif val == 'CosmicTime':
label = r'$\rm{Age}\,\left[ \rm{Gyr} \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Gyr,units.Unit_K])
val = nb.CosmicTime(units=out_units.UnitTime)
return val,label
elif val == 'Redshift':
label = r'$\rm{Redshift}$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Gyr,units.Unit_K])
val = nb.Redshift()
return val,label
elif val == 'ScalingFactor':
label = r'$\rm{a}$'
val = nb.tstar
return val,label
else:
label = r'%s'%val
print "val = %s"%val
exec("val = %s"%val)
return val, label
#######################################
# MakePlot
#######################################
def MakePlot(files,opt):
# some inits
colors = pt.Colors(n=len(files))
datas = []
# read files
for file in files:
nb = Nbody(file,ftype=opt.ftype)
################
# units
################
# define local units
unit_params = pt.do_units_options(opt)
nb.set_local_system_of_units(params=unit_params)
# define output units
# nb.ToPhysicalUnits()
if opt.forceComovingIntegrationOn:
nb.setComovingIntegrationOn()
################
# apply options
################
nb = nb.select(1)
nb = pt.do_reduc_options(nb,opt)
nb = pt.do_center_options(nb,opt)
nb = pt.do_cmd_options(nb,opt)
nb = pt.do_info_options(nb,opt)
nb = pt.do_display_options(nb,opt)
################
# some info
################
print "---------------------------------------------------------"
nb.localsystem_of_units.info()
nb.ComovingIntegrationInfo()
print "---------------------------------------------------------"
if opt.obs=='SFR':
opt.x = "CosmicTime"
opt.ylabel = r'$\rm{Sfr}\,\left[ M_\odot/yr \right]$'
opt.xlabel = r'$\rm{Time}\,\left[ Gyr \right]$'
if opt.obs=='MSTARS':
opt.x = "CosmicTime"
opt.ylabel = r'$\rm{Stellar\,\,Mass}\,\left[ M_\odot \right]$'
opt.xlabel = r'$\rm{Time}\,\left[ Gyr \right]$'
if opt.obs=='MSTARSz':
opt.x = "Redshift"
opt.ylabel = r'$\rm{Stellar\,\,Mass}\,\left[ M_\odot \right]$'
opt.xlabel = r'$\rm{Redshift}$'
if opt.obs=='MSTARSa':
opt.x = "ScalingFactor"
opt.ylabel = r'$\rm{Stellar\,\,Mass}\,\left[ M_\odot \right]$'
opt.xlabel = r'$\rm{a}$'
################
# get values
################
x,xlabel = get_value_and_label(nb,opt.x)
# do the 1d histogram
G = libgrid.Generic_1d_Grid(opt.xmin,opt.xmax,opt.nx)
y = G.get_MassMap(x,nb.mass) # mass weighted histogram
x = G.get_r()
if opt.obs=='SFR':
y_out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_yr,units.Unit_K])
y = y * nb.localsystem_of_units.convertionFactorTo(units.Unit_Ms)
dt = (x[1:]-x[:-1])*1e9 # Gyrs to yrs
dMdt = y[:-1] / dt
y = dMdt
x = x[1:]
if opt.obs=='MSTARS':
y_out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_yr,units.Unit_K])
y = y * nb.localsystem_of_units.convertionFactorTo(units.Unit_Ms)
y = add.accumulate(y)
if opt.obs=='MSTARSz':
y_out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_yr,units.Unit_K])
y = y * nb.localsystem_of_units.convertionFactorTo(units.Unit_Ms)
y = y[::-1]
y = add.accumulate(y)
y = y[::-1]
if opt.obs=='MSTARSa':
y_out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_yr,units.Unit_K])
y = y * nb.localsystem_of_units.convertionFactorTo(units.Unit_Ms)
y = add.accumulate(y)
data = pt.DataPoints(x,y,color=colors.get(),label=file)
datas.append(data)
#if file == files[0]:
# xmin,xmax,ymin,ymax = pt.SetLimits(opt.xmin,opt.xmax,opt.ymin,opt.ymax,x,y,opt.log)
for d in datas:
pt.plot(d.x,d.y,color=d.color)
# set limits
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
# plot axis
pt.SetAxis(xmin,xmax,ymin,ymax,opt.log)
pt.xlabel(opt.xlabel,fontsize=pt.labelfont)
pt.ylabel(opt.ylabel,fontsize=pt.labelfont)
pt.grid(False)
if opt.legend:
pt.LegendFromDataPoints(datas)
########################################################################
# MAIN
########################################################################
if __name__ == '__main__':
files,opt = parse_options()
pt.InitPlot(files,opt)
pt.pcolors
MakePlot(files,opt)
pt.EndPlot(files,opt)

Event Timeline

c4science · Help