One other way to generate velocities is to use the velocity moments (as compared to the full distribution functions). For a rather large range of models (spherical, axisymmetric, triaxial, multi-components), it is possible to derive the first two velocity moments along each principal axis (e.g. R, Z for cylindrical coordinates). We can then simply assume any simple functional form which shares the same moments to generate velocities consistently with the assumed dynamical model.

A common practice in this context is to assume a Gaussian distribution with its mean being the first velocity moment (the mean velocity in this direction) and its root mean square being its second (centred) velocity moment, namely the stellar velocity dispersion along that direction.

For each particle, it is then easy to generate the corresponding velocities:
 first get the position of the particle (e.g. R, z, theta)
 then draw randoms numbers using the mean velocities and dispersions along these directions at that specific position. One thing to remember is to reject velocities which are above the local escape velocity ($\sqrt2 \phi).

To illustrate this method, we attach here a python routine which includes modules to make an N body realization of a so-called MGE (Multi-Gaussian Expansion) model. This model is just a representation for a luminosity or mass model using a sum of tri-dimensuonal Gaussians; This simplifies the derivation of the gravitational potential, and allows the analytical derivation of the projected luminosity for any viewing angle. See Monnet et al. 1992, A&A 253, 366 and Emsellem et al. 1994, A&A 285, 723 for details (one application of this formalism to the dynamical modeling of galaxies can be found in Emsellem et al. 1999, MNRAS 303, 495).

So the steps to generate a set of particles with their positions and velocities is rather simple. First set up the parameters for an MGE model of your choice. This is done using the MGE class and specifying a specific ascii file where the parameters are provided (an example of such an ascii file is provided below). Then you just generate the positions/velocities using the function init_nbody.

This would look like (# are comments): # use the mge module

%run mge.py

# reads the ascii file and initialise the model

model = MGE("MGE_N3377.mge")

# make the N body realization

model.init_nbody(nbody=100000, Nquad=100, Rcut=10000, Zcut=10000, kSatoh=1.)

This will provide 100 000 particles with a cut-off in Z and R of 10 kpc and a Satoh parameter (setting the anisotropy of your model) of 1 (meaning here = isotropy). Details can be found in the papers mentioned above, or in the python routine (and otherwise, use emails...).

An example of an MGE ascii file can be found here:

MGE_N3377.mge

Ascii file including the MGE parameters for the model of a galaxy (here: NGC 3377).

and the python module can be found there:

mge.py

This is a python module which allows you to make an N body realization consistent with a Multi-Gaussian Expansion model.