Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65537736
pSFRvsTime
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Jun 4, 12:31
Size
7 KB
Mime Type
text/x-python
Expires
Thu, Jun 6, 12:31 (2 d)
Engine
blob
Format
Raw Data
Handle
18089986
Attached To
rGTOOLS Gtools
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
Log In to Comment