Programming Assignment Checklist: N-Body Simulation

Goals

Learn about a scientific computing application.
Learn about reading input using the StdIn library and printing formatted output using the StdOut library.
Learn about plotting graphics using the StdDraw library.
Learn about using the command line to redirect standard input to read from a file.
Use arrays.

Frequently Asked Questions

What preparation do I need before beginning this assignment? Read Sections 1.4 and 1.5 of the textbook.

Any advice on completing this programming assignment? Warning: this program is more involved than the previous ones, and you should budget more time accordingly. The best advice we can give you is to carefully test, debug, and re-test your code as you write it. Do not attempt to write the whole program at once—if you do, then you will have no idea where to find the error if the program doesn't work. Proactive testing will save you enormous amounts of time in the long run. Trust us!

What's the easiest way to copy the nbody directory from the COS 126 ftp site to my computer?

Safari. Click the ftp link above to mount the ftp site, then drag-and-drop the nbody folder from the ftp mount to the desired location.
Firefox. Install the DownThemAll! plugin to support downloading all of the links contained in a webpage and use it.
Internet Explorer. Click the ftp link above, then click Page and select Open FTP Site in Windows Explorer from the dropdown menu. Drag the desired folder to your desktop.
Alternative. Click the ftp link above, download the file nbody.zip, and unzip the contents.

DrJava doesn't let me redirect standard input. What should I do instead? Use the command line, as described in lecture. Here are instructions for configuring your system to use the command line: [ Windows · Mac · Linux ]

Do I have to follow the assignment specifications for reading input from the command line and standard input, and writing output to standard drawing and standard output? Yes, or you will lose a substantial number of points.

Where can I find the APIs for StdIn, StdOut, StdDraw, and StdAudio? They are on p. 716–718 of the textbook. They are also available online in this Java cheatsheet.

How do I plot a picture using StdDraw? The command StdDraw.picture(x, y, filename) plots the image in the given filename (either JPEG, GIF, or PNG format) on the canvas, centered on (x, y).

My computer doesn't display the animation smoothly. Is there anything I can do? Use StdDraw.show(30) to turn on the animation mode of standard draw. Call it once at the end of each time step, not after each drawing command.

Can I combine Steps 1, 2, and 3 into one massive loop? No! You must do Step 1 (compute all of the forces) before Steps 2 and 3 (moving and drawing them); otherwise, you will be computing the forces at time t using the positions of some of bodies at time t and others at time t + Δt. Of course, the three steps comprise the body of the main simulation loop.

I draw the bodies, but nothing appears on the screen. Why? Use StdDraw.setXscale() and StdDraw.setYscale() to change the coordinate system to use the physics coordinates instead of the unit box.

My bodies repel each other. Why don't they attract each other? Make sure that you get the sign right when you apply Newton's law of gravitation: (x[j] - x[i]) vs. (x[i] - x[j]). Note that Δx and Δy can be positive or negative so do not use Math.abs(). Do not consider changing the universal gravitational constant G to patch your code!

What is Δx? It's the change in x-coordinate between two bodies (not between a body at time t and at time t + Δt).

How many steps should my program simulate if T is not a multiple of ΔT? Simulate the universe until the t ≥ T. For example, if T is 50,000 and ΔT is 25,000, you should simulate the system for two steps. If T is 60,000 and ΔT is 25,000, you should simulate the system for three steps (until time 75,000).

When I compile NBody.java, it says "cannot resolve symbol StdDraw." Any thoughts? Be sure you have either StdDraw.java or StdDraw.class in the current directory.

How do I submit the extra credit? Submit the extra credit file in the Additional Files dropbox window. Document what you did in your readme.txt file.

How do I change my directory to the H: drive from the Windows Command Prompt? Type H: at the command prompt. Then cd to the appropriate directory.

Everything works until I add StdAudio.play("2001.mid") then everything hangs. What do I do? On some machines there may be a race problem between sending things to the drawing window and sending things to the sound card. Try moving StdAudio.play() to a place in the program after your initial calls to StdDraw.setXscale() and StdDraw.setYscale().

I can play MP3s using my favorite media player, but no sound plays in Java. What could be wrong? If you are running Windows XP, be sure that the audio stream that Java uses is not muted via Start -> Programs -> Accessories -> Multimedia -> Volume Control -> Wave Out.

Testing and Debugging

Inserting print statements is a good way to trace what your program is doing. Our input files were created with the following StdOut.printf() statements.

