Page MenuHomec4science
Files F91125265

pget_Luminosity
No OneTemporary
Actions

File Metadata

Created
Fri, Nov 8, 04:18

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")
parser.add_option("--forceComovingIntegrationOn",
action="store_true",
dest="forceComovingIntegrationOn",
default = False,
help="force the model to be in in comoving integration")
parser.add_option("--forceComovingIntegrationOff",
action="store_true",
dest="forceComovingIntegrationOff",
default = False,
help="force the model not to be in in comoving integration")
parse = pt.add_units_options(parser)
parser = pt.add_ftype_options(parser)
parser = pt.add_cmd_options(parser)
parser = pt.add_center_options(parser)
parser = pt.add_select_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)
################
# units
################
# define local units
unit_params = pt.do_units_options(opt)
nb.set_local_system_of_units(params=unit_params)
HubbleParam = unit_params['HubbleParam']
# define output units
# nb.ToPhysicalUnits()
if opt.forceComovingIntegrationOn:
nb.setComovingIntegrationOn()
if opt.forceComovingIntegrationOff:
nb.setComovingIntegrationOff()
################
# apply options
################
nb = pt.do_select_options(nb,opt)
nb = pt.do_cmd_options(nb,opt)
nb = pt.do_center_options(nb,opt)
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)

Event Timeline

c4science · Help