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!

How do I get the files needed to complete the assignment? There are two approaches: download a single compressed .zip file and then decompress it on your computer to the right location; or, grab the folders themselves using the File Transfer Protocol (FTP). If you've never used either of these tools, follow one of the methods below. (You'll need to do a similar process each week from now on.)

1. Zip files on Windows. Download the nbody.zip file listed on the assignment page. In Explorer, view the directory that you downloaded the zip file to. Right-click on the zip file and "Extract All..." which will ask you to choose the extraction location. (Don't double click the zip file, since it browses the contents but does not extract them.) A new subdirectory nbody will be created at the extraction location, with the files inside of it.

2. Zip files on OS X. Download the nbody.zip file listed on the assignment page. In Finder, view the directory that you downloaded the zip file to. Some browsers will automatically have extracted the zip file and created an nbody folder for you with the extracted files inside of it. (Some also will delete the original .zip file.) Otherwise if you see nbody.zip but no nbody folder, double-click on the .zip file to extract its contents. Finally, drag or copy-and-paste the nbody folder to your desired location.

3. FTP on Windows. Click Run (Windows-R) and enter explorer ftp://ftp.cs.princeton.edu/pub/cs126/ which will let you browse our FTP site. Right-click on nbody (the folder for this week) and select "Copy to Folder...". A new subdirectory nbody will be created at the target, with the files inside of it.

4. FTP on OS X. Either (a) in Safari, navigate to ftp://ftp.cs.princeton.edu/pub/cs126/nbody or (b) in Finder, press Command-K and enter this as the Server Address. If prompted for a name and password, choose to Connect as a Guest. This lets you browse our FTP site. Drag-and-drop the nbody folder to your desired location on your computer to make a copy of it.

5. Linux. Execute wget -r ftp://ftp.cs.princeton.edu/pub/cs126/nbody at the command line.

• You can also use any other zip extractor or FTP client if you have a preferred one.

• After following one of these sets of instructions, check to see that the files including earth.gif and planets.txt are indeed there.
  

• Note: all the instructions below assume your NBody.java file will be in the same directory as planets.txt and all the other files.

DrJava doesn't let me redirect standard input. What should I do instead? DrJava does not support redirection or piping; instead, you must use the command line.

Do I have to follow the assignment specifications for reading input from command-line arguments 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.

I did not use the installer. How do I install stdlib.jar? We highly recommend using the Windows installer or the Mac OS X installer. If you choose not to, you can refer to the standard libraries, which includes instructions on downloading and using them.

When I compile NBody.java, it says "cannot resolve symbol StdDraw." Any thoughts? You do not have our standard libraries properly installed. See the question above. If you tried to install stdlib.jar and still get this message, please consult a staff member (preceptor or lab TA) to fix this.

I do not have my own computer, so I am using a university computer cluster machine. How do I access the COS126 standard libraries? Refer to the instructions on using the standard libraries. The lab TAs can help you with this.

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).

Do I need to write my own functions (static methods)? It's not required. If you want to modularize your program, here are two possibilites.

// return Euclidean distance between (x1, y1) and (x2, y2)
private static double distance(double x1, double y1, double x2, double y2)

// return the magnitude of the gravitational force between two bodies
// of mass m1 and m2 that are distance r apart
private static double force(double m1, double m2, double r)

My computer doesn't display the animation smoothly. Is there anything I can do? Use void StdDraw.show(int msec) to use the animation mode of standard draw and pause for msec milliseconds. Call it once at the end of each time step (frame), not after each drawing command. A typical delay time would be between 10 and 50 milliseconds, but use whatever looks good on your computer.

Can I combine Steps 1, 2, and 3 for each body? You must separate Step 1 (compute all of the forces) from Step 2 (moving 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.

My animation looks fine, but the numbers are off a little bit when I follow the tests below, what could be causing this? Read the question just above — this is the most common cause. Double-check that you followed the assignment instructions exactly in the order of how positions and velocities are updated.

I draw the bodies, but they do not appear 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. Since we want it centered on the origin with a "square radius" of R, the minimum of each axis should be –R and the maximum should be +R.

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 as long as t < T. For example, if T is 50,000 and ΔT is 25,000, you should simulate the system for two steps (t = 0 and 25,000). If T is 60,000 and ΔT is 25,000, you should simulate the system for three steps (t = 0, 25,000 and 50,000).

Does my program need to plot anything if T equals 0? No, you do not need to plot anything.

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.

Everything works until I add StdAudio.play("2001.mid") then everything hangs. What do I do? On some machines there may be a timing conflict (known in computer science as a "race") between 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.

When I add the music, I get this weird error. How can I fix it? java.util.prefs.WindowsPreferences (init). Could not open/create prefs root node Software\JavaSoft\Prefs at root 0x80000002. Windows RegCreateKeyEx(...) returned error code 5. Try running your command prompt as administrator (right-click on the shortcut icon). Then use it as before and try again. If this doesn't help, contact the course staff.

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",
                   px[i], py[i], vx[i], vy[i], mass[i], image[i]);
}

% 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.

• Get the files, by following the instructions above.

• 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 those 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 px[i], py[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.

  Next to your code for reading from standard input, parse the two command-line inputs and store their values for later.

  Do not even think of continuing until you have tested and debugged this part.

• Between the code for your input and printed output, create the "time" loop. Its body will perform a single time step of the simulation, and increase the current time counter by Δt. Steps 3, 2 and 1 (below) will go inside this time loop. Use the information you parsed from the command-line inputs to help set up this loop. Check that your code still compiles and runs.

• Before entering the time loop, set the scale of the drawing window using the radius of the universe. The center of the drawing window should correspond to the center of the universe which is at the origin (0, 0).

• Step 3. Get a basic stationary picture working. Use StdDraw.picture(x, y, file) to center a picture of the background starfield on the origin. Then, write a loop to display the N bodies. If all goes correctly, you should see the four stationary planets and the sun. Test that it works. Finally, to stop the image from flickering, use the StdDraw.show(msec) method to turn on "animation mode." Test this as well; furthermore, see if it works on a universe other than planets.txt.
  

• 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. 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.

Reviewing your Program

• Review your header, comments, and style, keeping in mind the guidelines from previous weeks.

• Do you list the dependencies for each program in the header?

• Are numerical constants used in a way that is easy for the reader to understand?

• Upon submission, does checkstyle give you any important warnings?

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, (p_x, p_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 (p_x + Δt v_x, p_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.

• Are there analytic solutions known for N-body systems? Yes, Sundman and Wang developed global analytical solutions using convergent power series. However, the series converge so slowly that they are not useful in practice. See The Solution of the N-body Problem for a brief history of the problem.

• 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.