Generating initial conditions from a distribution function¶

It is possible to generate initial conditions for ergodic spherical models by sampling a distribution function. To this end, pNbody provides a submodule called DF. The method is strongly inspired by the method described in Errani 2022.

The list of models currently implemented is given in the following table:

script name	model description
ic_plummer	isolated Plummer sphere
ic_nfw	isolated NFW model
ic_hernquist	isolated Hernquist model
ic_gen_2_slopes	isolated 2-slopes generic model
ic_nfw+plummer	NFW model combined with a Plummer sphere
ic_gen_2_slopes+plummer	2-slopes generic model combined with a Plummer sphere
ic_homogeneous_box	homogeneous box

The documentation for each these scripts with examples is given below.

For a given mass density profile provided by a python function, the corresponding distribution function is obtained using the Eddigton equation. By default, through the process, the potential, the cumulative mass, as well as the total mass are computed by numerical integration. Note that those quantities can also be given directly by the appropriate function if analytically defined as in the example given below:

# import modules
import numpy as np
from pNbody import *
from pNbody.DF import DistributionFunction
from pNbody.mass_models import plummer

# define the gravitational constant
G = 43000

# define some parameters
Mtot = 1e-5     # total Plummer mass
a    = 0.1      # Plummer scale length

Rmin = 1e-2     # minimal radius
Rmax = 10       # maximal radius

N    = 1e5      # number of particles
NR   = 1e4      # number of radius bins
NE   = 1e4      # number of energy bins
Ndraw= 1e6      # number of particles to draw at each loop


# define the density, potential, cumulative mass and total mass
fctDensity        = lambda x:plummer.Density(M=Mtot, a=a, r=x, G=G)
fctPotential      = lambda x:plummer.Potential(M=Mtot, a=a, r=x, G=G)
fctCumulativeMass = lambda x:plummer.CumulativeMass(M=Mtot, a=a, r=x, G=G)
TotalMass         = plummer.TotalMass(M=Mtot, a=a, G=G)

# create the distribution function object
DF = DistributionFunction(Rmin=Rmin,Rmax=Rmax,Rinf=np.inf,NR=NR,NE=NE,
                          fctDensity=fctDensity,
                          fctPotential=fctPotential,
                          fctCumulativeMass=fctCumulativeMass,
                          TotalMass=TotalMass,G=G)

Once the DF object is defined, we can compute it, i.e., effectively solving the Eddington equation, clean it, i.e., removing negative values and compute the maximum likelihood for a given radius:

DF.computeDF()
DF.clean()
DF.computeMaxLikelihood()

Next, we can draw N particles from the distribution and retrieve the sampled positions, velocities and masses:

DF.sample(N,Ndraw)

pos = DF.pos
vel = DF.vel
mass= DF.mass

Finally, it is possible to create an Nbody object:

nb = Nbody(status='new',p_name="plummer.hdf5",pos=pos,vel=vel,mass=mass,ftype='swift')

Scripts¶

`ic_hernquist`¶

usage: ic_hernquist [-h] [-o OUTPUTFILENAME] [-t FTYPE] [--irand INT]
                    [--ptype INT] [--plot] [-N INT] [-a A] [--rho0 RHO0]
                    [-b B] [--Rmin RMIN] [--Rmax RMAX] [--NE NE] [--NR NR]
                    [--Ndraw NDRAW] [--mass MASS] [-q]

generate an hernquist model at equilibrium

options:
  -h, --help         show this help message and exit
  -o OUTPUTFILENAME  Name of the output file
  -t FTYPE           Type of the output file
  --irand INT        random seed
  --ptype INT        particle type
  --plot             plot dynamical time
  -N INT             total number of particles
  -a A               a parameter
  --rho0 RHO0        density at critical radius
  -b B               b parameter
  --Rmin RMIN        Rmin
  --Rmax RMAX        Rmax
  --NE NE            number of energy bins
  --NR NR            number of radius bins
  --Ndraw NDRAW      number of particles to draw at each loop
  --mass MASS        particle mass in solar mass
  -q, --quiet        quiet mode

Examples:
--------
ic_hernquist -N 730000 --Rmin 1e-3 --Rmax 100 --rho0 0.1 -a 0.1  -t swift -o hernquist.hdf5

`ic_plummer`¶

usage: ic_plummer [-h] [-o OUTPUTFILENAME] [-t FTYPE] [--irand INT]
                  [--ptype INT] [--plot] [-N INT] [-a A] [--Mtot MTOT]
                  [--Rmin RMIN] [--Rmax RMAX] [--NE NE] [--NR NR]
                  [--Ndraw NDRAW] [--mass MASS] [-q]

generate a plummer model at equilibrium

