Units in gravitational $N$ body simulations I am trying to write a code in Python to simulate $N$ bodies interacting through gravity. In particular I am trying to see whether a system of particles with random initial positions and zero velocity will fall into a viral equilibrium.
I understand that in $N$ body simulations, it is advisable to set $G=1$ and $M=1$. However, having done this what then are my units of time or length or energy?
 A: Kepler's Third Law is
$$\frac{T^2}{r^3} = \frac{4 \pi^2}{GM}$$
or equivalently
$$
\frac G{4\pi^2} = \frac{r^3}{T^2M}
$$
If you're modeling our solar system you can set 
\begin{align}
G &= 4\pi^2\frac{\text{AU}^3}{\text{year}^2 M_\odot}
\end{align}
where "AU" is the astronomical unit and $M_\odot$ is the mass of the sun. Then if you put an object with mass 1 (the sun) at the origin and a negligibly-massive object in a circular orbit at radius 1 (a.u., where earth orbits) you should find a period of 1 (year).
Here's a code snippet from a similar program of mine:
G = 4 * np.pi**2        ## in solar masses, AU, years
msun, au, year = np.ones(3)
km  = 1.0 / 1.5e8       ## one kilometer
m   = km  / 1e3
kg  = 1.0 / 2e30        ## one kilogram
sec = 1e-7 / pi         ## one second
joule = kg * (m / sec)**2

Now if I have a quantity whose dimension is energy, I can divide by joule to get a value in those units, even though in my strange "natural" units energies will have values near unity.
A: Most of the time, scientific computer code is written in such a way that variables have no "knowledge" of the units they are intended to represent. (Of course, you could be arbitrarily sophisticated in the way you write your program, e.g. by defining classes that keep track of dimensionality and used units, and then use these classes to define your variables. But that is the exception, in my experience.)
Assume, for generality, that you have bodies of different masses. You may choose one of them, called, say, $M_0 \ne 0$, to be the scale you use to measure masses. So, the mass of object $1$ may now be measured in multiples of $M_0$, e.g. you could have $M_1 = 42 M_0$, and likewise $M_2 = 0.01 M_0$. It is then common to set $M_0=1$, both on paper as well as in the program, and consequently drop the term $M_0$ altogether. Similarly, you can set $G=1$. This is often done, but in what follows, keep the initial situation with $M_0$ still present in mind.
So, $R$ being the distance between $M_1$ and $M_2$, you may write for the potential energy $V$, either on paper or in your code,
$$
V = - \frac{42 \cdot 0.01}{R}
$$
Now, $R$ may as well be measured as multiples of some characteristic length $R_0$ which you set to unity, so that all that remains for $V$ is a real number.
Your question, if I understood it correctly, is now how one can make sense of a calculated result of, say, $V = - 0.42$ (as in the case where, by coincidence or design, $R=1$.)
To do this, it may be helpful to review what has been done to reach this point: By convention, we agreed to measure masses as multiples of $M_0$, gravitational coupling strengths as multiples of $G$, and distances as multiples of $R_0$. However, let's not set $M_0=G=R_0=1$ this time, and consequently let's not drop those constants from our expressions.
So, instead of
$$
V = - G \frac{M_1 M_2}{R}
$$
one can write without loss of generality (see e.g. the definitions of $M_1$ and $M_2$ above)
$$
\left(v \cdot V_0\right) = - \left(g \cdot G\right) \frac{ \left(m_1 \cdot M_0\right) \left(m_2 \cdot M_0\right)}{\left(r \cdot R_0\right)}
$$
where all of $g=1$, $m_1 = 42$, $m_2 = 0.01$, $v$, and $r$ are dimensionless real numbers, and are the quantities you probably encounter in your program. The potential energy is now measured as a multiple of some characteristic energy $V_0$, i.e. $V=v \cdot V_0$.
So, what your computer program may calculate is actually the dimensionless quantity $v=-0.42$, which measures the multiples of $V_0$, so that the potential energy, with the correct dimensions, is $V=-0.42 V_0$. What remains is to understand what $v=-1$, i.e. $V=-V_0$, means in terms of, say, SI units. You would achieve $v=-1$ for example for $m_1=m_2=g=r=1$, so that
$$
-V_0 = - G \frac{M_0^2}{R_0}
$$
Depending on the values chosen for $M_0$ and $R_0$ (and $G$, if that is regarded a "choice"), $V_0$ is now fixed and known as well.
If you wanted to produce a graph of the dependence of the potential energy on the spatial separation of two fixed masses, you would the probably label the x-axis either $r$ or $R/R_0$, and the y-axis either $v$ or $V / \left(G \frac{M_0^2}{R_0}\right)$.
I hope this helps. Please feel free to ask for clarifications.
