Page MenuHomec4science
Files F121736172

pCylindricalProfile
No OneTemporary
Actions

File Metadata

Created
Sun, Jul 13, 11:45

pCylindricalProfile
View Options

#!/usr/bin/env python
from optparse import OptionParser
import Ptools as pt
from pNbody import *
from pNbody import units
from pNbody import ctes
from pNbody import thermodyn
import string
from scipy import optimize
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_postscript_options(parser)
parser = pt.add_ftype_options(parser)
parser = pt.add_reduc_options(parser)
parser = pt.add_center_options(parser)
parser = pt.add_select_options(parser)
parser = pt.add_cmd_options(parser)
parser = pt.add_display_options(parser)
parser = pt.add_info_options(parser)
parser = pt.add_limits_options(parser)
parser = pt.add_log_options(parser)
parse = pt.add_units_options(parser)
parser.add_option("--rmax",
action="store",
dest="rmax",
type="float",
default = 50.,
help="max radius of bins",
metavar=" FLOAT")
parser.add_option("--nr",
action="store",
dest="nr",
type="int",
default = 32,
help="number of bins in r",
metavar=" INT")
parser.add_option("--forceComovingIntegrationOn",
action="store_true",
dest="forceComovingIntegrationOn",
default = False,
help="force the model to be in in comoving integration")
parser.add_option("--curve",
action="store_true",
dest="curve",
default = False,
help="draw curve instead of histogram")
parser.add_option("-o",
action="store",
type="string",
dest="o",
default = 'sdens',
help="value to plot")
parser.add_option("--AdaptativeSoftenning",
action="store_true",
dest="AdaptativeSoftenning",
default = False,
help="AdaptativeSoftenning")
parser.add_option("--ErrTolTheta",
action="store",
dest="ErrTolTheta",
type="float",
default = 0.5,
help="Error tolerance theta",
metavar=" FLOAT")
parser.add_option("--eps",
action="store",
dest="eps",
type="float",
default = 0.1,
help="smoothing length",
metavar=" FLOAT")
parser.add_option("--nt",
action="store",
dest="nt",
type="int",
default = 32,
help="number of division in theta",
metavar=" INT")
(options, args) = parser.parse_args()
pt.check_files_number(args)
files = args
return files,options
#######################################
# MakePlot
#######################################
def MakePlot(files,opt):
# some inits
colors = pt.Colors(n=len(files))
datas = []
# read files
for file in files:
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)
# define output units
# nb.ToPhysicalUnits()
if opt.forceComovingIntegrationOn:
nb.setComovingIntegrationOn()
################
# apply options
################
nb = pt.do_reduc_options(nb,opt)
nb = pt.do_select_options(nb,opt)
nb = pt.do_center_options(nb,opt)
nb = pt.do_cmd_options(nb,opt)
nb = pt.do_info_options(nb,opt)
nb = pt.do_display_options(nb,opt)
################
# some info
################
print "---------------------------------------------------------"
nb.localsystem_of_units.info()
nb.ComovingIntegrationInfo()
print "---------------------------------------------------------"
# grid division
rc = 1
f = lambda r:log(r/rc+1.)
fm = lambda r:rc*(exp(r)-1.)
###############################
# compute physical quantities
###############################
##########################################
# Cylindrical_2drt_Grid
##########################################
if opt.o=='sdens':
G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=1,g=f,gm=fm)
x,t = G.get_rt()
y = G.get_SurfaceDensityMap(nb)
# comoving conversion
if nb.isComovingIntegrationOn():
print " converting to physical units (a=%5.3f h=%5.3f)"%(nb.atime,nb.hubbleparam)
x = x*nb.atime/nb.hubbleparam # length conversion
y = y/nb.atime**2*nb.hubbleparam # surface density conversion
# output units
out_units_x = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
out_units_y = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
fx = nb.localsystem_of_units.convertionFactorTo(out_units_x.UnitLength)
fy = nb.localsystem_of_units.convertionFactorTo(out_units_y.UnitSurfaceDensity)
x = x*fx
y = y*fy
# set labels
xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
ylabel = r'$\rm{Surface\,\,Density}\,\left[ M_{\odot}/kpc^2 \right]$'
x,y = pt.CleanVectorsForLogX(x,y)
x,y = pt.CleanVectorsForLogY(x,y)
datas.append(pt.DataPoints(x, y, color=colors.get(),label=None,tpe='points'))
if opt.o=='mass':
G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=1,g=f,gm=fm)
x,t = G.get_rt()
y = G.get_MassMap(nb)
y = sum(y,axis=1)
# comoving conversion
if nb.isComovingIntegrationOn():
print " converting to physical units (a=%5.3f h=%5.3f)"%(nb.atime,nb.hubbleparam)
x = x*nb.atime/nb.hubbleparam # length conversion
y = y/nb.hubbleparam # mass conversion
# output units
out_units_x = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
out_units_y = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
fx = nb.localsystem_of_units.convertionFactorTo(out_units_x.UnitLength)
fy = nb.localsystem_of_units.convertionFactorTo(out_units_y.UnitDensity)
x = x*fx
y = y*fy
# set labels
xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
ylabel = r'$\rm{Mass}\,\left[ M_{\odot} \right]$'
x,y = pt.CleanVectorsForLogX(x,y)
x,y = pt.CleanVectorsForLogY(x,y)
datas.append(pt.DataPoints(x, y, color=colors.get(),label=None,tpe='points'))
if opt.o=='imass':
G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=1,g=f,gm=fm)
x,t = G.get_rt()
y = G.get_MassMap(nb)
y = sum(y,axis=1)
y = add.accumulate(y)
# comoving conversion
if nb.isComovingIntegrationOn():
print " converting to physical units (a=%5.3f h=%5.3f)"%(nb.atime,nb.hubbleparam)
x = x*nb.atime/nb.hubbleparam # length conversion
y = y/nb.hubbleparam # mass conversion
# output units
out_units_x = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
out_units_y = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
fx = nb.localsystem_of_units.convertionFactorTo(out_units_x.UnitLength)
fy = nb.localsystem_of_units.convertionFactorTo(out_units_y.UnitDensity)
x = x*fx
y = y*fy
# set labels
xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
ylabel = r'$\rm{Mass}\,\left[ M_{\odot} \right]$'
x,y = pt.CleanVectorsForLogX(x,y)
x,y = pt.CleanVectorsForLogY(x,y)
datas.append(pt.DataPoints(x, y, color=colors.get(),label=None,tpe='points'))
if opt.o=='sigmaz':
# comoving conversion
if nb.isComovingIntegrationOn():
print " converting to physical units (a=%5.3f h=%5.3f)"%(nb.atime,nb.hubbleparam)
x0 = nb.pos
v0 = nb.vel
Hubble = ctes.HUBBLE.into(nb.localsystem_of_units)
pars = {"Hubble":Hubble,"HubbleParam":nb.hubbleparam,"OmegaLambda":nb.omegalambda,"Omega0":nb.omega0}
a = nb.atime
Ha = cosmo.Hubble_a(a,pars=pars)
nb.pos = x0*nb.atime/nb.hubbleparam # length conversion
nb.vel = v0*sqrt(a) + x0*Ha*a # velocity conversion
# output units
out_units_x = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
out_units_y = units.UnitSystem('local',[units.Unit_km,units.Unit_Ms,units.Unit_s,units.Unit_K])
fx = nb.localsystem_of_units.convertionFactorTo(out_units_x.UnitLength)
fy = nb.localsystem_of_units.convertionFactorTo(out_units_y.UnitVelocity)
nb.pos = nb.pos*fx
nb.vel = nb.vel*fy
G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=1,g=f,gm=fm)
x,t = G.get_rt()
y = G.get_SigmaValMap(nb,nb.Vz())
y = sum(y,axis=1)
# set labels
xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
ylabel = r'$\sigma_z\,\left[ \rm{km}/\rm{s} \right]$'
x,y = pt.CleanVectorsForLogX(x,y)
x,y = pt.CleanVectorsForLogY(x,y)
datas.append(pt.DataPoints(x, y, color=colors.get(),label=None,tpe='points'))
if opt.o=='vcirc':
G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=opt.nt,z=0,g=f,gm=fm)
x,t = G.get_rt()
# compute tree
nb.getTree(force_computation=True,ErrTolTheta=opt.ErrTolTheta)
# compute accel
Accx,Accy,Accz = G.get_AccelerationMap(nb,eps=opt.eps,UseTree=True,AdaptativeSoftenning=opt.AdaptativeSoftenning)
Ar = sqrt(Accx**2+Accy**2)
Ar = sum(Ar,1)/opt.nt
dPhit = Ar
y = libdisk.Vcirc(x,dPhit)
# comoving conversion
if nb.isComovingIntegrationOn():
print " converting to physical units (a=%5.3f h=%5.3f)"%(nb.atime,nb.hubbleparam)
x = x*nb.atime/nb.hubbleparam # length conversion
#y = y/nb.hubbleparam # velocity
print "!!! this is not defined for comobiles coords, please, do it !!!"
sys.exit()
# output units
out_units_x = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
out_units_y = units.UnitSystem('local',[units.Unit_km,units.Unit_Ms,units.Unit_s,units.Unit_K])
fx = nb.localsystem_of_units.convertionFactorTo(out_units_x.UnitLength)
fy = nb.localsystem_of_units.convertionFactorTo(out_units_y.UnitVelocity)
x = x*fx
y = y*fy
# set labels
xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
ylabel = r'$\rm{Circular\,Velocity}\,\left[ \rm{km}/\rm{s} \right]$'
x,y = pt.CleanVectorsForLogX(x,y)
x,y = pt.CleanVectorsForLogY(x,y)
datas.append(pt.DataPoints(x, y, color=colors.get(),label=None,tpe='points'))
if opt.o=='Lum':
# comoving conversion
if nb.isComovingIntegrationOn():
print " converting to physical units (a=%5.3f h=%5.3f)"%(nb.atime,nb.hubbleparam)
x0 = nb.pos
v0 = nb.vel
Hubble = ctes.HUBBLE.into(nb.localsystem_of_units)
pars = {"Hubble":Hubble,"HubbleParam":nb.hubbleparam,"OmegaLambda":nb.omegalambda,"Omega0":nb.omega0}
a = nb.atime
Ha = cosmo.Hubble_a(a,pars=pars)
nb.pos = x0*nb.atime/nb.hubbleparam # length conversion
nb.vel = v0*sqrt(a) + x0*Ha*a # velocity conversion
# output units
out_units_x = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
out_units_y = units.UnitSystem('local',[units.Unit_km,units.Unit_Ms,units.Unit_s,units.Unit_K])
fx = nb.localsystem_of_units.convertionFactorTo(out_units_x.UnitLength)
fy = nb.localsystem_of_units.convertionFactorTo(out_units_y.UnitVelocity)
nb.pos = nb.pos*fx
nb.vel = nb.vel*fy
G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=1,g=f,gm=fm)
x,t = G.get_rt()
y = reshape(G.get_SurfaceValMap(nb,nb.luminosity_spec()*1e10),opt.nr) # in Lsol/(unitlength)
print "the surface density units is wrong here...!!!"
# set labels
xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
ylabel = r'$\rm{Luminosity}\left[ L_{\odot}/\rm{pc}^2 \right]$'
if opt.log:
x,y = pt.CleanVectorsForLogX(x,y)
x,y = pt.CleanVectorsForLogY(x,y)
datas.append(pt.DataPoints(x, y, color=colors.get(),label=None,tpe='points'))
##################
# plot all
##################
for d in datas:
pt.plot(d.x,d.y,color=d.color)
# set limits and draw axis
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
pt.SetAxis(xmin,xmax,ymin,ymax,opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
pt.grid(False)
########################################################################
# MAIN
########################################################################
if __name__ == '__main__':
files,opt = parse_options()
pt.InitPlot(files,opt)
pt.pcolors
MakePlot(files,opt)
pt.EndPlot(files,opt)

Event Timeline

c4science · Help