Page MenuHomec4science
Files F86377159

pget_VirialRadius
No OneTemporary
Actions

File Metadata

Created
Sun, Oct 6, 03:39

pget_VirialRadius
View Options

#!/usr/bin/python
from optparse import OptionParser
import Ptools as pt
from pNbody import *
from pNbody import units
from pNbody import ctes
import string
from pNbody import units,io,ctes
from scipy.optimize import bisect as bisection
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser.add_option("--dX",
action="store",
dest="dX",
type="float",
default = 200,
help="over density",
metavar=" STRING")
parse = pt.add_units_options(parser)
parser = pt.add_ftype_options(parser)
(options, args) = parser.parse_args()
files = args
return files,options
########################################################################
# MAIN
########################################################################
if __name__ == '__main__':
files,opt = parse_options()
X = opt.dX
# define local units
unit_params = pt.do_units_options(opt)
system_of_units = units.Set_SystemUnits_From_Params(unit_params)
G=ctes.GRAVITY.into(system_of_units)
H = ctes.HUBBLE.into(system_of_units)
HubbleParam = unit_params['HubbleParam']
rhoc = pow(H,2)*3/(8*pi*G)
rhoX = rhoc*X * unit_params['Omega0']
print "rhoX (code units, dX=%g)"%opt.dX,rhoX
# output system of units (the mass units is the hydrogen mass)
Unit_atom = ctes.PROTONMASS.into(units.cgs)*units.Unit_g
Unit_atom.set_symbol('atom')
out_units = units.UnitSystem('local',[units.Unit_cm,Unit_atom,units.Unit_s,units.Unit_K])
funit = system_of_units.convertionFactorTo(out_units.UnitDensity)
atime = 1.0
print " converting to physical units (a=%5.3f h=%5.3f)"%(atime,HubbleParam)
rhoXu = rhoX/atime**3*HubbleParam**2 *funit
print "rhoX (atom/cm^3)",rhoXu
print "log10rhoX (atom/cm^3)",log10(rhoXu)
###############################
# now, deal with files
###############################
if len(files)==0:
sys.exit()
file = files[0]
nb = Nbody(file,ftype=opt.ftype)
Rmax = nb.rxyz().max()
############################
# find rX using bissectrice
############################
rs = nb.rxyz()
def getRes(r):
nb_s = nb.selectc(rs<r)
M = sum(nb_s.mass)
V = 4/3.*pi*r**3
return M/V - rhoX
rX = bisection(getRes, 0.1, Rmax, args = (), xtol = 0.001, maxiter = 400)
nb_s = nb.selectc(nb.rxyz()<rX)
MX = sum(nb_s.mass)
V = 4/3.*pi*rX**3
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Msol,units.Unit_s,units.Unit_K])
fL = system_of_units.convertionFactorTo(out_units.UnitLength)
fM = system_of_units.convertionFactorTo(out_units.UnitMass)
print
print "Virial radius : r%d = %g [kpc/h]"%(X,rX*fL)
print "Virial mass : M%d = %g [Msol/h]"%(X,MX*fM)
print
print "Virial radius : r%d = %g [kpc]"%(X,rX*fL/HubbleParam)
print "Virial mass : M%d = %g [Msol]"%(X,MX*fM/HubbleParam)

Event Timeline

c4science · Help