Homec4science
Diffusion pNbody (master) Browse

pNbody/Doc/rst/GeneratingMassProfiles.rst
20b0c82b73acmaster

Last Change
Raw File

GeneratingMassProfiles.rst
View Options

Generating mass profiles
************************
.. currentmodule:: pNbody.ic
A serie of routines dedictated to the generation of mass profiles
are provided by :doc:`IcModule`.
They can be devided according to their geometry:
.. table:: cubic distribution
=============== ===============
function name description
=============== ===============
box particles in a box
=============== ===============
.. table:: axi-symetrical (disk) distribution
=============== ===============
function name description
=============== ===============
homodisk homogeneous disk
kuzmin kuzmin disk
expd exponnential disk
miyamoto_nagai Miyamoto-Nagai disk
=============== ===============
.. table:: spherical distribution
=============== ===============
function name description
=============== ===============
homosphere homogeneous sphere
plummer phlummer sphere
nfw Navarro-Frenk-White profile
hernquist Hernquist profile
dl2 .. math:: \rho \sim (1-\epsilon (r/r_{\rm{max}})^2)
isothm Isothermal profile
pisothm Pseudo-Isothermal profile
shell shell
burkert Burkert profile
nfwg Modified Navarro-Frenk-White profile
generic2c generic two slope model
generic_alpha .. math:: \rho \sim (r+e)^a
generic_Mr spherical generic model
=============== ===============
Finally, some routines generate special object used
for rendering:
.. table:: primer objects
=============== ===============
function name description
=============== ===============
line a simple line
square a square
circle a circle
grid a 2d grid
cube a cube
sphere a sphere
=============== ===============
Each of these routines returns an **Nbody** object.
Examples:
=========
>>> from pNbody import ic
>>> nb=ic.box(1000,1,1,1)
>>> nb.info()
-----------------------------------
particle file : ['box.dat']
ftype : 'Nbody_binary'
mxntpe : 6
nbody : 1000
nbody_tot : 1000
npart : [1000, 0, 0, 0, 0, 0]
npart_tot : [1000, 0, 0, 0, 0, 0]
mass_tot : 0.9999907
byteorder : 'little'
pio : 'no'
len pos : 1000
pos[0] : array([ 0.78348911, 0.03045347, 0.34180355], dtype=float32)
pos[-1] : array([ 0.14677997, -0.05187676, 0.11189017], dtype=float32)
len vel : 1000
vel[0] : array([ 0., 0., 0.], dtype=float32)
vel[-1] : array([ 0., 0., 0.], dtype=float32)
len mass : 1000
mass[0] : 0.001
mass[-1] : 0.001
len num : 1000
num[0] : 0
num[-1] : 999
len tpe : 1000
tpe[0] : 0
tpe[-1] : 0
tnow : 0.0
label : binary
dt : 0.0
>>>
>>> nb.display()
Generic spherical distributions:
================================
In spherical cases, its possible to generate any mass profile, using a
generic function (:func:`generic_Mr`). The function takes as argument a vector that
defines the cumulative mass as a function of the radius:
.. math:: M(r)=4\pi\int_0^r \rho(r')r'^2 dr'
for example, an homogeneous sphere may be obtained by:
>>> from pNbody import ic
>>> from numpy import *
>>> rmax=10
>>> dr=0.01
>>> n = 10000
>>> rs = arange(0,rmax,dr)
>>> Mrs = rs**3
>>> nb = ic.generic_Mr(n,rmax=rmax,R=rs,Mr=Mrs,ftype='gadget')
The verctor ``rs`` is not necessary linear. This is usefull
to avoid resolution problems, if for example, the mass profile is steep
in some regions. There, intervals between radial bins may be reduced.
In this more funny example, we generate a radial profile corresponding to the following function:
.. math:: \rho(r) \sim \frac{1}{2}\left[1+\cos(4r) \right] \exp(-r)
As before, we use the function :func:`generic_Mr`::
>>> from pNbody import ic
>>> from numpy import *
>>> n = 1000000
>>> rmax = 10.
>>> dr = 0.01
>>> rs = arange(0,rmax,dr)
>>> rho = 0.5*(1+cos(4*rs) )*exp(-rs)
>>> Mrs = add.accumulate(rho*rs**2 * dr)
>>> nb = ic.generic_Mr(n,rmax=rmax,R=rs,Mr=Mrs,ftype='gadget')
It is possible to check the result by extractig the accumulated mass directely from the **Nbody** object,
using :doc:`LibgridModule`::
>>> from pNbody import libgrid
>>> # create a grid object and extract density and mass profile
>>> G = libgrid.Spherical_1d_Grid(rmin=0,rmax=rmax,nr=512,g=None,gm=None)
>>> r = G.get_r()
>>> rho = G.get_DensityMap(nb)
>>> mr = add.accumulate(G.get_MassMap(nb))
>>> # plot results
>>> import pylab as pt
>>> pt.plot(rs,Mrs)
>>> pt.plot(r,mr)
>>> pt.xlabel("radius")
>>> pt.ylabel("accumulated mass")
>>> pt.show()
If everything goes well, we obtain a nice correspondance between the imposed and
the computed mass profile, as seen in the graph below.
.. image:: ../images/mass_profile.png
Generic cubic distributions:
============================
Its possible generate a cubic box with the density varying along one axix. This is obtained
by using the :func:`generic_Mx` function. As for :func:`generic_Mr`, we have to define
the variation of the cumulative mass, in this case, the x axis. In the following example,
we set a two level transition along the x axis, following the function:
.. math:: \rho(x) = \rho_1 + \frac{\rho_2-\rho_1}{ (1+\exp(-2(x-x_1)/\sigma))(1+\exp(2(x-x_2)/\sigma)) }
This gives:
>>> from pNbody import ic
>>> from numpy import *
>>> xs = arange(0,1,1./10000.)
>>> n = 100000
>>> rho1 = 1
>>> rho2 = 10
>>> x1 = 0.25
>>> x2 = 0.75
>>> sigma = 0.025
>>> rho = rho1 + (rho2-rho1)/( (1+exp(-2.*(xs-x1)/sigma))*(1+exp(2.*(xs-x2)/sigma)))
>>> Mxs = add.accumulate(rho)
>>> nb = ic.generic_Mx(n,1.0,x=xs,Mx=Mxs,name='tbox.dat',ftype='gadget')
Let's display now the model:
>>> nb.translate([-0.5,-0.5,-0.5])
>>> nb.show(view='xy',size=(0.6,0.6))
The two density levels are well seen in projection.
.. image:: ../images/two_levels.png
Spherical distribution with variable resolution:
================================================
In some circumstances, it may be insterested to generate a profile, with
a variable resolution. In all previous examples, the mass of the particles where
identical. Here, we show a case where the mass of the particles varies with the radius.
This allows to increase the number of particles in the central regrions, where resolution
is needed, and to decrease it in the outer regions.
>>> nb_halo = GenericDistribution(n_halo,dmdr_fct_halo,(rs_halo,gamma_halo),rmax_halo,ftype='gadget')
c4science · Help