options:
  -h, --help         show this help message and exit
  -o OUTPUTFILENAME  Name of the output file
  -t FTYPE           Type of the output file
  --irand INT        random seed
  --ptype INT        particle type
  --plot             plot dynamical time
  -N INT             total number of particles
  -a A               a parameter
  --Mtot MTOT        Total mass in 1e10 Msol
  --Rmin RMIN        Rmin
  --Rmax RMAX        Rmax
  --NE NE            number of energy bins
  --NR NR            number of radius bins
  --Ndraw NDRAW      number of particles to draw at each loop
  --mass MASS        particle mass in solar mass
  -q, --quiet        quiet mode

Examples:
--------
ic_plummer   -N 100000 --Mtot 1e-5 -a 0.1 --Rmax 10 -t swift -o plummer.hdf5

`ic_nfw`¶

usage: ic_nfw [-h] [-o OUTPUTFILENAME] [-t FTYPE] [--irand INT] [--ptype INT]
              [--plot] [-N INT] [--rs RS] [--rho0 RHO0] [--c C] [--M200 M200]
              [--Rmin RMIN] [--Rmax RMAX] [--NE NE] [--NR NR] [--Ndraw NDRAW]
              [--mass MASS] [-q]

generate an nfw model at equilibrium

options:
  -h, --help         show this help message and exit
  -o OUTPUTFILENAME  Name of the output file
  -t FTYPE           Type of the output file
  --irand INT        random seed
  --ptype INT        particle type
  --plot             plot dynamical time
  -N INT             total number of particles
  --rs RS            rs parameter
  --rho0 RHO0        density at critical radius
  --c C              concentration parameter
  --M200 M200        virial mass in Msol
  --Rmin RMIN        Rmin
  --Rmax RMAX        Rmax
  --NE NE            number of energy bins
  --NR NR            number of radius bins
  --Ndraw NDRAW      number of particles to draw at each loop
  --mass MASS        particle mass in solar mass
  -q, --quiet        quiet mode

Examples:
--------
ic_nfw       -N 730000 --Rmin 1e-3 --Rmax 100  --rho0 0.1   --rs 0.1    -t swift -o nfw.hdf5
ic_nfw       -N 730000 --Rmin 1e-3 --Rmax 100  --c    8     --M200 1e10 -t swift -o nfw.hdf5

`ic_gen_2_slopes`¶

usage: ic_gen_2_slopes [-h] [-o OUTPUTFILENAME] [-t FTYPE] [--irand INT]
                       [--plot] [-N INT] [--c C] [--M200 M200] [--rho0 RHO0]
                       [--rs RS] [--alpha ALPHA] [--beta BETA] [--Rmin RMIN]
                       [--Rmax RMAX] [--NE NE] [--NR NR] [--Ndraw NDRAW]
                       [--ptype INT] [--mass MASS] [--r_cutoff R_CUTOFF]
                       [--power_cutoff POWER_CUTOFF] [-q]

Generates a generalized 2-slopes model at equilibrium.

The generalized 2-slope model is parametrized either by:

  rs    : the scale radius, in kpc
  rho0  : the density at the scale radius, in Msol/kpc^3
  a     : the inner slope
  b     : the outer slope

or:

  M200  : the virial mass, in Msol
  c     : the concentration parameter
  a     : the inner slope
  b     : the outer slope

options:
  -h, --help            show this help message and exit
  -o OUTPUTFILENAME     Name of the output file
  -t FTYPE              Type of the output file
  --irand INT           Random seed
  --plot                Plot dynamical time
  -N INT                Total number of particles
  --c C                 concentration parameter
  --M200 M200           virial mass in Msol
  --rho0 RHO0           Gen2slopes : density at rs
  --rs RS               Gen2slopes : rs parameter
  --alpha ALPHA         Gen2slopes : alpha parameter
  --beta BETA           Gen2slopes : beta parameter
  --Rmin RMIN           Rmin
  --Rmax RMAX           Rmax
  --NE NE               number of energy bins
  --NR NR               number of radius bins
  --Ndraw NDRAW         number of particles to draw at each loop
  --ptype INT           particle type
  --mass MASS           particle mass in solar mass
  --r_cutoff R_CUTOFF   Density cutoff radius
  --power_cutoff POWER_CUTOFF
                        Power of the cutoff
  -q, --quiet           quiet mode

Examples:
--------
ic_gen_2_slopes --alpha 1 --beta 3 --Rmin 1e-3 --Rmax 30 --r_cutoff 30  --rho0 0.007 --r_s 0.5  --ptype 1--mass 1000 -t swift -o gen_2_slopes.hdf5

ic_gen_2_slopes --alpha 1 --beta 3 --Rmin 1e-3 --Rmax 30 --r_cutoff 30 --M200 1e9 --c 17  --mass 1000  -t swift -o gen_2_slopes.hdf5

