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 <https://arxiv.org/abs/2210.01131>`_.


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

.. program-output:: ic_hernquist -h

``ic_plummer``
..............

.. program-output:: ic_plummer -h


``ic_nfw``
..........

.. program-output:: ic_nfw -h


``ic_nfw+plummer``
..................

.. program-output:: ic_nfw+plummer -h



``ic_homogeneous_box``
......................

.. program-output:: ic_homogeneous_box -h