Velocity distribution in Plummer's models and others mass distributions The Plummer's sphere is an model for the mass density in a globular cluster of stars. For an $N$-body simulation I have initialized the position of $N$ masses with a Monte-Carlo technique but cannot find a way of initializing the velocity initial conditions.
Is there a simple function that given a position in a Plummers sphere assigns a velocity to a given mass? Lots of sites list the velocity for a circular orbit but is this a good approximation to a globular cluster and how should it be treated off the $x$-$y$ plane?
 A: The Plummer model has a potential of the form
$$
\Phi(r)=-\frac{1}{\sqrt{r^2+1}}
$$
(obviously ignoring all constants). Equating the above with the kinetic energy, you get
$$
\frac12v_e^2+\Phi(r)=0\to v_e(r)=\sqrt{2}\left(r+1\right)^{-1/4}
$$
This velocity is the maximum velocity you can have at a radius $r$, so we must have that $0\leq v\leq v_{e}$. In order to get the velocity, we need to know its distribution function; this can be found in this book. It reduces down to
$$
g(v)dv\propto\left(-E\right)^{7/2}v^2dv
$$
where $E=-v_e^2+\frac12v^2$. If we let $q=v/v_e$, then the above becomes
$$
g(q)=\left(1-q\right)^{7/2}q^2
$$
with $0\leq q\leq1$. It then becomes a matter of using either rejection sampling (which I generally do not recommend) or the Newton root-finding method (which I usually do recommend) to find a fit. Then you just randomize the angles $\theta,\,\phi$ and get your x,y,z velocities:
vel = x*sqrt(2.0)*(1.0+r(i)*r(i))**(-0.25)
theta = acos(rand(-1.0,1.0))
phi = rand(0.0, 2.0*pi)
vx(i) = vel*sin(theta)*cos(phi)
vy(i) = vel*sin(theta)*sin(phi)
vz(i) = vel*cos(theta)

You can find more information in this great pdf book.
