Page MenuHomec4science

gget_values
No OneTemporary

File Metadata

Created
Mon, Nov 25, 23:10

gget_values

#!/usr/bin/env python
'''
Get the values of different observables as a function of others
Yves Revaz
mar fev 28 19:06:48 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_units_options(parser)
parser.add_option("--MeanWeight",
action="store",
dest="MeanWeight",
type="float",
default=2,
help="MeanWeight")
parser.add_option("--GammaCte",
action="store",
dest="GammaCte",
type="float",
default=5/3.,
help="GammaCte")
parser.add_option("--Rho",
action="store",
dest="Rho",
type="float",
default=None,
help="Rho : density")
parser.add_option("--U",
action="store",
dest="U",
type="float",
default=None,
help="U : spec. energy")
parser.add_option("--A",
action="store",
dest="A",
type="float",
default=None,
help="A : entropy")
parser.add_option("--P",
action="store",
dest="P",
type="float",
default=None,
help="P : pressure")
parser.add_option("--T",
action="store",
dest="T",
type="float",
default=None,
help="T : temperature")
parser.add_option("-o","--obs",
action="store",
dest="obs",
type="string",
default=None,
help="observable to compute")
(options, args) = parser.parse_args()
return options
#############################
# main
#############################
# get options
options = parse_options()
GammaCte = options.GammaCte
MeanWeight = options.MeanWeight
Rho = options.Rho
U = options.U
A = options.A
P = options.P
T = options.T
obs = options.obs
units = get_units_options(options)
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;
UnitDensity_in_cgs = UnitMass_in_g / pow(UnitLength_in_cm, 3);
UnitPressure_in_cgs = UnitMass_in_g / UnitLength_in_cm / pow(UnitTime_in_s, 2);
UnitCoolingRate_in_cgs = UnitPressure_in_cgs / UnitTime_in_s;
UnitEnergy_in_cgs = UnitMass_in_g * pow(UnitLength_in_cm, 2) / pow(UnitTime_in_s, 2);
# convert into user units
Boltzmann = BOLTZMANN /UnitEnergy_in_cgs
ProtonMass = PROTONMASS/UnitMass_in_g
MuMh = MeanWeight*ProtonMass
Av = AV/UnitMass_in_g/UnitLength_in_cm**5*UnitTime_in_s**2
Bv = BV/UnitLength_in_cm**3
pars = {'k':Boltzmann,'mh':ProtonMass,'mu':MeanWeight,'gamma':GammaCte,'a':Av,'b':Bv}
##################################################
# COMPUTE OBSERVABLE
##################################################
if obs == None:
pass
elif obs == 'Prt':
if Rho==None or T==None:
print "you must specify Rho and T"
sys.exit()
else:
print vw.Prt(Rho,T,pars)
elif obs == 'Trp':
if Rho==None or P==None:
print "you must specify Rho and P"
sys.exit()
else:
print vw.Trp(Rho,P,pars)
elif obs == 'Art':
if Rho==None or T==None:
print "you must specify Rho and T"
sys.exit()
else:
print vw.Art(Rho,T,pars)
elif obs == 'Tra':
if Rho==None or A==None:
print "you must specify Rho and A"
sys.exit()
else:
print vw.Tra(Rho,A,pars)
elif obs == 'Urt':
if Rho==None or T==None:
print "you must specify Rho and T"
sys.exit()
else:
print vw.Urt(Rho,T,pars)
elif obs == 'Tru':
if Rho==None or U==None:
print "you must specify Rho and U"
sys.exit()
else:
print vw.Tru(Rho,U,pars)
elif obs == 'Pra':
if Rho==None or A==None:
print "you must specify Rho and A"
sys.exit()
else:
print vw.Pra(Rho,A,pars)
elif obs == 'Arp':
if Rho==None or P==None:
print "you must specify Rho and P"
sys.exit()
else:
print vw.Arp(Rho,P,pars)
elif obs == 'Pru':
if Rho==None or U==None:
print "you must specify Rho and U"
sys.exit()
else:
print vw.Pru(Rho,U,pars)
elif obs == 'Urp':
if Rho==None or P==None:
print "you must specify Rho and P"
sys.exit()
else:
print vw.Urp(Rho,P,pars)
elif obs == 'Ura':
if Rho==None or A==None:
print "you must specify Rho and A"
sys.exit()
else:
print vw.Ura(Rho,A,pars)
elif obs == 'Aru':
if Rho==None or U==None:
print "you must specify Rho and U"
sys.exit()
else:
print vw.Aru(Rho,U,pars)

Event Timeline