`ic_nfw+plummer`¶

usage: ic_nfw+plummer [-h] [-o OUTPUTFILENAME] [-t FTYPE] [--irand INT]
                      [--plot] [-N INT] [--rs RS] [--rho0 RHO0] [--c C]
                      [--M200 M200] [-a A] [--Mtot MTOT] [--Rmin RMIN]
                      [--Rmax RMAX] [--NE NE] [--NR NR] [--Ndraw NDRAW]
                      [--ptype1 INT] [--ptype2 INT] [--mass1 MASS1]
                      [--mass2 MASS2] [-q]

generate an nfw+plummer model at equilibrium

options:
  -h, --help         show this help message and exit
  -o OUTPUTFILENAME  Name of the output file
  -t FTYPE           Type of the output file
  --irand INT        random seed
  --plot             plot dynamical time
  -N INT             total number of particles
  --rs RS            NFW rs parameter
  --rho0 RHO0        NFW density at critical radius
  --c C              concentration parameter
  --M200 M200        virial mass in Msol
  -a A               a parameter
  --Mtot MTOT        Total mass in solar mass
  --Rmin RMIN        Rmin
  --Rmax RMAX        Rmax
  --NE NE            number of energy bins
  --NR NR            number of radius bins
  --Ndraw NDRAW      number of particles to draw at each loop
  --ptype1 INT       particle type of component 1
  --ptype2 INT       particle type of component 2
  --mass1 MASS1      particle 1 mass in solar mass
  --mass2 MASS2      particle 2 mass in solar mass
  -q, --quiet        quiet mode

Examples:
--------
ic_nfw+plummer --Rmin 1e-3 --Rmax 30  --rho0 0.007 --rs 0.5   --Mtot 1e5 -a 0.1 --ptype1 4 --ptype2 1  --mass1 1000 --mass2 1000  -t swift -o nfw+plummer.hdf5
ic_nfw+plummer --Rmin 1e-3 --Rmax 100 --M200 1e11  --c  17    --Mtot 1e9 -a 1   --ptype1 4 --ptype2 1  -N 100000                  -t swift -o nfw+plummer.hdf5
ic_nfw+plummer --Rmin 1e-3 --Rmax 100 --M200 1e11  --c  17    --Mtot 1e9 -a 1   --ptype1 4 --ptype2 1  --mass1 1e4   --mass2 1e5  -t swift -o nfw+plummer.hdf5

`ic_gen_2_slopes+plummer`¶

usage: ic_gen_2_slopes+plummer [-h] [-o OUTPUTFILENAME] [-t FTYPE]
                               [--irand INT] [--plot] [-N INT] [--c C]
                               [--M200 M200] [--rho0 RHO0] [--rs RS]
                               [--alpha ALPHA] [--beta BETA] [-a A]
                               [--Mtot MTOT] [--Rmin RMIN] [--Rmax RMAX]
                               [--NE NE] [--NR NR] [--Ndraw NDRAW]
                               [--ptype1 INT] [--ptype2 INT] [--mass1 MASS1]
                               [--mass2 MASS2] [--r_cutoff R_CUTOFF]
                               [--power_cutoff POWER_CUTOFF] [-q]

Generates a generalized 2-slopes + Plummer model at equilibrium.

The model is composed of a  generalized 2-slopes model (the dark halo) as well as a Plummer model (the stellar counterpart).

The Plummer is parametrized by :

  Mtot : the total Plummer mass, in Msol 
  a    : the scale radius, equivalent to the half light radius when projected, in kpc

The generalized 2-slope model is parametrized either by:

  rs    : the scale radius, in kpc
  rho0  : the density at the scale radius, in Msol/kpc^3
  a     : the inner slope
  b     : the outer slope

or:

  M200  : the virial mass, in Msol
  c     : the concentration parameter
  a     : the inner slope
  b     : the outer slope

options:
  -h, --help            show this help message and exit
  -o OUTPUTFILENAME     Name of the output file
  -t FTYPE              Type of the output file
  --irand INT           Random seed
  --plot                Plot dynamical time
  -N INT                Total number of particles
  --c C                 concentration parameter
  --M200 M200           virial mass in Msol
  --rho0 RHO0           Gen2slopes : density at rs
  --rs RS               Gen2slopes : rs parameter
  --alpha ALPHA         Gen2slopes : alpha parameter
  --beta BETA           Gen2slopes : beta parameter
  -a A                  Plummer a parameter
  --Mtot MTOT           Plummer total mass in solar mass
  --Rmin RMIN           Rmin
  --Rmax RMAX           Rmax
  --NE NE               number of energy bins
  --NR NR               number of radius bins
  --Ndraw NDRAW         number of particles to draw at each loop
  --ptype1 INT          particle type of component 1
  --ptype2 INT          particle type of component 2
  --mass1 MASS1         particle 1 mass in solar mass
  --mass2 MASS2         particle 2 mass in solar mass
  --r_cutoff R_CUTOFF   Density cutoff radius
  --power_cutoff POWER_CUTOFF
                        Power of the cutoff
  -q, --quiet           quiet mode

