litt: S.J. Aarseth, M. Henon and R. Wielen (1974), Astron. and Astrophys. 37, p. 183.
Calculating the coordinates is easiest in STRUCTURAL units; conversion to VIRIAL units will be performed below.
Recipe for scaling to the proper system of units:
Since G = M = 1, if we switch from a coordinate system with length unit r_old to a coordinate system with length unit r_new , the length units simply scale by a factor C = r_new / r_old . Consequently, the coordinate values of physical quantities such as positions should transform inversely to maintain the same coordinate-invariant meaning. Similarly, the square of the velocities should transform inversely proportional to the positions, since GM = 1 (cf. a relation such as v*v = G*M/r ).
To summarize: If
r_unit(new) = C * r_unit(old) ,
then
pos(new) = (1/C) * pos(old)
and
vel(new) = sqrt(C) * vel(old) .
the position coordinates are determined by inverting the cumulative mass-radius relation, with the cumulative mass drawn randomly from [0, mfrac]; cf. Aarseth et al. (1974), eq. (A2).
the velocity coordinates are determined using von Neumann's rejection technique, cf. Aarseth et al. (1974), eq. (A4,5). First we take initial values for x, the ratio of velocity and escape velocity (q in Aarseth et al.), and y, as a trick to enter the body of the while loop.
If the mass fraction and radius fraction are both 1.0 then particles will be sprinkled in all over space. If mfrac < 1.0 or rfrac < 1.0 then each particle is constrained to lie within both the radial and (cumulative) mass bound. For example, if rfrac( mfrac ) > rfrac then rfrac is the limiting factor, but if rfrac( mfrac ) < rfrac then mfrac limits the extent of the Plummer realization.
-i