Page MenuHomec4science
Files F103939836

gplot_r-masses
No OneTemporary
Actions

File Metadata

Created
Wed, Mar 5, 13:56

gplot_r-masses
View Options

#!/home/revaz/local/bin/python
'''
Plot mass/mass fraction/total mass of each components, as a function of radius
Yves Revaz
mar sep 26 11:06:07 CEST 2006
'''
from numarray import *
from Nbody import *
import SM
import string
import sys
import os
from libjeans import *
from Nbody.libutil import histogram
from optparse import OptionParser
from Gtools import *
from Gtools import vanderwaals as vw
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = add_postscript_options(parser)
parser = add_color_options(parser)
parser = add_limits_options(parser)
parser = add_units_options(parser)
parser = add_reduc_options(parser)
parser = add_center_options(parser)
parser.add_option("-t",
action="store",
dest="ftype",
type="string",
default = None,
help="type of the file",
metavar=" TYPE")
parser.add_option("--nbins",
action="store",
dest="nbins",
type="int",
default=25,
help="number of bins")
parser.add_option("--cgs",
action="store_true",
dest="cgs",
default = 0,
help="invert into cgs units")
parser.add_option("--fmass",
action="store_true",
dest="fmass",
default = 0,
help="plot mass fraction instead of masses")
parser.add_option("--dmmass",
action="store",
dest="dmmass",
type="float",
default = 0.,
help="dark matter mass in r=xmax")
(options, args) = parser.parse_args()
if options.colors!=None:
exec("options.colors = array([%s])"%(options.colors))
if len(args) == 0:
print "you must specify a filename"
sys.exit(0)
files = args
return files,options
#############################
# graph
#############################
# get options
files,options = parse_options()
ps = options.ps
col = options.colors
xmin = options.xmin
xmax = options.xmax
ymin = options.ymin
ymax = options.ymax
cgs = options.cgs
ftype = options.ftype
reduc = options.reduc
center = options.center
nbins = options.nbins
fmass = options.fmass
dmmass= options.dmmass
localsystem = Set_SystemUnits_From_Options(options)
#######################################
# define output system of unit
#######################################
outputunits = UnitSystem('mks',[Unit_kpc, Unit_Ms, Unit_yr , Unit_K, Unit_mol, Unit_C])
#######################################
# open sm
#######################################
g = Graph_Init(ps)
Graph_SetDefaultsGraphSettings(g)
colors = Graph_SetColorsForFiles(files,col)
#######################################
# LOOP
#######################################
# read files
for file in files:
nb = Nbody(file,ftype=ftype)
# center the model
if center=='hdcenter':
print "centering %s"%(center)
nb.hdcenter()
elif center=='histocenter':
print "centering %s"%(center)
nb.histocenter()
if reduc!= None:
print "reducing %s"%(reduc)
nb = nb.reduc(reduc)
#########
# sph
#########
# mass in r
nb1 = nb.select('sph')
xt,yt = nb1.mr(nb=nbins,rm=xmax)
dx = xt[2]-xt[1]
x1 = concatenate((xt[1:],xt[-1]+dx))
y1 = yt - concatenate((0,yt[:-1]))
x1 = x1[:-1]
y1 = y1[:-1]
#########
# sticky
#########
# mass in r
nb2 = nb.select('sticky')
xt,yt = nb2.mr(nb=nbins,rm=xmax)
dx = xt[2]-xt[1]
x2 = concatenate((xt[1:],xt[-1]+dx))
y2 = yt - concatenate((0,yt[:-1]))
x2 = x2[:-1]
y2 = y2[:-1]
#########
# stars
#########
# mass in r
nb3 = nb.select('stars')
xt,yt = nb3.mr(nb=nbins,rm=xmax)
dx = xt[2]-xt[1]
x3 = concatenate((xt[1:],xt[-1]+dx))
y3 = yt - concatenate((0,yt[:-1]))
x3 = x3[:-1]
y3 = y3[:-1]
if fmass:
# fraction de masse
yt = y1+y2+y3
y1 = y1/yt
y2 = y2/yt
y3 = y3/yt
else:
# output units
FACT = PhysCte(1,localsystem.UnitDic['kg'])
FACT = FACT.into(outputunits)
y1 = y1*FACT
y2 = y2*FACT
y3 = y3*FACT
tot = sum(y1)+sum(y2)+sum(y3)+dmmass
print
print "TotCompMass 1 (sph) = %7.1e = %5.1f%%"%(sum(y1),100*sum(y1)/tot)
print "TotCompMass 2 (sticky) = %7.1e = %5.1f%%"%(sum(y2),100*sum(y2)/tot)
print "TotCompMass 3 (stars) = %7.1e = %5.1f%%"%(sum(y3),100*sum(y3)/tot)
print "TotCompMass 4 (dm ) = %7.1e = %5.1f%%"%(dmmass,100*dmmass/tot)
print "TotMass = %7.1e = %5.1f%%"%(tot,100*tot/tot)
x = concatenate((x1,x2,x3))
y = concatenate((y1,y2,y3))
if file == files[0]:
xmin,xmax,ymin,ymax = Graph_SetLimits(g,xmin,xmax,ymin,ymax,x,y)
Graph_DrawBox(g,xmin,xmax,ymin,ymax,log)
g.ctype(colors[file])
# plot points
g.ltype(2)
g.histogram(x1,y1)
g.ltype(1)
g.histogram(x2,y2)
g.ltype(0)
g.histogram(x3,y3)
g.ltype(0)
g.ctype(0)
if log == 'xy' or log == 'yx':
g.xlabel('log R')
g.ylabel('log Mass')
elif log == 'x':
g.xlabel('log R')
g.ylabel('Mass')
elif log == 'y':
g.xlabel('R')
g.ylabel('log Mass')
else:
g.xlabel('R')
g.ylabel('Mass')
# -- end ---
Graph_End(g,ps)

Event Timeline

c4science · Help