2. N-Body Simulation

The key to composing a large program is decomposing it into smaller steps that can be implemented and tested incrementally. Here is the decomposition we recommend for this assignment.

0. Start with comments
   
1. Parse command-line arguments
   
2. Read universe from standard input
   
3. Initialize standard drawing
   
4. Play music on standard audio
   
5. Simulate the universe (main time loop)
   

• A Calculate net forces
  
• B. Update velocities and positions
  
• C. Draw universe to standard drawing
  

6. Print universe to standard output

These are the order in which the components will appear in your code. Although your final code will appear in order 1–6, we recommend implementing these steps in the order 0, 1, 2, 6, 3, 4, 5, 5B, 5C, 5A. Why? Because this order will facilitate testing and debugging along the way.

Step 0: Start with comments

This outline should help you organize your program.

// Step 1. Parse command-line arguments.

// Step 2. Read the universe from standard input.

// Step 3. Initialize standard drawing.

// Step 4. Play music on standard audio.

// Step 5. Simulate the universe (main time loop).

  // Step 5A. Calculate net forces.
  // Step 5B. Update velocities and positions.
  // Step 5C. Draw universe to standard drawing.

// Step 6. Print universe to standard output.

Step 1: Parse command-line arguments

• Parse the two command-line arguments Τ and Δt and store their values. Name the first command-line argument tau or stoppingTime since you should not begin a Java variable name with an uppercase letter.
• Print the variables to make sure you parsed them correctly (and in the specified order). For example:

 java-introcs NBody 10 1					java-introcs NBody 157788000.0 25000.0 < planets.txt
T = 10.0							T = 1.57788E8
dt = 1.0							dt = 25000.0

• Remove the print statements after you confirm you are parsing the command-line arguments correctly.

Step 2: Read universe from standard input

• Read the number of bodies n as an int from standard input.
• Read the radius of the universe radius as a double from standard input.
• Create six (6) parallel arrays, each of length n, to store the six (6) pieces of information characterizing a body. Each array px[i], py[i], vx[i], vy[i], and mass[i] is a double arrays and usused to store the current position (\(x\)- and \(y\)-coordinates), velocity (\(x\)- and \(y\)-components), and mass of particle i, respectively; let image[i] be a String that represents the filename of the image used to display particle i.
• Read data associated with each body and store it in parallel arrays.

Hint!
Recall the Students.java exercise from your precept.

• Run your program with planets.txt.

java-introcs NBody 157788000.0 25000.0 < planets.txt

Step 6: Print universe to standard output

• Print the universe in the specified format.
• Our input files were created with the following StdOut.printf() statements.

StdOut.printf("%d\n", n);
StdOut.printf("%.2e\n", radius);
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]);
}

Per the course collaboration policy, copying or adapting code that is not yours is permitted only if it comes from the course materials. If you do so, you must cite any code that you copy or adapt (with the exception of code that is included with the assignment). So, you are free to copy or adapt this code fragment, with or without citation.

• Run your program with planets.txt. You should see exactly the output below.

java-introcs NBody 0.0 25000.0 < 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

• Test your program on another input file.

java-introcs NBody 0.0 1.0 < 3body-zero-gravity.txt
3
5.12e+02
 0.0000e+00  0.0000e+00  1.0000e+00  1.0000e+00  1.0000e-30    earth.gif
 1.2800e+02  0.0000e+00  2.0000e+00  1.0000e+00  1.0000e-40    venus.gif
 0.0000e+00  1.2800e+02  1.0000e+00  2.0000e+00  1.0000e-50     mars.gif

• Submit your code to TigerFile and click the Check Submitted Files button. If everything is correct so far (reading and printing the universe), your code should pass Tests 1–7. Do not continue to the next step until it does.

Step 3: Initialize standard drawing

• The default \(x\)- and \(y\)-scale supports coordinates between 0 and 1. Set the scale of the standard drawing window so that (0, 0) is at the center, (−radius, −radius) is at the lower-left corner, and (+radius, +radius) is at the upper-right corner.

  Hint!
  Recall StdDraw.setXscale() and StdDraw.setYscale()

• Call StdDraw.enableDoubleBuffering() to enable double buffering and support animation.

Step 4: Play music on standard audio

• Call StdAudio.playInBackground("2001.wav") to play the background music.

Step 5: Simulate the universe (the big time loop)

• Create the time simulation loop. Each iteration will perform a single time step of the simulation, starting at \(t = 0.0\), and continuing in \(Δt\) increments as long as \(t < Τ\). This code goes between the code where you read the input and the code where you print the output.
• Test your program by printing the time \(t\) at each time step.

java-introcs NBody 23.0 2.5 < planets.txt		java-introcs NBody 25.0 2.5 < planets.txt
t = 0.0							t = 0.0
t = 2.5							t = 2.5
t = 5.0							t = 5.0
t = 7.5							t = 7.5
t = 10.0						t = 10.0
t = 12.5						t = 12.5
t = 15.0						t = 15.0
t = 17.5						t = 17.5
t = 20.0						t = 20.0
t = 22.5						t = 22.5

• Once you are confident that your time loop is correct, remove the code that prints the time \(t\) at each time step.

Steps 5B, 5C, and 5A (below) will go inside this time loop.

FAQ
In which order should I implement these 3 sub-steps? 5B, 5C, 5A because calculating forces is hardest.

FAQ
Can I interleave steps 5A, 5B, and 5C? No. Not only is it bad design, but it ruins the physics: you need the position of all bodies at time \(t\) to calculate the forces between each pair of bodies, and updating some positions to time \(t + ∆t\) would mess up the forces for time \(t\).

