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 Part 1: Operating System Questions

At the terminal/command prompt, Java is not compiling, saying it "cannot resolve symbol StdIn" (or StdDraw, etc) even though I used the introcs installer. Why is this? To use the standard libraries, someone has to tell the Java compiler and runtime where to find them. The introcs installer should have configured DrJava to do this automatically when compiling or running in DrJava. But when using the Command Prompt/Terminal, use java-introcs and javac-introcs in place of java and javac. These are modified commands that tell Java where to find the standard libraries. With the older, more intrusive installers, it may be the case that java/javac were automatically reconfigured to do this.

I used the installer and neither of these commands worked, why is that? Check with a lab TA or preceptor first (you may simply be running things in the wrong directory). If you moved any of the files that the installer added, the environment variables and DrJava preferences set up by the installer will no longer be accurate. If you think you did this, you may reinstall. If this was not the issue, it may be an installer bug which you should diagnose through the lab TAs and on Piazza.

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.

I am a power user and installed everything myself. Where do I get stdlib.jar? Here on the booksite (scroll to top of page) and make sure you add it to both the classpath of DrJava and your CLASSPATH environment variable.

I do not have my own computer, so I am using a university computer cluster machine, and cannot install there. How do I access the COS126 standard libraries? It may work automatically but if not, directly download copies of Stdin.java, StdOut.java, StdAudio.java and StdDraw.java into your assignment work directory (the same folder as NBody.java) from the booksite (scroll to top of page). Compile them and then you can use them. (You'll have to copy them into future weeks' assignment folders too.) The lab TAs can help you with this.

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 Prompt/Terminal.

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.

Frequently Asked Questions Part 2: Programming Questions

What preparation do I need before beginning this assignment? Read Sections 1.4 and 1.5 of the textbook. You also should get familiar with using the command-line.

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!

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.

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

May I write my own functions (static methods)? It's not required and not expected for this assignment, but it is permitted if you want. Here are two possibile ways to modularize your program.

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

May I write my own objects/classes? There are organizational advantages to using objects, which we'll explore later. However, for this assignment, for full credit, you must adhere to the parallel array data structure specified in the assignment. Many well-used programming languages do not allow data abstraction and it is good to be able to write readable, well-organized code subject to that limitation.

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.

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

• Review the input format specified on the assignment page. Here are two sample input files, planets.txt and 3body.txt that conform to this format. Your program should be able to read either one of these files from standard input.

  Read in the data file 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. Does it work? Try running

  % java-introcs NBody 0 1 < planets.txt 
  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
  

  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?

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

• Avoid hard-wired, unexplained constants.
  • For example, in a program about generating and counting digits where 9 is used in several places, define a variable MAX_DIGIT = 9; and change all relevant 9s to MAX_DIGIT.

  • An all-capital name is often understood by readers to indicate a constant that does not depend on the input to the program.

  • This makes your program easy to understand and easy to modify if you want to use a different base (say, binary instead of decimal).

  • 0, 1 or 2 used in an obvious manner is OK, but anything else should either be named, or have a comment explaining where it comes from. This is especially true if the same constant is used more than once.

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.