Page MenuHomec4science
Files F100639310

gplot_rho-A
No OneTemporary
Actions

File Metadata

Created
Sat, Feb 1, 09:51

gplot_rho-A
View Options

#!/usr/bin/env python
'''
Plot entropie of the model as a function of radius
assuming an ideal gas
Yves Revaz
Thu Feb 23 15:00:11 CET 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_log_options(parser)
parser.add_option("-t",
action="store",
dest="ftype",
type="string",
default = None,
help="type of the file",
metavar=" TYPE")
parser.add_option("--cgs",
action="store_true",
dest="cgs",
default = 0,
help="invert into cgs units")
parser.add_option("--gamma",
action="store",
dest="gamma",
type="float",
default=5/3.,
help="adiabatic index")
parser.add_option("--mu",
action="store",
dest="mu",
type="float",
default=2,
help="mean molecular mass")
parser.add_option("--av",
action="store",
dest="av",
type="float",
default=None,
help="Av consante (in cgs)")
parser.add_option("--bv",
action="store",
dest="bv",
type="float",
default=None,
help="Bv consante (in cgs)")
(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
cgs = options.cgs
ftype = options.ftype
gamma = options.gamma
mu = options.mu
av = options.av
bv = options.bv
units = get_units_options(options)
#######################################
# open sm
#######################################
g = Graph_Init(ps)
Graph_SetDefaultsGraphSettings(g)
colors = Graph_SetColorsForFiles(files,col)
#######################################
# LOOP
#######################################
# read files
for file in files:
nbody = Nbody(file,ftype=ftype)
nbody = nbody.select('gas')
nbody.hdcenter()
x = nbody.rho.astype(Float)
y = nbody.u.astype(Float)
# ici, il faut convertire toute les constantes dans le user units
# en utilisant la variable units
# ctes.convert_ctes(units)
UnitLength_in_cm = units[0]
UnitMass_in_g = units[1]
UnitVelocity_in_cm_per_s = units[2]
UnitTime_in_s = UnitLength_in_cm / UnitVelocity_in_cm_per_s
UnitEnergy_in_cgs = UnitMass_in_g*UnitVelocity_in_cm_per_s**2
BOLTZMANN = BOLTZMANN/UnitEnergy_in_cgs
PROTONMASS = PROTONMASS/UnitMass_in_g
if av != None:
AV = av
if bv != None:
BV = bv
AV = AV/UnitMass_in_g/UnitLength_in_cm**5*UnitTime_in_s**2
BV = BV/UnitLength_in_cm**3
# convert into cgs
if cgs:
x = Density2cgs(x,units)
y = EnergySpec2cgs(y,units)
BOLTZMANN = BOLTZMANN * UnitEnergy_in_cgs
PROTONMASS = PROTONMASS* UnitMass_in_g
AV = AV *UnitMass_in_g *UnitLength_in_cm**5 /UnitTime_in_s**2
BV = BV *UnitLength_in_cm**3
# define params
pars = {"k":BOLTZMANN,"mh":PROTONMASS,"mu":mu,"gamma":gamma,"a":AV,"b":BV}
# convert energy in temperature
y = vw.Aru(x,y,pars)
# use log
if log != None:
x,y = Graph_UseLog(x,y,log)
if file == files[0]:
xmin,xmax,ymin,ymax = Graph_SetLimits(g,xmin,xmax,ymin,ymax,x,y)
g.box()
# plot points
g.ctype(colors[file])
g.points(x,y)
###############################
# add P=0 limit
###############################
x = arange(xmin,xmax,(xmax-xmin)/1000)
if log=='x' or log=='xy' or log=='yx':
x = 10**x
y = vw.Aneg(x,pars)
x,y = Graph_UseLog(x,y,log)
g.ctype(128)
g.ltype(0)
g.connect(x,y)
g.ctype(0)
g.ltype(0)
###############################
# add A limit
###############################
x = arange(xmin,xmax,(xmax-xmin)/1000)
if log=='x' or log=='xy' or log=='yx':
x = 10**x
y = vw.Amin(x,pars)
x,y = Graph_UseLog(x,y,log)
g.ctype(192)
g.ltype(0)
g.connect(x,y)
g.ctype(0)
g.ltype(0)
###############################
# add rho limit
###############################
rl = vw.Rl(pars)
if log=='x' or log=='xy' or log=='yx':
rl = log10(rl)
g.ctype(192)
g.ltype(0)
g.relocate(rl,ymin)
g.draw(rl,ymax)
g.ctype(0)
g.ltype(0)
###############################
# add extremas
###############################
T1 = arange(1e-5,vw.Tc(pars),0.1)
T2 = arange(1e-5,vw.Tc(pars),0.1)
x1,x2,x3,c = vw.Extremas(T1,pars)
A1 = vw.Art(x1,T1,pars)
A2 = vw.Art(x2,T2,pars)
x1,A1 = Graph_UseLog(x1,A1,log)
x2,A2 = Graph_UseLog(x2,A2,log)
g.ctype(64)
g.ltype(0)
g.connect(x1,A1)
g.connect(x2,A2)
g.ctype(0)
g.ltype(0)
g.ctype(0)
if log == 'xy' or log == 'yx':
g.xlabel('log rho')
g.ylabel('log A')
elif log == 'x':
g.xlabel('log rho')
g.ylabel('A')
elif log == 'y':
g.xlabel('rho')
g.ylabel('log A')
else:
g.xlabel('rho')
g.ylabel('A')
# -- end ---
Graph_End(g,ps)

Event Timeline

c4science · Help