.. _mockimgs_preparation-label:

Setup the N-body model
======================

We describe here the quantities that must be stored in the N-body model in order
to compute a surface brightness image.
We assume here a ``swift`` hdf5 format, while other formats can be used.

The N-body model must contain the following fields:

* positions for each stellar particle, in :math:`\rm{kpc}` (``nb.pos``)
* velocities for each stellar particle, in :math:`\rm{km/s}` (``nb.vel``)
* mass for each stellar particle, in :math:`10^{10}\,\rm{M_\odot}` (``nb.mass``). The mast must corresponds to the initial SSP mass.
* age for each stellar particle, in :math:`\rm{Gyr}` (``nb.age``)
  An old stellar particle will have a typical age of :math:`12\,\rm{Gyr}` and a recently formed particles of :math:`2\,\rm{Gyr}`

* metallicity for each stellar particle ([M/H])  (``nb.mh``)

  .. math::
  
    [M/H] = \log_{10}   \frac{M_{\rm{M}}/M_{\rm{H}}}{ \left( M_{\rm{M}}/M_{\rm{H}}\right)_{\odot}  } 
  
  where :math:`M_{\rm{M}}/M_{\rm{H}}` is the mass fraction of metals and usually denoted by :math:`Z`
  and :math:`\left( M_{\rm{M}}/M_{\rm{H}}\right)_{\odot}` is the solar metal fraction and is set to 0.02.

* an estimation of the size for each stellar particle, in :math:`\rm{kpc}` (``nb.rsp``). 


.. warning::

  All quantities that must be provided must be in proper coorinates (not in comoving ones) and must
  be free of Hubble parameter (:math:`h`)!
  
  
  
The following lines show how to get a snapshot of the appropriate format using 
`pNbody <http://obswww.unige.ch/~revaz/pNbody>`_::

  import numpy as np
  from pNbody import *
  from astropy import units as u
  
  # define units
  u_Length   = 1* u.kpc
  u_Mass     = 10**10 * u.M_sun
  u_Velocity = 1* u.km/u.s
  u_Time     = u_Length/u_Velocity 
  toMsol     = u_Mass.to(u.M_sun).value
  
  n = 10
  outputfile = "snapshot.hdf5"
  
  x = np.random.random(n)             # <<<  x in kpc
  y = np.random.random(n)             # <<<  y in kpc
  z = np.random.random(n)             # <<<  z in kpc
  pos = np.transpose(np.array([x, y, z]))
  
  vx = np.random.random(n)            # <<<  vx in km/s
  vy = np.random.random(n)            # <<<  vy in km/s
  vz = np.random.random(n)            # <<<  vz in km/s
  vel = np.transpose(np.array([vx, vy, vz]))  
  
  mass = np.ones(n)                   # <<<  mass in 10^10 solar mass
  
  # create the pNbody object
  nb = Nbody(status='new',pos=pos,vel=vel,mass=mass,ftype='swift')
  # set all particles to stellar particles
  nb.set_tpe(4)
  
  nb.age                =   np.random.random(n)       # <<< age of the stellar particles in Gyr
  nb.mh                 =   np.random.random(n)       # <<< metallicity  [M/H]
  
  # add units
  nb.UnitLength_in_cm         = u_Length.to(u.cm).value
  nb.UnitMass_in_g            = u_Mass.to(u.g).value        
  nb.UnitVelocity_in_cm_per_s = u_Velocity.to(u.cm/u.s).value
  nb.Unit_time_in_cgs         = u_Time.to(u.s).value
  
  nb.hubblefactorcorrection      = False
  nb.comovingtoproperconversion  = False
  nb.atime                       = 1
  
  nb.rename(outputfile)
  nb.write() 
    




If not provided, an estimation for the size of stellar particles can be obtained using the following command
(see :ref:`mockimgs_scripts-label` for mode details)::

  mockimgs_sb_addfields --do_not_compute_ages --do_not_compute_magnitudes snapshot.hdf5 -o snapshot.hdf5 

Note that the input can be also the output (here ``snapshot.hdf5``).
The size may be retrieved via the variable ``nb.rsp`` once ``snapshot.hdf5`` is opened with ``pNbody``.

.. warning::

  The following command may sometimes fail if two stellar particles are too close to each other, which should normally
  not happen. However, if this is the case, the workaround is to "shake" very slightly the particles before saving
  the hdf5 file. For example, in the previous script, add::
    
    x = x + np.random.uniform(-0.5,0.5,nb.nbody)*eps
    y = y + np.random.uniform(-0.5,0.5,nb.nbody)*eps
    z = z + np.random.uniform(-0.5,0.5,nb.nbody)*eps
  
  where `eps` must take a value much smaller than the typical spacial resolution of the system.

Check the model
===============

To check that the created model contains all the necessary fields, you can print the content of the fields 
related to the stellar particles using the  ``sw_getDescriptionFields`` command::

  sw_getDescriptionFields snapshot.hdf5 

which should return::
  
  PartType4

  Coordinates
  Masses
  ParticleIDs
  SmoothingLengths
  StellarAge
  StellarMetallicity
  Velocities

Another way of checking is to run::

  gpy snapshot.hdf5
  
and successively ask for the existence of different arrays::
  
  nb.has_array("pos")
  nb.has_array("vel")
  nb.has_array("mass")
  nb.has_array("num")  
  nb.has_array("rsp")
  nb.has_array("age")
  nb.has_array("mh")
    
  
the output should be::
  
  In [1]: nb.has_array("pos")
  Out[1]: True
  
  In [2]: nb.has_array("vel")
  Out[2]: True
  
  In [3]: nb.has_array("mass")
  Out[3]: True
  
  In [4]: nb.has_array("num")
  Out[4]: True
  
  In [5]: nb.has_array("rsp")
  Out[5]: True
  
  In [6]: nb.has_array("age")
  Out[6]: True
  
  In [7]: nb.has_array("mh")
  Out[7]: True
  

  
  

``mockimgs_sb_addfields``
.........................

.. program-output:: mockimgs_sb_addfields -h