Page MenuHomec4science
Files F101685039

gplot_T-P
No OneTemporary
Actions

File Metadata

Created
Wed, Feb 12, 18:19

gplot_T-P
View Options

#!/usr/bin/env python
'''
Plot temperature 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()
z = 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:
z = Density2cgs(z,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}
# compute temperature from energy
x = vw.Tru(z,y,pars)
# compute pression from energy
y = vw.Pru(z,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)
###############################
# critical point
###############################
Tc = log10(vw.Tc(pars))
Pc = log10(vw.Pc(pars))
g.ctype(64)
g.relocate(Tc,Pc)
g.putlabel(5,'X')
g.ctype(0)
###############################
# add extremas
###############################
T = arange(1e-5,vw.Tc(pars),0.1)
x1,x2,x3,c = vw.Extremas(T,pars)
x1 = vw.Prt(x1,T,pars)
x2 = vw.Prt(x2,T,pars)
g.ctype(64)
g.ltype(0)
T,x1 = Graph_UseLog(T,x1,log)
g.connect(T,x1)
T = arange(1e-5,vw.Tc(pars),0.1)
T,x2 = Graph_UseLog(T,x2,log)
g.connect(T,x2)
g.ctype(0)
g.ltype(0)
'''
###############################
# add T limit
###############################
zmin = -10
zmax = vw.Rc(pars)
z = arange(zmin,zmax,(zmax-zmin)/1000)
if log=='x' or log=='xy' or log=='yx':
z = 10**z
x = vw.Tmin(z,pars)
y = vw.Prt(z,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)
'''
g.ctype(0)
if log == 'xy' or log == 'yx':
g.xlabel('log T')
g.ylabel('log P')
elif log == 'x':
g.xlabel('log T')
g.ylabel('P')
elif log == 'y':
g.xlabel('T')
g.ylabel('log P')
else:
g.xlabel('T')
g.ylabel('P')
# -- end ---
Graph_End(g,ps)

Event Timeline

c4science · Help