StdOut.printf("%d\n", N);
StdOut.printf("%.2e\n", R);
for (int i = 0; i < N; i++) {
    StdOut.printf("%11.4e %11.4e %11.4e %11.4e %11.4e %12s\n",
                   rx[i], ry[i], vx[i], vy[i], mass[i], image[i]);
}

Here are our results for a few sample inputs.

% java NBody 0.0 25000.0 < planets.txt           // zero steps
5
2.50e+11
 1.4960e+11  0.0000e+00  0.0000e+00  2.9800e+04  5.9740e+24    earth.gif
 2.2790e+11  0.0000e+00  0.0000e+00  2.4100e+04  6.4190e+23     mars.gif
 5.7900e+10  0.0000e+00  0.0000e+00  4.7900e+04  3.3020e+23  mercury.gif
 0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  1.9890e+30      sun.gif
 1.0820e+11  0.0000e+00  0.0000e+00  3.5000e+04  4.8690e+24    venus.gif

% java NBody 25000.0 25000.0 < planets.txt       // one step
5
2.50e+11
 1.4960e+11  7.4500e+08 -1.4820e+02  2.9800e+04  5.9740e+24    earth.gif
 2.2790e+11  6.0250e+08 -6.3860e+01  2.4100e+04  6.4190e+23     mars.gif
 5.7875e+10  1.1975e+09 -9.8933e+02  4.7900e+04  3.3020e+23  mercury.gif
 3.3087e+01  0.0000e+00  1.3235e-03  0.0000e+00  1.9890e+30      sun.gif
 1.0819e+11  8.7500e+08 -2.8329e+02  3.5000e+04  4.8690e+24    venus.gif

% java NBody 50000.0 25000.0 < planets.txt       // two steps
5
2.50e+11
 1.4959e+11  1.4900e+09 -2.9640e+02  2.9799e+04  5.9740e+24    earth.gif
 2.2790e+11  1.2050e+09 -1.2772e+02  2.4100e+04  6.4190e+23     mars.gif
 5.7826e+10  2.3945e+09 -1.9789e+03  4.7880e+04  3.3020e+23  mercury.gif
 9.9262e+01  2.8198e-01  2.6470e-03  1.1279e-05  1.9890e+30      sun.gif
 1.0818e+11  1.7499e+09 -5.6660e+02  3.4998e+04  4.8690e+24    venus.gif

% java NBody 60000.0 25000.0 < planets.txt       // three steps
5
2.50e+11
 1.4958e+11  2.2349e+09 -4.4460e+02  2.9798e+04  5.9740e+24    earth.gif
 2.2789e+11  1.8075e+09 -1.9158e+02  2.4099e+04  6.4190e+23     mars.gif
 5.7752e+10  3.5905e+09 -2.9682e+03  4.7839e+04  3.3020e+23  mercury.gif
 1.9852e+02  1.1280e+00  3.9705e-03  3.3841e-05  1.9890e+30      sun.gif
 1.0816e+11  2.6248e+09 -8.4989e+02  3.4993e+04  4.8690e+24    venus.gif


% java NBody 31557600.0 25000.0 < planets.txt    // one year
5
2.50e+11
 1.4959e+11 -1.6531e+09  3.2949e+02  2.9798e+04  5.9740e+24    earth.gif
-2.2153e+11 -4.9263e+10  5.1805e+03 -2.3640e+04  6.4190e+23     mars.gif
 3.4771e+10  4.5752e+10 -3.8269e+04  2.9415e+04  3.3020e+23  mercury.gif
 5.9426e+05  6.2357e+06 -5.8569e-02  1.6285e-01  1.9890e+30      sun.gif
-7.3731e+10 -7.9391e+10  2.5433e+04 -2.3973e+04  4.8690e+24    venus.gif

Possible Progress Steps

These are purely suggestions for how you might make progress. You do not have to follow these steps. If you get stumped or frustrated on some portion of the assignment, you should not hesitate to consult a preceptor.

Copy the nbody directory from the COS 126 ftp site to your computer in a new directory nbody. In particular, be sure you have StdDraw.class, StdIn.class, StdOut.class, StdAudio.class, the image files, and the sample input files are in your working directory.
Read Section 1.5. Review DeluxeBouncingBall.java for help with standard drawing. There are many methods in StdDraw, but for this assignment, you won't need more than the ones that are used in DeluxeBouncingBall.java. Review Students.java for help with standard input. Compile and execute both of these programs to make sure your system is configured properly. For DeluxeBouncingBall.java you will also need the sound files laser.wav and pop.wav. For Students.java you can use the file students.txt as input.
Read in the data file planets.txt from standard input, storing the number of bodies in an integer variable N, the radius of the universe in a floating-point variable R, and the remaining information in six parallel arrays: Let rx[i], ry[i], vx[i], vy[i], and mass[i] be floating-point numbers which store the current position (x- and y-coordinates), velocity (x- and y-components), and mass of body i; let image[i] be a string which represents the filename of the image used to display body i.
Print the information back out in the specified format. Right now, you are printing the information to aid in debugging; later, you will use this code to print out the state of the universe at the end of the simulation.
Do not even think of continuing until you have tested and debugged this part.
Step 3. Plot the background starfield.jpg image. Note that StdDraw.picture(x, y, file) centers the picture on (x, y). Test that it works. Now, write a loop to plot the N bodies. If all goes correctly, you should see the four stationary planets and the sun.

