Homec4science
Diffusion Gtools (master) Browse

gtools/pscripts/pSfrvsTime
dfeae5220b74master

Last Change
Raw File

pSfrvsTime
View Options

#!/usr/bin/env 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 == "SFR":
opt.integrate = True
opt.derive = True
opt.interpolate=None
opt.obs = "mnewstars"
opt.rf = 2
opt.rfdx = 10.
opt.ylabel = r'$\rm{Sfr}\,\left[ M_\odot/yr \right]$'
opt.xlabel = r'$\rm{Time}\,\left[ Gyr \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)
if opt.obs == "MSTARS":
opt.integrate = False
opt.derive = False
opt.obs = "mstars"
opt.ylabel = r'$\rm{Mass}\,\left[ 10^6 \,\,M_\odot \right]$'
opt.xlabel = r'$\rm{Time}\,\left[ Gyr \right]$'
opt.ux = local_units.convertionFactorTo(out_units_TIME.UnitTime)
opt.uy = local_units.convertionFactorTo(out_units_MASS.UnitMass)
if opt.obs == "MNEWSTARS":
opt.integrate = True
opt.derive = False
opt.obs = "mnewstars"
opt.ylabel = r'$\rm{Mass}\,\left[ 10^6 \,\,M_\odot \right]$'
opt.xlabel = r'$\rm{Time}\,\left[ Gyr \right]$'
opt.ux = local_units.convertionFactorTo(out_units_TIME.UnitTime)
opt.uy = local_units.convertionFactorTo(out_units_MASS.UnitMass)
# 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,mstars,mnewstars,sfr = pt.io.read_ascii(file,[0,1,2,3,4])
data = {}
data['t'] = t
data['dt'] = dt
data['mstars'] = mstars
data['mnewstars'] = mnewstars
data['sfr'] = sfr
data['mstars'] = data['mstars']
data['mnewstars'] = data['mnewstars']
data['T'] = data['t']
data['t'] = data['t']
data['dt'] = data['dt']
data['sfr'] = sfr
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)
c4science · Help