for (i = 0; i < nSt; i++) {
j = (int) (x[i] * scale);
k = (int) (y[i] * scale);
printf("%d %d^", j, k);
}
printf("\n");
The nBody program.
The program is currently configured to run on a Linux system, and these instructions assume that that is the case. This is a very CPU intensive application, so the faster your system, the better your results. It was developed on a 64 bit system, but will run on a 32 bit system. You'll need to have the GNU C development suite installed, including gcc and make. The debugger, gdb, is also recommended, as, of course, is a familiarity with C and programming in general.
These programs are designed with the idea that you, the user, will edit the code to enhance/tweak/modify it to make it do your bidding. I have tried to comment it so that what it's doing should be clear.
What this program does is take an initial arrangement of stars, specified by newRun.c, and allow them to gravitationally interact, using standard Newtonian gravity, for the number of iterations specified in orbit.c.
The simple story (as coded in smallCluster.c). For each iteration, each star looks at every other star and calculates the distance between them. Then the gravitational force between the stars is calculated by, you guessed it:
Fg = G m1 m2 / r2.
Of course, we need to allow for the vector components, so the real equation is, in three dimensions:
Fgx = G m1 m2 X / r3.One then sums these forces into the total force, and then calculates the acceleration for a given star, again calculating for all three dimensions:
ax = F / m(i).The change in velocity is:
vx = vx + (ax dt / 2)Finally, the position changes by:
x = xo + vx.over the three dimensions.
Simple, huh? With this the planets move, with considerable accuracy, as Feynman shows in volume 1 of his physics lectures.
Stars aren't so simple. The problem is when two stars happen to get very close to one another, the force becomes quite large. Then the velocity change is greatly effected, to, and this passes on to the position. Bingo. Suddenly these stars are quite far apart, and they have high velocity. In reality, as the stars separate, their mutual gravity slows them down, but not here. We only calculate gravity each dt, not continuously, as nature does. We've added fictitious energy into the stem. Something must be done.
It should also be noted that stars very far from the star you are currently moving have only a minimal impact on your star. This is inefficient.
Orbit deals with these issues by only looking at stars that are within a cube that is accDst * 2 in size. These stars are given the Newtonian treatment, all other stars are lumped into a single mass at the cluster's center, and then treated as a single star. It also checks the distance between the two stars in the cube, and if it's smaller than minRadius, it puts them into a group of close stars (*group). This small group is the processed with the Newtonian treatment, except that dt is modified so the stars would only move at most about 10 degrees in their orbit per dt.
How it works:
NewRun.c creates as many stars as you desire (nSt) and gives them an initial position and velocity. In practice, users back up newRun.c to another file, say originalNewRun.c and copy a version of newRun.c that configures the simulation to their current project.
Orbit.c moves the stars. Each iteration begins by clearing the acceleration and hash variables, updating each star's position, and calculating the amount of mass is between the center of the star and a series of radii from it (nRadii). This latter quantity is used for stars greater than accDis from the star whose motion is being calculated (current star).
Also calculated at this time is a series of hashes that place a star in the three coordinates, x, y, and z. Later, these hashes will be searched for stars close to current star. Not having to check all of the stars for each current stars is a helpful speed up.
The iteration continues by searching the hashes for stars near the current star and saving them into an array, near.
For all stars in near, the mutual distances are calculated and stored in the two dimensional r array. If the stars are closer than minRadius, they are placed in the group array. After all of the stars not in the group array are processed, the group array as a whole is processed, using only the stars in the group, and making the time interval as small as needs be to prevent adding the factious energy mentioned above.
When the contributions from all the nearby stars have been added up, all other star are treated as if they were at the center of the cluster. These are then added to the force as a whole, and the star is ready to be moved.
Stars near the cluster's center are all treated as if they are close together.
Every iteration/iterCt times, orbit calls ckMomentum. This calculates the virial theorem for the cluster, which is a good measure of how it's doing. If the (potential energy) / (kinetic energy) is not 2, something's amiss. There are also quite a few diagnostic routines that can be enable by commenting out one or the other and recompiling.
All of the units used by nBody are MKS. I know, one can rationalize the units, (e.g. G == 1) and avoid a multiplication, but this can be done after everything else is shown to be correct.
Unlike most programs you find on the internet, nBody expects that you're going to edit the code itself to change parameters. If you mess up, you can always download it again and start over! So a short description of nBody's components is in order.
nBody.c: The program that moves the stars and generates the results.
nBody.h: Variables common to two or more modules, system includes.
nBodyDisplay.py: A python program that displays star movement projected onto the X-Y plane. You'll need to have output formatted by orbit.c with this snippet for it to work:
// Convert the positions into accurate integers for the python visualizer.
for (i = 0; i < nSt; i++) {
j = (int) (x[i] * scale);
k = (int) (y[i] * scale);
printf("%d %d^", j, k);
}
printf("\n");
N Body simulation of 1000 stars after 1,000,000 interations.
plotSmallCluster.py: This plots the results of smallCluster.c much as nBodyDisplay.py plots the results of nBody.
Makefile: Compile everything, and make an executable called nBody. Just type "make" to do this.
newRun.c: This sets up the stars, and assigns each a position, velocity and mass.
readme: This document.
output directory: Every N iterations, write out the status of the program and eliminate any systematic drift (all velocities should sum to zero). This is done by updateRun, and the results are stored here.
runBigNBody: Runs the nb simulation. Errors are put in err, data in nbOut. Snapshots of the X and Y coordinates of the stars at given intervals are written to the output directory.
smallCluster.c: Apply the simple method to all stars. Don't check for close encounters. O n^2, but fast enought for less than 100 stars. It's currently set up for 16 stars. You can compile it with:
gcc -o smallC -Wall -g smallCluster.c -lm
This will create a debuggable executable called smallC.
Use this program to accurately simulate small numbers (<100) of stars.
SmallCluster simulation of 16 stars.
Currently the program only runs for a fixed number of iterations. Enhance it so that if the user starts it with a data file name, it reads in that file and restarts the program.