Programming Assignment Checklist: N-Body Simulation

Goals

• Learn about a scientific computing application.
• Learn about reading from standard input using StdIn.
• Learn about redirecting standard input to read from a file.
• Learn about plotting using Turtle graphics.
• Learn about arrays.

Frequently Asked Questions

The submission system says my program exceeds the 80 character limit, but I don't see which line it's on. Can you help? The program WidthChecker.java and helper library CharStdIn.java print out all lines longer than a command line input N. Please note that all tabs are converted to 4 spaces during submission.

My computer it too slow to display the animation smoothly. Is there anything I can do? Here are a few suggestions. First, be sure that you are only calling Turtle.pause once at the end of each time step, instead of after each Turtle.spot command. You might also try increasing the delay parameter in Turtle.pause to prevent hoarding system resources. Or do several updates, say 20, before calling Turtle.pause.

Can I combine Steps 1, 2, and 3 into one massive loop? No! This will simultaneously screw up the physics and make your code harder to understand and debug.

I draw the planets, but nothing appears on the screen. Why? Make sure that the drawing coordinates are between 0 and 512. If not, rescale them so that they are. Also, be sure that you call Turtle.pause at the end of each time step.

How do I scale the bodies so that they stay inside the 512-by-512 Turtle graphics window? The first number in the data file is the radius R of the universe. You can assume all of the x and y coordinates of the planets remain between -R and R. Your job is to compute scaled and translated coordinates that are always between 0 and 512. Compute turtleX and turtleY with ax + b and ay + b for appropriate constants a and b. Note: you should not actually change the x and y coordinates of each planet (this will ruin the physics); instead use the temporary variables turtleX and turtleY to temporarily store the scaled coordinates you intend to plot.

I'm confused about all of the &Delta t / 2 terms. Do I need to worry about them to get the physics right? No! The update formulas for velocity and position already take this into account.

What should I use for the initial velocity in the leapfrog method? Use the value from the input file. As a technicality, the leapfrog method should be initialized with the velocity at time t = -Δt / 2, so we'll assume this is the value in the input file. In real codes, special care must be made to deal with this.

I'm a physicist. Why should I use the leapfrog method instead of the formula I derived in high school? The leapfrom 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 have completely different long-term behavior. For example, the earth would slowly spiral in to (or away from) the sun. For these reasons, symplectic methods are extremely popular for N-body calculations. You asked!

My planets 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. Note that Δx and Δy can be positive or negative. Do not even think of changing G to -G.

How should I compute x²? The simplest way is x*x. The ^ operator does not mean exponentiation in Java.

I can't compile Turtle.java. Any thoughts? Be sure you are running Java 1.4 and not some older version.

Input, Output, and Testing

Input. Copy the nbody directory from the COS 126 ftp site to your computer.

• Windows NT using IE. Click the ftp link above, then right click the nbody folder, and select Copy To Folder.

• OS X using Safari. Click the ftp link above to mount the ftp site, then drag-and-drop the nbody folder from the ftp mount.

• Arizona Unix cluster. Type the following command:

  cp -pr /u/cos126/files/nbody .
  

Compilation.  Your program must be named NBody.java. The capitalization is important. Compile your program with:

javac NBody.java

Execution.  To redirect standard input from a file, execute your program with:

java NBody < planets.txt

Submission

readme.txt. Use the following readme file template and answer any questions.

Submission.  Submit NBody.java. Don't forget to hit the "Run Script" button on the submission system to test that it compiles cleanly.

Possible Progress Steps

These are purely suggestions for how you might make progress. You do not have to follow these steps. 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-rest your code as you write it. Do not attempt to write the whole program at once - that is not how programming works since when the program doesn't work, you will have no idea where to find the errror. We promise that proactive testing will save you enormous amounts of time in the long run. Trust us! Also, if you get stumped or frustrated on some portion of the assignment, you should not hesitate to consult a preceptor.

• Review the Bouncing Ball program from lecture and the Student Precept List program from Section 2.5. Compile and execute it to make sure your system is configured properly.

• Copy the files from /u/cos126/files/nbody into your working directory. In particular, be sure you have Turtle.java, StdIn.java, and some sample input files.

• Read in the data file planets.txt using StdIn.java. To test that you read it in correctly, print the information back out using System.out.println. Do not even think of continuing until you have checked that you read in the data correctly.

• Now, read in the data file and store the information in several arrays. Let px[i], py[i], vx[i], vy[i], and mass[i] be real numbers which store the current position (x and y coordinates), velocity (x and y components), and mass of planet i. Let image[i] be a string which represents the filename of the GIF image used to display planet i. Print the information back out. Again, test that you read it in corretly by printing some of it back out.

• Step 3. Plot the background starfield.jpg image. Recall that Turtle.spot centers the image at the current location, so use Turtle.fly to pre-position the turtle. Test that it works. Now, write a loop to plot the N bodies using Turtle graphics. Be sure that you scale and translate the data so that it fits within the 512-by-512 Turtle graphics window. Don't forget to call Turtle.pause(20) after the loop to allow the turtle some some time rest and display the picture on the screen. Test it. 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 stuff you wrote in Step 3.) Since we haven't yet incorporated gravity, assume the acceleration acting on each planet is 0. Add an outer loop to repeat Steps 1 and 2 forever. Test it. 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 planet i. First, write a loop to initialize all the net forces to 0.0. Then write two nested for loops to calculate the net force exerted by planet j on planet i. Add these values to fx[i] and fy[i], but skip the case when i = 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.

• At the beginning of the simulation (outside of any loop), use Turtle.grunt to play the 2001 theme. Test it.

Extra Credit

One excellent idea is to add the Asteroid Belt, the Kuiper Belt, or Comets as in this N-body applet. Here is data for our solar system that the applet uses. You will need to to do some calculations to turn the data into our format and find appropriate images.

Enrichment

• Here's some information about our solar system.

• Here's a wealth of information on N-body simulation [pdf].

• Here's some of the N-body physics.

• 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. We will revisit one such approach in the final project.

• The same methods are also widely used in molecular dynamics, except that the 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 accout, e.g., to account for the covalent bonds in diamond crystals. To simulate bulk pieces of matter (10^23 atoms) scientists often use periodic boundary conditions since it's not possible to simulate nearly that many atoms.