Homec4science
Diffusion pNbody (master) Browse

pNbody/examples/plots/star_formation_rate/plot_star_formation_rate.py
179d38abe82dmaster

Last Change
Raw File

plot_star_formation_rate.py
View Options

#!/usr/bin/env python
'''
Plot value computed in a spherical grid, as a function of the radius
Yves Revaz
Wed Jan 9 09:29:46 CET 2008
'''
import Ptools as pt
from pNbody import *
from pNbody import myNumeric
from pNbody import libdisk
from pNbody import libgrid
from pNbody import profiles
from pNbody import cosmo
import string
import sys
import os
from numpy import *
from optparse import OptionParser
from scipy.optimize import leastsq
from scipy.optimize import bisection
#from mpfit import *
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)
parser = pt.add_legend_options(parser)
parser.add_option("--Rmax",
action="store",
dest="Rmax",
type="float",
default = 50.,
help="max radius of bins",
metavar=" FLOAT")
parser.add_option("--nR",
action="store",
dest="nR",
type="int",
default = 32,
help="number of bins in r",
metavar=" INT")
parser.add_option("--nt",
action="store",
dest="nt",
type="int",
default = 32,
help="number of bins in r",
metavar=" INT")
parser.add_option("--nz",
action="store",
dest="nz",
type="int",
default = 65,
help="number of bins in z",
metavar=" INT")
parser.add_option("--ErrTolTheta",
action="store",
dest="ErrTolTheta",
type="float",
default = 0.5,
help="Error tolerance theta",
metavar=" FLOAT")
parser.add_option("--fct",
action="store",
dest="fct",
type="string",
default = None,
help="transformation function ex : 'r**3' should be of the form M(r)",
metavar=" STR")
parser.add_option("--fctm",
action="store",
dest="fctm",
type="string",
default = None,
help="inverse transformation function ex : 'r**(1/3.)'",
metavar=" STR")
parser.add_option("--AdaptativeSoftenning",
action="store_true",
dest="AdaptativeSoftenning",
default = False,
help="AdaptativeSoftenning")
parser.add_option("--eps",
action="store",
dest="eps",
type="float",
default = 0.28,
help="smoothing length",
metavar=" FLOAT")
parser.add_option("--dt",
action="store",
dest="dt",
type="float",
default = 10,
help="max star age in Myr = inerval time",
metavar=" FLOAT")
(options, args) = parser.parse_args()
if len(args) == 0:
print "you must specify a filename"
sys.exit(0)
files = args
return files,options
#######################################
# units
#######################################
# gear system of units
params = {}
params['UnitVelocity_in_cm_per_s'] = 20725573.785998672
params['UnitMass_in_g'] = 1.989e+43
params['UnitLength_in_cm'] = 3.085e+21
gear_units = units.Set_SystemUnits_From_Params(params)
out_units = units.UnitSystem('local',[units.Unit_pc,units.Unit_Msol,units.Unit_yr,units.Unit_K])
tokms = gear_units.convertionFactorTo(out_units.UnitVelocity)
toMsol = gear_units.convertionFactorTo(out_units.UnitMass)
toMsolpc2 = gear_units.convertionFactorTo(out_units.UnitSurfaceDensity)
topc2 = gear_units.convertionFactorTo(out_units.UnitSurface)
toyr = gear_units.convertionFactorTo(out_units.UnitTime)
toMyr= toyr/1e6
#######################################
# MakePlot
#######################################
def MakePlot(files,opt):
# some inits
colors = pt.Colors(n=len(files))
datas = []
# read files
for i,file in enumerate(files):
nb = Nbody(file,ftype=opt.ftype)
# apply options
nb = pt.do_reduc_options(nb,opt)
nb = pt.do_select_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)
#######################
# set the grid scaling
#######################
if opt.fct!=None:
print "g =lambda r:(%s)"%opt.fct
exec("g =lambda r:(%s)"%opt.fct)
else:
g = None
if opt.fctm!=None:
print "gm =lambda r:(%s)"%opt.fctm
exec("gm =lambda r:(%s)"%opt.fctm)
else:
gm = None
##############################
# using a cylindrical rt grid
##############################
G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.Rmax,nr=opt.nR,nt=2,z=0,g=g,gm=gm)
# radius
R,t = G.get_rt()
# star formation rate map
val = nb.Sfr(opt.dt/toMyr) * toMsol/toyr # star formation rate, Msol/yr per particles
Sfr = G.get_ValMap(nb,val/nb.mass)
Sfr = sum(Sfr,axis=1)
# surface
S = G.get_SurfaceMap(None)
S = sum(S,axis=1)*topc2
Sfr=Sfr/S
R,Sfr=pt.CleanVectorsForLogX(R,Sfr)
R,Sfr=pt.CleanVectorsForLogY(R,Sfr)
datas.append(pt.DataPoints(R, Sfr, color=colors.next(),label=i))
# plot
for d in datas:
pt.plot(d.x,d.y,color=d.color)
# set limits and draw axis
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
xlabel = r'$\rm{Radius}\,[\rm{kpc}]$'
ylabel = r'$\rm{Star\,Formation\,\,Rate}\,[\rm{M_{\odot}/yr/pc^2}]$'
# plot axis
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(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)
c4science · Help