Step 5B: Update the velocities and positions

• Since we have not implemented Step 5A yet, the forces \(F_x\) and \(F_y\) are 0, so we can initially assume that the accelerations \(a_x\) and \(a_y\) are zero and the velocities are constant.
• Write a loop to calculate the new velocity and position for each particle.
  • Update the velocity of each body: \(v_x = v_x + a_x∆t\) and \(v_y = v_y + a_y∆t\) (note for now \(a_x\) and \(a_y\) are 0).
  • Update the position of each body: \(p_x = p_x + v_x∆t\) and \(p_y = p_y + v_y∆t\).
• Test on an artificial universe, where it is easy to check the expected results by hand.

java-introcs NBody 1 1 < 3body-zero-gravity.txt
3
5.12e+02
 1.0000e+00  1.0000e+00  1.0000e+00  1.0000e+00  1.0000e-30    earth.gif
 1.3000e+02  1.0000e+00  2.0000e+00  1.0000e+00  1.0000e-40    venus.gif
 1.0000e+00  1.3000e+02  1.0000e+00  2.0000e+00  1.0000e-50     mars.gif

java-introcs NBody 192 1 < 3body-zero-gravity.txt
3
5.12e+02
 1.9200e+02  1.9200e+02  1.0000e+00  1.0000e+00  1.0000e-30    earth.gif
 5.1200e+02  1.9200e+02  2.0000e+00  1.0000e+00  1.0000e-40    venus.gif
 1.9200e+02  5.1200e+02  1.0000e+00  2.0000e+00  1.0000e-50     mars.gif

• Submit your code via TigerFile and click the Check Submitted Files button. If everything is correct so far (ignoring the gravitational forces and drawing), your code should pass Tests 1–10. Do not continue to the next step until it does.

Step 5C: Draw universe to standard drawing

• Draw each particle at its current position to standard drawing, and repeat this process at each time step until the designated stopping time. By displaying this sequence of snapshots (or frames) in rapid succession, you will create the illusion of movement.
• After each time step, draw the background image starfield.jpg; redraw all the particles in their new positions; and control the animation speed.
  • Use StdDraw.picture(x, y, filename) to draw the background image starfield.jpg, centered at (0, 0).
  • Then, write a loop to display the n particles.
  • Call StdDraw.show() to display results on screen.
  • Call StdDraw.pause(20) to control animation speed, which produces about 50 frames per second.
• Run the program with planets.txt. You should now see the four planets moving off the screen in a straight line, with constant velocity.

java-introcs NBody 157788000.0 25000.0 < planets.txt

• If the planets are flickering, be sure you are using double buffering, following the template in BouncingBallsDeluxe.java. In particular, call StdDraw.enableDoubleBuffering() once at the beginning of the program, and call StdDraw.show() and StdDraw.pause(20) at the end of each time step (frame), not after each call to StdDraw.picture().

• Run your program on another input file, such as kaleidoscope.txt, which produces this animation:

java-introcs NBody 157788000.0 25000.0 < kaleidoscope.txt

• Submit your code via TigerFile and click the Check Submitted Files button. If everything is correct so far (ignoring the gravitational forces), your code should pass Tests 1–17. Do not continue to the next step until it does.

Step 5A: Calculate the forces

• This step is the most difficult one to implement. Let’s review:
  • Pairwise force. Recall this diagram from above:

• In Java, you can represent the gravitational constant using scientific notation: 6.67e-11 (or 6.67E-11)

• Net force. The principle of superposition says that the net force acting on a particle in the \(x\)- or \(y\)-direction is the sum of the pairwise forces acting on the particle in that direction. Calculating the force (between all pairs of bodies at time \(t\)) involves adding all the pairwise forces. For example, the forces on earth (at time \(t\)) are: $$ \vec{F}_{earth} = \vec{F}_{mars \to earth} + \vec{F}_{mercury \to earth} + \vec{F}_{sun \to earth} + \vec{F}_{venus \to earth} $$

• To calculate the net force acting on each body, you will need two additional arrays fx[i] and fy[i] to store the net force acting on body i.

• The net force (arrays fx[i] and fy[i]) must be initialized to 0 at each time step.

• Use a nested loop pattern 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 (you don’t want to calculate the force of body i on itself).

  Hint!
  Before calculating the net forces, can you implement a nested loop pattern that enumerates all pairs of bodies? For example: what pairs are not printed?

n == 4                  n == 5
0-1 0-2 0-3             0-1 0-2 0-3 0-4
1-0 1-2 1-3             1-0 1-2 1-3 1-4
2-0 2-1 2-3             2-0 2-1 2-3 2-4
3-0 3-1 3-2             3-0 3-1 3-2 3-4
                        4-0 4-1 4-2 4-3

• Once you have these values computed, go back to Step 5B and use them to compute the acceleration (instead of assuming it is zero).
  • Calculate the acceleration \( (a_x, a_y) \) at time \(t\) using the net force computed in Step 5A and Newton’s second law of motion: \( a_x = \frac{F_x}{m} \) and \( a_y = \frac{F_y}{m} \).
  • The code from Step 5B should then compute the new velocities and positions.
• Test your program on several data files.
• Submit your code via TigerFile and click the Check Submitted Files button. Your code should now pass all of the tests!

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!

Make sure you are using the correct index variables. It’s very easy to confuse an i with a j, especially when copying and pasting Java statements.

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 using 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 + v_xΔt, p_y+v_yΔt) \) 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.

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.

Kepler’s Orrery uses an n-body simulator to compose and play music.