Homec4science
Diffusion Gtools (master) Browse

gtools/Gtools/pscripts/pget_Luminosity
77f210e77b30master

Last Change
Raw File

pget_Luminosity
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("--tnow",
action="store",
dest="tnow",
type="float",
default = None,
help="desired time",
metavar=" FLOAT")
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()
if len(files)==0:
sys.exit()
file = files[0]
nb = Nbody(file,ftype=opt.ftype)
unit_params = pt.do_units_options(opt)
HubbleParam = unit_params['HubbleParam']
nb.set_local_system_of_units(params=unit_params)
nb = nb.select(1)
Rmax = nb.rxyz().max()
############################
# find tR using bissectrice
############################
lmax=90
nb.L = nb.Luminosity(tnow=opt.tnow)
Ltot = sum(nb.L)
rs = nb.rxyz()
def getL(r):
nb_s = nb.selectc(rs<r)
L = sum(nb_s.L)
return lmax - 100*L/Ltot
Rt = bisection(getL, 0.0, 50, args = (), xtol = 0.1, maxiter = 400)
nb_s = nb.selectc(rs<Rt)
L = sum(nb_s.L)
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Msol,units.Unit_s,units.Unit_K])
fL = nb.localsystem_of_units.convertionFactorTo(out_units.UnitLength)
print
print "tidal radius : Rt = %g [kpc/h]"%(Rt*fL)
print "tidal radius : Rt = %g [kpc]"%(Rt*fL/HubbleParam)
print "Luminosity in tidal radius : L = %g [1e6 Lsol]"%(L/1e6)
print "Total luminosity : Ltot = %g [1e6 Lsol]"%(Ltot/1e6)
c4science · Help