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.
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_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_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
