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 \(\rm{kpc}\) (nb.pos)

• velocities for each stellar particle, in \(\rm{km/s}\) (nb.vel)

• mass for each stellar particle, in \(10^{10}\,\rm{M_\odot}\) (nb.mass). The mast must corresponds to the initial SSP mass.

• age for each stellar particle, in \(\rm{Gyr}\) (nb.age) An old stellar particle will have a typical age of \(12\,\rm{Gyr}\) and a recently formed particles of \(2\,\rm{Gyr}\)

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

  \[[M/H] = \log_{10} \frac{M_{\rm{M}}/M_{\rm{H}}}{ \left( M_{\rm{M}}/M_{\rm{H}}\right)_{\odot} }\]

  where \(M_{\rm{M}}/M_{\rm{H}}\) is the mass fraction of metals and usually denoted by \(Z\) and \(\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 \(\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 (\(h\))!

The following lines show how to get a snapshot of the appropriate format using 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 Scripts 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

usage: mockimgs_sb_addfields [-h] [-o OUTPUTFILENAME] [--nngb NNGB]
                             [--do_not_compute_ages] [--do_not_compute_rsp]
                             [--do_not_compute_magnitudes] [--ftype FTYPE]
                             [FILE ...]

positional arguments:
  FILE                  a file name

options:
  -h, --help            show this help message and exit
  -o OUTPUTFILENAME     Name of the output file
  --nngb NNGB           Number of neighbouring particles to consider to
                        compute RSP(==HSML)
  --do_not_compute_ages
                        do not compute ages
  --do_not_compute_rsp  do not compute rsp
  --do_not_compute_magnitudes
                        do not compute compute_magnitudes
  --ftype FTYPE         type of file