Simulation of orbiting bodies I am writing a computer program to simulate orbiting bodies such as planets and stars.
I wish to have a starting point in which a number of bodies are randomly scattered around a central heavy body. Each body (except the central one) is sent off at the correct orbiting velocity but in a random direction.
I am using an x,y,z coordinate system.
For each orbiting body, I can establish the starting position and orbiting velocity.
I want each orbiting body to start moving in a random direction (but on the right plane) for which I need to calculate the initial x, y and z velocity components. This bit I cannot figure out how to do.
Does anyone have any suggestions?  
 A: Do I correctly restate your question as : "I want to get a vector of random orientation in the plane, but with a given magnitude" ? Then you should take an angle according to a uniform distribution in $[0, 2\pi)$, and then calculate $x$ and $y$ components with the magnitude you want to impose.
A: The task posed is to create the initial velocity as if the body under consideration was alone in the system with the fixed center of mass $M$ and has a circular orbit in a random plane. Of course this initial order will be beautifully destroyed in the first steps of the simulation.
For the randomly generated vector $\vec x_0=r\hat s$, $\|\hat s\|=1$, generate a likewise random, but normalized orthogonal vector $\hat n$. Then the operator $(\hat n\times)$ acts like a complex unit $J$, $J^2=-I$, in the plane orthogonal to $\hat n$. The circular orbit can thus be written as
$$
\vec x(t)=e^{J\omega t}\vec x_0
$$
using the matrix exponential. Now insert into the gravity equation
$$
-ω^2\vec x(t)=\ddot{\vec x}(t)=-GM\frac{\vec x(t)}{\|\vec x(t)\|^3}
\\~\\
\implies ω=\sqrt{\frac{GM}{r^3}}.
$$
This determination now can be translated into the initial velocity
$$
\dot{\vec x}(t)=\omega J\vec x(t)\implies 
\color{blue}{\vec v_0=\sqrt{\frac{GM}{r^3}}(\hat n\times\vec x_0)}=\sqrt{\frac{GM}{r}}(\hat n\times\hat s)
$$
The initial speed can be read off as $\|\vec v_0\|=ωr=\sqrt{\frac{GM}{r}}$.