Examples:
--------
ic_gen_2_slopes+plummer --alpha 1 --beta 3 --Rmin 1e-3 --Rmax 30 --r_cutoff 30  --rho0 0.007 --r_s 0.5 --Mtot 1e5 -a 0.1 --ptype1 4 --ptype2 1 --mass1 1000 --mass2 1000 -t swift -o gen_2_slopes+plummer.hdf5 

ic_gen_2_slopes+plummer --alpha 0 --beta 3 --Rmin 1e-3 --Rmax 30 --r_cutoff 30 --rho0 0.0008 --r_s 0.5 --Mtot 1e5 -a 0.1 --ptype1 4 --ptype2 1 --mass1 125 --mass2 125 -t swift -o 2s0_3_0008_05+plummer.hdf5 

ic_gen_2_slopes+plummer --alpha 1 --beta 3 --Rmin 1e-3 --Rmax 30 --r_cutoff 30 --M200 1e9 --c 17 --Mtot 1e5 -a 0.1 --ptype1 4 --ptype2 1 --mass1 1000 --mass2 10000 -t swift -o gen_2_slopes+plummer.hdf5 
ic_gen_2_slopes+plummer --alpha 0 --beta 3 --Rmin 1e-3 --Rmax 30 --r_cutoff 30 --M200 1e9 --c 17 --Mtot 1e5 -a 0.1 --ptype1 4 --ptype2 1 --mass1 1000 --mass2 10000 -t swift -o gen_2_slopes+plummer.hdf5 

ic_gen_2_slopes+plummer --alpha 1 --beta 3 --Rmin 1e-3 --Rmax 30 --r_cutoff 30 --M200 1e9 --c 17 --Mtot 1e5 -a 0.1 --ptype1 4 --ptype2 1 -N 1000000 -t swift -o gen_2_slopes+plummer.hdf5 
ic_gen_2_slopes+plummer --alpha 0 --beta 3 --Rmin 1e-3 --Rmax 30 --r_cutoff 30 --M200 1e9 --c 17 --Mtot 1e5 -a 0.1 --ptype1 4 --ptype2 1 -N 1000000 -t swift -o gen_2_slopes+plummer.hdf5

`ic_homogeneous_box`¶

usage: ic_homogeneous_box [-h] [-o OUTPUTFILENAME] [-t FTYPE] [--lJ FLOAT]
                          [--irand INT] [--hydro] [--ptype INT] [-N INT] [-q]

generate initial an homogeneous box

options:
  -h, --help         show this help message and exit
  -o OUTPUTFILENAME  Name of the output file
  -t FTYPE           Type of the output file
  --lJ FLOAT         Jeans wavelength
  --irand INT        random seed
  --hydro            consider gas
  --ptype INT        particle type
  -N INT             total number of particles
  -q, --quiet        quiet mode

Examples:
--------
ic_homogeneous_box -N 1024                                             -o output.hdf5
ic_homogeneous_box -N 1024 --lJ=0.5                                    -o output.hdf5
ic_homogeneous_box -N 1024 --lJ=0.5 --hydro                            -o output.hdf5
ic_homogeneous_box -N 1024 --lJ=0.5 --hydro -t swift                   -o output.hdf5
ic_homogeneous_box -N 1024 --lJ=0.5 --hydro -t swift --hydro           -o output.hdf5
ic_homogeneous_box -N 1024 --lJ=0.5 --hydro -t swift --hydro --ptype 2 -o output.hdf5

Generating initial conditions from a distribution function¶

Scripts¶

`ic_hernquist`¶

`ic_plummer`¶

`ic_nfw`¶

`ic_gen_2_slopes`¶

`ic_nfw+plummer`¶

`ic_gen_2_slopes+plummer`¶

`ic_homogeneous_box`¶

Table of Contents

Previous topic

Next topic

This Page

Generating initial conditions from a distribution function¶

Scripts¶

ic_hernquist¶

ic_plummer¶

ic_nfw¶

ic_gen_2_slopes¶

ic_nfw+plummer¶

ic_gen_2_slopes+plummer¶

ic_homogeneous_box¶

`ic_hernquist`¶

`ic_plummer`¶

`ic_nfw`¶

`ic_gen_2_slopes`¶

`ic_nfw+plummer`¶

`ic_gen_2_slopes+plummer`¶

`ic_homogeneous_box`¶