Test it on another data file.
Step 2. Write a loop to calculate the new velocity and position for each body. (This code goes before the plotting code you wrote in Step 3.) Since we haven't yet incorporated gravity, assume the acceleration acting on each body is zero. Add an outer loop to repeat Steps 2 and 3 until the time reaches T. You should now see the four planets moving off the screen in a straight line, with constant velocity. Test it on another data file.
Step 1. Now, calculate the net force acting on each body. (This code goes before the stuff you wrote in Step 2.) You will need two additional arrays fx[i] and fy[i] to store the net force acting on body i. First, initialize all the net forces to 0.0. Then write two nested for loops to calculate the net force exerted by body j on body i. Add these values to fx[i] and fy[i], but skip the case when i equals j. Once you have these values computed, go back to Step 2 and use them to compute the acceleration (instead of assuming it is zero). Test your program on several data files.
Finally, add the music. At the beginning of the simulation (outside of any loop), use StdAudio.play("2001.mid") to play the 2001 theme. Test it.

Enrichment

What is the music in 2001.mid? It's the fanfare to Also sprach Zarathustra by Richard Strauss. It was popularized as the key musical motif in Stanley Kubrick's 1968 film 2001: A Space Odyssey.
I'm a physicist. Why should I use the leapfrog method instead of the formula I derived in high school? In other words, why does the position update formula use the velocity at the updated time step rather than the previous one? Why not use the 1/2 a t^2 formula? The leapfrog method is more stable for integrating Hamiltonian systems than conventional numerical methods like Euler's method or Runge-Kutta. The leapfrog method is symplectic, which means it preserves properties specific to Hamiltonian systems (conservation of linear and angular momentum, time-reversibility, and conservation of energy of the discrete Hamiltonian). In contrast, ordinary numerical methods become dissipative and exhibit qualitatively different long-term behavior. For example, the earth would slowly spiral into (or away from) the sun. For these reasons, symplectic methods are extremely popular for N-body calculations in practice. You asked!
Here's a more complete explanation of how you should interpret the variables. The classic Euler method updates the position uses the velocity at time t instead of using the updated velocity at time t + Δt. A better idea is to use the velocity at the midpoint t + Δt / 2. The leapfrog method does this in a clever way. It maintains the position and velocity one-half time step out of phase: At the beginning of an iteration, (r_x, r_y) represents the position at time t and (v_x, v_y) represents the velocity at time t - Δt / 2. Interpreting the position and velocity in this way, the updated position (r_x + Δt v_x, r_y + Δt v_y). uses the velocity at time t + Δt / 2. Almost magically, the only special care needed to deal with the half time-steps is to initialize the system's velocity at time t = -Δt / 2 (instead of t = 0.0), and you can assume that we have already done this for you. Note also that the acceleration is computed at time t so that when we update the velocity, we are using the acceleration at the midpoint of the interval under consideration.
Here are some interesting two-body systems, perhaps relevant to the extra credit. Here is beautiful 21-body system in a figure-8 - reproducing this one will definitely earn you extra credit.
Here is a website that generates music using an N-body simulator!
Here's a wealth of information on N-body simulation [pdf].
N-body simulations play a crucial role in our understanding of the universe. Astrophysicists use it to study stellar dynamics at the galactic center, stellar dynamics in a globular cluster, colliding galaxies, and the formation of the structure of the Universe. The strongest evidence we have for the belief that there is a black hole in the center of the Milky Way comes from very accurate N-body simulations. Many of the problems that astrophysicists want to solve have millions or billions of particles. More sophisticated computational techniques are needed.
The same methods are also widely used in molecular dynamics, except that the heavenly bodies are replaced by atoms, gravity is replaced by some other force, and the leapfrog method is called Verlet's method. With van der Waals forces, the interaction energy may decay as 1/R^6 instead of an inverse square law. Occasionally, 3-way interactions must be taken into account, e.g., to account for the covalent bonds in diamond crystals.