Page MenuHomec4science

gphase
No OneTemporary

File Metadata

Created
Thu, Apr 24, 17:41
#!/usr/bin/env python
'''
Extract and plot massfraction of different phases contained in the
output Gadget file called by default "phase.txt".
Yves Revaz
jeu jun 15 14:44:34 CEST 2006
'''
from numarray import *
from Nbody import *
import SM
import string
import sys
import os
from Nbody.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 = add_postscript_options(parser)
parser = add_color_options(parser)
parser = add_limits_options(parser)
parser = add_log_options(parser)
(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
log = options.log
#######################################
# open sm
#######################################
g = Graph_Init(ps)
Graph_SetDefaultsGraphSettings(g)
colors = Graph_SetColorsForFiles(files,col)
#######################################
# LOOP
#######################################
# read files
for file in files:
Step,Time,GasPart,SphPart,StickyPart = io.read_phase(file)
x = Time
y3 = StickyPart/GasPart[0]
y2 = SphPart/GasPart[0]
y1 = (GasPart[0] - (StickyPart+SphPart))/GasPart[0]
y2 = y1 + y2
y3 = y2 + y3
y3[0] = 0
if file == files[0]:
xmin,xmax,ymin,ymax = Graph_SetLimits(g,xmin,xmax,ymin,ymax,x,y3)
g.box()
ymin = 0
ymax = 1
# coloriage
x = concatenate(([0],x ,[x[-1]]))
y1 = concatenate(([0],y1,[0]))
y2 = concatenate(([0],y2,[0]))
y3 = concatenate(([0],y3,[0]))
g.ctype(65)
g.shade(0,x,y3)
g.ctype(192)
g.shade(0,x,y2)
g.ctype(165)
g.shade(0,x,y1)
g.ctype(0)
# plot points
g.ctype(colors[file])
g.connect(x[1:-1],y1[1:-1])
g.connect(x[1:-1],y2[1:-1])
# labels
obs = 'mass fraction'
g.ctype(0)
g.xlabel('T')
g.ylabel('%s'%obs)
g.ctype(0)
if log == 'xy' or log == 'yx':
g.xlabel('log T')
g.ylabel('log %s'%obs)
elif log == 'x':
g.xlabel('log T')
g.ylabel('%s'%obs)
elif log == 'y':
g.xlabel('T')
g.ylabel('log %s'%obs)
else:
g.xlabel('T')
g.ylabel('%s'%obs)
# -- end ---
Graph_End(g,ps)

Event Timeline