Page MenuHomec4science
Files F116253370

pSphericalProfile
No OneTemporary
Actions

File Metadata

Created
Thu, Jun 5, 20:00

pSphericalProfile
View Options

#!/usr/bin/python3
###########################################################################################
# package: Gtools
# file: pSphericalProfile
# brief:
# copyright: GPLv3
# Copyright (C) 2019 EPFL (Ecole Polytechnique Federale de Lausanne)
# LASTRO - Laboratory of Astrophysics of EPFL
# author: Yves Revaz <yves.revaz@epfl.ch>
#
# This file is part of Gtools.
###########################################################################################
from optparse import OptionParser
import Ptools as pt
from pNbody import *
from pNbody import units
from pNbody import ctes
from pNbody import libgrid, cosmo, libdisk
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)
parser = pt.add_units_options(parser)
parser = pt.add_legend_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(
"--forceComovingIntegrationOff",
action="store_true",
dest="forceComovingIntegrationOff",
default=False,
help="force the model not 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='density',
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")
(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()
if opt.forceComovingIntegrationOff:
nb.setComovingIntegrationOff()
################
# 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 specific stuffs
nb.rho = nb.mass/nb.rsp**3
#nb = nb.selectc(nb.rho<1e-6)
################
# some info
################
print("---------------------------------------------------------")
nb.localsystem_of_units.info()
nb.ComovingIntegrationInfo()
print("---------------------------------------------------------")
# grid division
rc = 1
def f(r): return np.log(r / rc + 1.)
def fm(r): return rc * (np.exp(r) - 1.)
###############################
# compute physical quantities
###############################
##########################################
# Spherical_1d_Grid
##########################################
if opt.o == 'density':
G = libgrid.Spherical_1d_Grid(
rmin=0, rmax=opt.rmax, nr=opt.nr, g=f, gm=fm)
x = G.get_r()
y = G.get_DensityMap(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**3 * nb.hubbleparam**2 # density conversion
# output units
out_units_x = units.UnitSystem(
'local', [units.Unit_kpc, units.Unit_Ms, units.Unit_Myr, units.Unit_K])
Unit_atom = ctes.PROTONMASS.into(units.cgs) * units.Unit_g
Unit_atom.set_symbol('atom')
out_units_y = units.UnitSystem(
'local', [units.Unit_cm, Unit_atom, 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.UnitDensity)
x = x * fx
y = y * fy
# set labels
xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
ylabel = r'$\rm{Density}\,\left[ \rm{atom/cm^3} \right]$'
x, y = pt.CleanVectorsForLogX(x, y)
x, y = pt.CleanVectorsForLogY(x, y)
datas.append(
pt.DataPoints(
x,
y,
color=colors.get(),
label=file,
tpe='points'))
if opt.o == 'mass':
G = libgrid.Spherical_1d_Grid(
rmin=0, rmax=opt.rmax, nr=opt.nr, g=f, gm=fm)
x = G.get_r()
y = G.get_MassMap(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.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=file,
tpe='points'))
if opt.o == 'imass':
G = libgrid.Spherical_1d_Grid(
rmin=0, rmax=opt.rmax, nr=opt.nr, g=f, gm=fm)
x = G.get_r()
y = G.get_MassMap(nb)
y = np.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=file,
tpe='points'))
if ((opt.o == 'sigmar') or (opt.o == "sigmat") or (opt.o == "sigmatheta") or (opt.o == "sigmaphi") or (opt.o == "sigmaz") or (opt.o == "beta") or (opt.o == "Vt") or (opt.o == "Vr") or (opt.o == "Vtheta") or (opt.o == "Vphi")):
# 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 * np.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.Spherical_1d_Grid(
rmin=0, rmax=opt.rmax, nr=opt.nr, g=f, gm=fm)
x = G.get_r()
if opt.o == "beta":
sigmaphi = G.get_SigmaValMap(nb, nb.Vphi())
sigmatheta = G.get_SigmaValMap(nb, nb.Vtheta())
st = np.sqrt(sigmaphi**2 + sigmatheta**2)
sr = G.get_SigmaValMap(nb, nb.Vr())
y = 1-((st**2)/(2.*sr**2))
ylabel = r'$\beta$'
else:
if opt.o == 'sigmaz':
y = G.get_SigmaValMap(nb, nb.Vz())
ylabel = r'$\sigma_z\,\left[ \rm{km}/\rm{s} \right]$'
if opt.o == 'sigmar':
y = G.get_SigmaValMap(nb, nb.Vr())
ylabel = r'$\sigma_r\,\left[ \rm{km}/\rm{s} \right]$'
if opt.o == 'sigmat':
sigmaphi = G.get_SigmaValMap(nb, nb.Vphi())
sigmatheta = G.get_SigmaValMap(nb, nb.Vtheta())
y = np.sqrt(sigmaphi**2 + sigmatheta**2)
ylabel = r'$\sigma_t\,\left[ \rm{km}/\rm{s} \right]$'
if opt.o == 'sigmatheta':
y = G.get_SigmaValMap(nb, nb.Vtheta())
ylabel = r'$\sigma_\theta\,\left[ \rm{km}/\rm{s} \right]$'
if opt.o == 'sigmaphi':
y = G.get_SigmaValMap(nb, nb.Vphi())
ylabel = r'$\sigma_\phi\,\left[ \rm{km}/\rm{s} \right]$'
if opt.o == 'Vt':
y = G.get_MeanValMap(nb, nb.Vt())
ylabel = r'$V_t\,\left[ \rm{km}/\rm{s} \right]$'
if opt.o == 'Vr':
y = G.get_MeanValMap(nb, nb.Vr())
ylabel = r'$V_r\,\left[ \rm{km}/\rm{s} \right]$'
if opt.o == 'Vtheta':
y = G.get_MeanValMap(nb, nb.Vtheta())
ylabel = r'$V_\theta\,\left[ \rm{km}/\rm{s} \right]$'
if opt.o == 'Vphi':
y = G.get_MeanValMap(nb, nb.Vphi())
ylabel = r'$V_\phi\,\left[ \rm{km}/\rm{s} \right]$'
# set xlabels
xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
#x, y = pt.CleanVectorsForLogX(x, y)
#x, y = pt.CleanVectorsForLogY(x, y)
datas.append(
pt.DataPoints(
x,
y,
color=colors.get(),
label=file,
tpe='points'))
if opt.o == 'vcirc':
G = libgrid.Spherical_1d_Grid(rmin=0, rmax=opt.rmax, nr=opt.nr, g=f, gm=fm)
if nb.isComovingIntegrationOn():
# enventually correct pos and mass from h and a
fp = nb.ConversionFactor(units=None,mode='pos')
fm = nb.ConversionFactor(units=None,mode='mass')
else:
fp = 1
fm = 1
# r
r = G.get_r()*fp
# M(r)
M = G.get_MassMap(nb)
M = np.add.accumulate(M)*fm
# Newton theorem
G=ctes.GRAVITY.into(nb.localsystem_of_units)
vc = np.sqrt(G*M/r)
x = r
y = vc
# 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'$V_{\rm c}\,\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=file,
tpe='points'))
if opt.o == 'tdyn':
G = libgrid.Spherical_1d_Grid(rmin=0, rmax=opt.rmax, nr=opt.nr, g=f, gm=fm)
if nb.isComovingIntegrationOn():
# enventually correct pos and mass from h and a
fp = nb.ConversionFactor(units=None,mode='pos')
fm = nb.ConversionFactor(units=None,mode='mass')
else:
fp = 1
fm = 1
# r
r = G.get_r()*fp
# M(r)
M = G.get_MassMap(nb)
M = np.add.accumulate(M)*fm
# Newton theorem
G=ctes.GRAVITY.into(nb.localsystem_of_units)
vc = np.sqrt(G*M/r)
x = r
y = 2*np.pi*r/vc
# 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_Myr, units.Unit_K])
fx = nb.localsystem_of_units.convertionFactorTo(
out_units_x.UnitLength)
fy = nb.localsystem_of_units.convertionFactorTo(
out_units_y.UnitTime)
x = x * fx
y = y * fy
# set labels
xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
ylabel = r'$t_{\rm dyn}\,\left[ \rm{Myr} \right]$'
x, y = pt.CleanVectorsForLogX(x, y)
x, y = pt.CleanVectorsForLogY(x, y)
datas.append(
pt.DataPoints(
x,
y,
color=colors.get(),
label=file,
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)
if opt.legend:
pt.LegendFromDataPoints(datas, opt.legend_loc)
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