Distribution of orbital velocities in a disk galaxy for N-body simulation? I'd like to write an N-body simulation in which I collide two disk galaxies. To give you an idea of the accuracy I'm trying to achieve, I'm aiming to make this my screensaver at 30fps on my work desktop. Think Euler's method, large epsilons in denominators to prevent blow up, all that lazy stuff. 
One place I don't want to be lazy, however, is in the initial construction of the galaxies I'm trying to collide. I want to design a galaxy that will at least sort of hold together by its own gravity until it smash it into another one. 
How might I go about selecting the orbital velocities for each of the bodies composing the galaxy? 
 A: You could probably expect that $\langle v_r\rangle=\langle v_z\rangle=0$--that is, most of your velocities will be in the $\hat\phi$ direction (assuming you are using cylindrical coordinates).1
For stable orbits, the kinetic and potential energies must be equal, so you should end up seeing that
$$
v(r)=\sqrt{-2\Phi(r)}
$$
where $\Phi$ is the gravitational potential energy (and there are a few choices). You also can add Normally-distributed values to the velocity components (as suggested here) to give a little perturbations to the orbits. It is pretty straight-forward to convert to another coordinate system from here, if need be.
I also discuss in this related question an algorithm to pick the velocities assuming a Plummer model (discussed in the first link in the post). I expect a similar algorithm could be developed for the alternative potentials, but haven't worked the math of it.

Note also that, while the system of equations for the $N$-body problem is pretty straight-forward,
\begin{eqnarray}
\mathbf v&=\frac{\mathrm d\mathbf x}{\mathrm d t}\\
\mathbf F&=m\frac{\mathrm d\mathbf v}{\mathrm dt}
\end{eqnarray}
writing it in code is actually quite hard for large $N$. The force $\mathbf F$ is computed between each pair of objects in the $N$-body system:
$$
\mathbf F_i=-\sum_{j\neq i}G\frac{\hat{\mathbf x}_{ij}}{\left|\mathbf x_{ij}\right|^2}
$$
where $\mathbf x_{ij}=\mathbf x_i-\mathbf x_j$ and $\hat{\mathbf{x}}_{ij}$ is the direction of the force between bodies $i$ and $j$. This means you automatically have at least an $\mathcal O(N^2)$ problem (at least in time, it's $\mathcal O(N)$ in memory).2
A real galaxy has $N\sim10^{11}$ stars, which is probably not possible for any modern computer to handle.3 I believe that $N\sim10^6$ is a good round number for modelling galaxies, but these can still take hours on computer clusters (depending on what the simulation does) and I'll go out on a limb and guess that you don't have a computer cluster available to you, so you may want to try $N\sim10^3$ or $N\sim10^4$ instead--but don't be surprised if this still takes a very long time!

1. It is probably easiest to use $r,\phi,z$ coordinates for generating the velocities, then transforming them back to Cartesian (assuming you are working in Cartesian).
2. There are tricky algorithms that can reduce this to $\mathcal O\left(n\log n\right)$, but this may be a bit much for your purposes.
3. Though there was the one trillion body simulation, they used the whole K-computer for something like a month straight (if I recall the news articles correctly). I doubt you have that type of resources available.
