Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F101347921
pAccretionvsTime
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
Sat, Feb 8, 07:06
Size
7 KB
Mime Type
text/x-python
Expires
Mon, Feb 10, 07:06 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
24140730
Attached To
rGTOOLS Gtools
pAccretionvsTime
View Options
#!/usr/bin/python
'''
Extract and plot sfr contained in the
output Gadget file called by default "sfr.txt".
Yves Revaz
Mon Feb 2 19:15:39 CET 2009
'''
import Ptools as pt
from numpy import *
from pNbody import *
import string
import sys
import os
from pNbody.libutil import histogram
from optparse import OptionParser
from Gtools import *
from Gtools import io
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_postscript_options(parser)
parser = pt.add_limits_options(parser)
parser = pt.add_log_options(parser)
parser = pt.add_cmd_options(parser)
parser = pt.add_legend_options(parser)
parse = pt.add_units_options(parser)
parser.add_option("-o",
action="store",
dest="obs",
type="string",
default = 'EnergyTot',
help="observable name",
metavar=" NAME")
parser.add_option("--relative",
action="store_true",
dest="rel",
default = 0,
help="plot relative value")
parser.add_option("--nc",
action="store",
dest="nc",
type="float",
default = 1000,
help="number of points per bin")
parser.add_option("--rf",
action="store",
dest="rf",
type="float",
default = 1,
help="reduction factor")
parser.add_option("--rfn",
action="store",
dest="rfn",
type="int",
default = None,
help="reduction factor (number of points)")
parser.add_option("--rfmn",
action="store",
dest="rfmn",
type="float",
default = None,
help="reduction factor (min)")
parser.add_option("--rfmx",
action="store",
dest="rfmx",
type="float",
default = None,
help="reduction factor (max)")
parser.add_option("--rfdx",
action="store",
dest="rfdx",
type="float",
default = None,
help="reduction factor (dx)")
parser.add_option("--integrate",
action="store_true",
dest="integrate",
default = 0,
help="integrate values")
parser.add_option("--interpolate",
action="store_true",
dest="interpolate",
default = 0,
help="interpolate values")
parser.add_option("--derive",
action="store_true",
dest="derive",
default = 0,
help="derive values")
(options, args) = parser.parse_args()
if len(args) == 0:
print "you must specify a filename"
sys.exit(0)
files = args
return files,options
#######################################
# MakePlot
#######################################
def MakePlot(files,opt):
# some inits
colors = pt.Colors(n=len(files))
linestyles = pt.LineStyles()
# define local units
unit_params = pt.do_units_options(opt)
local_units = units.Set_SystemUnits_From_Params(unit_params)
out_units_TIME = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Gyr,units.Unit_K])
out_units_SFR = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_yr,units.Unit_K])
out_units_MASS = units.UnitSystem('local',[units.Unit_kpc,1e6*units.Unit_Ms,units.Unit_yr,units.Unit_K])
opt.ux = 1 # output time units in code units
opt.uy = 1 # output mass units in code units
labels = []
opt.xlabel = r'$\textrm{Time}$'
opt.ylabel = r'$\textrm{%s}$'%(opt.obs)
if opt.obs == "MT":
opt.integrate = True
opt.derive = False
opt.obs = "dM"
opt.xlabel = r'$\rm{Time}\,\left[ Gyr \right]$'
opt.ylabel = r'$\rm{Mass}\,\left[ 10^6 \,\,M_\odot \right]$'
opt.ux = local_units.convertionFactorTo(out_units_TIME.UnitTime)
opt.uy = local_units.convertionFactorTo(out_units_MASS.UnitMass)
if opt.obs == "MT_th":
opt.integrate = True
opt.derive = False
opt.obs = "dM_th"
opt.xlabel = r'$\rm{Time}\,\left[ Gyr \right]$'
opt.ylabel = r'$\rm{Mass}\,\left[ 10^6 \,\,M_\odot \right]$'
opt.ux = local_units.convertionFactorTo(out_units_TIME.UnitTime)
opt.uy = local_units.convertionFactorTo(out_units_MASS.UnitMass)
if opt.obs == "DMDT":
opt.integrate = True
opt.derive = True
opt.interpolate=None
opt.obs = "dM"
opt.xlabel = r'$\rm{Time}\,\left[ Gyr \right]$'
opt.ylabel = r'$\rm{Sfr}\,\left[ M_\odot/yr \right]$'
opt.ux = local_units.convertionFactorTo(out_units_TIME.UnitTime)
opt.uy = local_units.convertionFactorTo(out_units_SFR.UnitMass)/local_units.convertionFactorTo(out_units_SFR.UnitTime)
opt.rf = 100
opt.rfdx = 10.
if opt.obs == "DMDT_th":
opt.integrate = True
opt.derive = True
opt.interpolate=None
opt.obs = "dM_th"
opt.xlabel = r'$\rm{Time}\,\left[ Gyr \right]$'
opt.ylabel = r'$\rm{Sfr}\,\left[ M_\odot/yr \right]$'
opt.ux = local_units.convertionFactorTo(out_units_TIME.UnitTime)
opt.uy = local_units.convertionFactorTo(out_units_SFR.UnitMass)/local_units.convertionFactorTo(out_units_SFR.UnitTime)
opt.rf = 100
opt.rfdx = 10.
# observable
opt.obs = string.split(opt.obs,',')
# observable
opt = pt.do_legend_options(opt)
datas=[]
# read files
for ifile,file in enumerate(files):
t,dt,NaccRemain,MRemain,dM,p,Macctot = pt.io.read_ascii(file,[0,1,2,3,4,5,6])
data = {}
data['t'] = t
data['dt'] = dt
data['AccreatedMass'] = MRemain[0]-MRemain
data['dM_th'] = dM
data['dMdt_th'] = dt
data['dM'] = Macctot
data['T'] = data['t']
Ti = arange(0,data['T'][-1],10)
for obs in opt.obs:
if opt.legend_txt!=None:
label = opt.legend_txt[ifile]
else:
label = file
d = pt.DataPoints(data['t'],data[obs],label=label)
d.color = colors.get()
datas.append(d)
# reduction
if opt.rf>1:
#for d in datas:
d.reduc(opt.rf,opt.rfn,opt.rfmn,opt.rfmx,opt.rfdx)
# integrate
if opt.integrate:
#for d in datas:
d.integrate()
# interpolate
if opt.interpolate:
#for d in datas:
d.interpolate(Ti)
d.x = d.xi
d.y = d.yi
# derive
if opt.derive:
#for d in datas:
d.derive()
# set final units
#for d in datas:
d.x = d.x * opt.ux
d.y = d.y * opt.uy
# plot points
##for d in datas:
pt.plot(d.x,d.y,color=d.color,label=d.label)
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
# labels
pt.xlabel(opt.xlabel,fontsize=pt.labelfont)
pt.ylabel(opt.ylabel,fontsize=pt.labelfont)
if opt.legend:
pt.LegendFromDataPoints(datas,opt.legend_loc)
if __name__ == '__main__':
files,opt = parse_options()
pt.InitPlot(files,opt)
#pt.figure(figsize=(8*2,6*2))
#pt.figure(dpi=10)
pt.pcolors
#fig = pt.gcf()
#fig.subplots_adjust(left=0.1)
#fig.subplots_adjust(right=1)
#fig.subplots_adjust(bottom=0.12)
#fig.subplots_adjust(top=0.95)
#fig.subplots_adjust(wspace=0.25)
#fig.subplots_adjust(hspace=0.02)
MakePlot(files,opt)
pt.EndPlot(files,opt)
Event Timeline
Log In to Comment