Homework 2: NBody

A. Goals

The goal of the NBody assignment is to familiarize you with arrays, iteration, graphics, and Java command-line arguments. The specific goals are to:

At the end of this assignment, you will have produced a visual simulation of the gravitational interactions among particles in two dimensions.

B. Background

AND YOU THOUGHT YOU KNEW KNEWTON

In 1687, in his famous Principia, Sir Isaac Newton formulated the principles governing the motion of two particles under the influence of their mutual gravitational attraction. However, Newton was unable to solve the problem for three particles. Indeed, in general, systems of three or more particles can only be solved numerically.

Your challenge in this assignment is to write a program to simulate the motion of a designated number of particles in the plane that are mutually affected by gravitational forces, and animate the results. Such methods are widely used in astronomy, semiconductors, and fluid dynamics to study complex physical systems. Scientists also apply the same techniques to other pairwise interactions including Coulombic, Biot-Savart, and van der Waals.

We provide all the formulas you need in the assignment. Don't worry if you don't fully understand them – this is neither a physics nor a math course. Instead, focus on the computer science task: building up a complicated result through a series of simple statements.

C. Your program

We are going to use Newton's universal law of gravitation and Java to simulate the manner in which celestial bodies interact. Since we are most familiar with the solar system, the primary goal is to create a simulation showing the planets (Mercury through Mars) orbiting around the Sun. Here's what your end result will look like.

Each body in the simulation has three characteristics that define its movement: its position, its velocity (the speed and direction of movement), and its acceleration (the rate and direction of change of the particle’s velocity).

At each timestep in your simulation, you will draw all particles in their proper location. Then, for each pair of particles A and B, you will find the gravitational force that B exerts on A based on their respective masses and the distance between them using Newton’s Law of Gravitation. Since there is a force acting on A, Newton’s Second Law of Motion predicts that A will accelerate (i.e., its velocity will change). You will then compute A's new velocity. After finding the new velocities of all of the particles, you will compute their new positions, which will be displayed on the simulation.

D. Getting started

Create a Java class called NBody with an empty main() function. Give your file a header comment.

Download and compile PennDraw.java into the same folder as NBody.java.

Download nbody_data.zip and decompress it into the same folder as NBody.java. Note, this means the CONTENTS of the zipped folder (including the images and solarsystem.txt files should be in the SAME directory as your NBody.java file, NOT in a subdirectory called nbody-data).

E. Advice

A. Refresher: Using command-line arguments

Study the following program.

A Java program can take arguments when it is executed from the command-line at the Interactions pane. These arguments are passed to the main() function as a String array called args; thus, the first command-line argument is args[0], the second is args[1], and so on.

public class CommandLineArguments {
    public static void main(String[] args) {
        System.out.println(args.length);
        System.out.println(args[0]);
        System.out.println(Double.parseDouble(args[1]));
        System.out.println(Integer.parseInt(args[2]));
    }
}

Compile this program (you may cut and paste it into DrJava if you like) and run it by clicking Run on the toolbar. What happens?

Open the Interactions pane and enter these commands at the prompt:

Make sure you understand each line of code, and why you get the output and runtime errors that you do for each set of command-line arguments.

B. Setting up command-line arguments

Your NBody program will take three command-line arguments:

Declare these variables and set their values to the first, second, and third command-line arguments given to your program. These variables will be used later in your NBody program.

C. Checkpoint

Compile NBody, and run it with the arguments 15778000.0, 25000.0, and solarSystem.txt. Your program should not crash. If it does, use System.out.println() to isolate and rectify the bug. Any debugging code that you add should be commented out when you submit.

The specifications for the simulation are provided in a data file that provides the list of all particles, their masses, initial positions, images, etc. By storing the specification in a separate file, we can write a simulator that can load and simulate a variety of scenarios.

A. Input file format

You will be provided with a text file of the following format.

5        numParticles 
2.50e+11 radius 
5.97400e+24  1.49600e+11  0.00000e+00  0.00000e+00  2.98000e+04    earth.gif
6.41900e+23  2.27900e+11  0.00000e+00  0.00000e+00  2.41000e+04     mars.gif
3.30200e+23  5.79000e+10  0.00000e+00  0.00000e+00  4.79000e+04  mercury.gif
1.98900e+30  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00      sun.gif
4.86900e+24  1.08200e+11  0.00000e+00  0.00000e+00  3.50000e+04    venus.gif
m[]          px[]         py[]         vx[]         vy[]             img[]

In this file,

B. Setting up inStream

The booksite's standard library provides the file In.java to support reading from a file. Study StudentsFileProcessor.java, which is contained in the nbody_data.zip you downloaded in Part 0. It is an example of reading from a file. Compile StudentsFileProcessor.java, and run it from the DrJava Interactions Pane with the argument students.txt:

java StudentsFileProcessor students.txt

In your NBody class, declare and initialise a variable, inStream, as below:

In inStream = new In(filename); // creates a variable inStream to read from the file

C. Reading values from the file

Now that inStream is initialised, you can read from it using the following function calls. These functions behave identically to those in StdIn.

boolean b = inStream.isEmpty();     // boolean value that is true if there are no more values, false otherwise
int     i = inStream.readInt();     // reads in an int from inStream
double  d = inStream.readDouble();  // reads in a double from inStream
boolean b = inStream.readBoolean(); // reads in a boolean from inStream
String  s = inStream.readString();  // reads in a string from inStream
String  s = inStream.readLine();    // reads in an entire line from inStream
String  s = inStream.readAll();     // reads in the entire file from inStream

inStream will start reading from the beginning of the file. Each time a function like inStream.readDouble() is called, inStream interprets the next number as a double (an error will occur if it cannot be parsed to a double). The next time a read function is called, inStream moves to the next item in the file.

For example, say that a file, sample.txt, is as follows:

4 5

The code snippet

In inStream = new In("sample.txt");
int x = inStream.readInt();
double y = inStream.readDouble();

will set x to 4 and y to 5.0.

In your NBody program, declare an integer variable numParticles and a double variable radius. Using the above functions, assign to numParticles and radius the first and second values in the file, respectively.

D. Reading arrays from the file

Declare double arrays m, px, py, vx, and vy, and String array img, and initialise them to length numParticles. Use the inStream functions and a loop to read the values in the numParticles lines in the file.

E. Checkpoint

Copy the code below and paste it at the end of your main() function.

System.out.printf("%d\n", numParticles);
System.out.printf("%.2e\n", radius);
for (int i = 0; i < numParticles; i++) {
    System.out.printf("%12.5e %12.5e %12.5e %12.5e %12.5e %12s\n",
    m[i], px[i], py[i], vx[i], vy[i], img[i]);
}

Compile your program and run it with the same arguments: 15778000.0, 25000.0, and solarSystem.txt. Note that by pressing ↑ you can recall commands previously typed into the Interactions pane. The output from the Interactions pane should exactly match the text in the file solarSystem.txt. If your program crashes or if the output is different, use System.out.println()statements to isolate and rectify the bug. Any debugging code that you add should be commented out (but not deleted) when you submit.

A. Drawing the background

Use

PennDraw.setXscale(-radius, radius);
PennDraw.setYscale(-radius, radius);

to set the coordinates of the simulation window.

Use PennDraw.picture() to draw the image starfield.jpg in the centre of the window.

B. Drawing particles

Write a loop that draws each particle in the correct position, on top of the background image. The position of the ith particle is (px[i], py[i]), and the image to be used is img[i].

C. Checkpoint

Compiling and running your program now using the same command-line arguments 15778000.0, 25000.0, and solarSystem.txt should show five particles (the Sun, Mercury, Venus, Earth, and Mars), in a line, stationary, on a starfield background image.

NBody still

The motion of particles is governed by a set of equations known as Newton's Laws of Motion and Gravitation. (At the scale of the galaxy, we can treat planets as particles.) The equations contain two basic parts: velocity, which is the speed and direction that a particle is moving in, and acceleration, which is the change in speed and direction. Particles always move in a straight line at a constant speed unless some outside force (e.g. friction, gravity, or a push) accelerates it. In our simulation, the only force is gravity. First, however, you will implement a universe without gravity and test to make sure that it works.

Although we think of time as advancing continuously and of objects as moving continuously, computers are not good at simulating constant change. Instead, we settle for letting time jump forward in discrete time steps of Δt (timeStep, which your program already reads as the second command-line argument). This variable is called the time quantum and is measured in seconds. Particles will simply jump from the old position to the new position at each time step. By setting Δt to a very small value, we can get an accurate simulation, but the number of computations will be large. Setting it very large will speed up the simulation, but it will not be very accurate. This simulation scheme is called the leapfrog finite difference approximation scheme. It numerically integrates Newton's equations, and is the basis for most astrophysical simulations of gravitational systems.

A. The time loop

Write a loop that surrounds the code that you have to draw the background and particles. You will need to declare a variable that represents the elapsed time of your simulation. The elapsed time increases by timeStep in each iteration, and the simulation continues as long as the elapsed time is less than simulationTime.

Add a call to PennDraw.enableAnimation(30) before your time loop to enable animation at 30fps. You may adjust the frame rate to suit your tastes or to help you debug (decrease the frame rate a lot to see how the particles are moving from frame to frame; increase it to get to the end of the simulation faster).

Within your time loop, add a call to PennDraw.advance() after drawing all the particles.

Place the given code from Section 2E after the time loop.

B. Update particle positions

Inside the time loop, before redrawing the particles, update the position of each of the numParticles particles. For the original values of position: (px, py) the resultant values should be (px + Δt vx, py + Δt vy).

C. Checkpoint

Compiling and running your program now with the same arguments should result in the Sun staying put in the center of the window, and the four planets moving straight up and off the screen. This is because the initial velocities specified in solarSystem.txt all point straight up, and we have not yet incorporated the acceleration.

If the animation flickers or smears, check that you are drawing the background only once per iteration of the time loop.

D. Update velocities

force diagram

Outer loop – For every particle A (represented by the Earth in the image above, you must compute the forces exerted on it by every other particle B (represented in the image above by the Sun). This outer loop will use the net force on particle A, Fx and Fy, to calculate the acceleration of A and update its velocity. See above image for the formulas on how to calculate Fx and Fy. We suggest that you declare variables for each calculation you see in the image. (Your variable names should respect the camel case naming conventions of Java variables rather than the the convential symbols used in physics.)

Inner loop – The net force exerted on particle A should be the sum of the forces exerted on it by each other particle B. Define variables to hold the x- and y-components of the force exerted by particle B on A. Hint: Each iteration of the inner loop should increase the net force Fx and Fy on particle A by the force exerted by particle B on A in the x- and y-directions. When the inner loop finishes for particle A, the values of Fx and Fy should be the net force components for particle A that will be used to calculate its acceleration. (Your variable names should respect the camel case naming conventions of Java variables rather than the the convential symbols used in physics.)

Δx and Δy – Be careful computing Δx and Δy. These are the distances in the x- and y-directions from the first planet to the second one, and can be negative. You should not take the absolute value of the differences between their x- and y-positions.

Gravity – Define a variable G for the gravitational constant so the formulas will be easier to read. Initialize it to the value 6.67e-11. When calculating the force between particle A and B, use the following order of operations to avoid NaN values: ((G * mA) / (d * d)) * mB

Acceleration – By Newton's Second Law, for each particle, the resultant acceleration components of that particle are ax = Fx / m and ay = Fy / m. Compute the x- and y-components of the particle's acceleration.

Velocity – Using the acceleration you just computed, update the velocity (vx, vy) to (vx + Δt ax, vy + Δt ay). This gives you the new direction and speed of each particle.

Order of computation – It is imperative to update the velocity of every particle before updating any positions. Even though your program computes these velocities one after the other, you have to maintain the illusion that all planets are moving at the same time. For instance, if you update the Earth's position before computing Mars' velocity, then you will end up computing Mars' velocity based on Earth's future position, rather than its current position.

E. Common Problems

Anti-Gravity – If you flip the order of your particles when you compute Δx and Δy, your particles will repel rather than retract each other. For the solarSystem.txt example, you might end up with something that looks like this:



Planets jump to the upper-left corner – Check your loop for computing the pairwise forces between particles. Are you computing the force exerted by a particle on itself (e.g. Earth's gravitational pull on itself)? If you are, do the computation by hand to see exactly what it would give you. That should help you understand why this is a bad idea, and why the planets jump the way they do. (Video is slowed for illustration.)



Planets go around the sun, then zoom off the lower left corner – You are not summing up the pairwise forces correctly. Run your program with a small value of simulationTime (e.g. 51000) so it runs for only two or three time steps. Print out the value of every variable before and after you change its value (or right after you first initialize it). Include enough information in the print statement so you can tell which variable you are printing out. Make sure to print all loop variables as well so you can tell from the output exactly when each print statement occured. Run the program, and examine the output. Ask yourself each time a variable changes whether its starting and new values make sense. (You don't need to work out the exact values of the variables to figure out if they make sense in this case.) Once you've found the problem, you may want to move the declarations of some variables farther inside your loops to the first point you need them. Declaring variables as late as possible is one strategy to avoid the type of bug you are encountering. (Video is slowed for illustration.)



Planets are flickering – You are either not calling PennDraw.advance() at all, or calling it after drawing each particle instead of after drawing all particles. This means each planet appears to move as soon as you draw it, instead of all the planets appearing to move in sync.

Planets are flickering, and only one shows at a time – You are redrawing the background before drawing each planet, instead of redrawing it once per timestep. As a result, you keep covering up each planet as soon as you draw it.

Planets "smear" accross the screen – You are drawing the background once at the beginning of the programming instead of once each time step. PennDraw does not erase the window automatically. In this program, we're never erasing it at all. Instead, we keep covering the old planet positions up with the starfield image, which happens to be exactly the same size as the window. If you don't keep redrawing the background before drawing the planets at their new positions, the old planets are still there, causing the smearing effect.

F. Checkpoint

Compile and run your program. The particles should now orbit the Sun, and at the end of the simulation, your program should print the final state of the system.

The following are the values you should obtain for the final state of the system with various values of simulationTime and timeStep. When testing, you should start with a single time step, making sure that it is working completely. Then try a small number of time steps, again making certain that it is working completely, before going on to a full year. It is much easier to diagnose problems for smaller time frames.

java NBody 25000.0 25000.0 solarSystem.txt (one step)

5
2.50e+11
5.97400e+24  1.49596e+11  7.45000e+08 -1.48201e+02  2.98000e+04    earth.gif
6.41900e+23  2.27898e+11  6.02500e+08 -6.38597e+01  2.41000e+04     mars.gif
3.30200e+23  5.78753e+10  1.19750e+09 -9.89331e+02  4.79000e+04  mercury.gif
1.98900e+30  3.30867e+01  0.00000e+00  1.32347e-03  0.00000e+00      sun.gif
4.86900e+24  1.08193e+11  8.75000e+08 -2.83294e+02  3.50000e+04    venus.gif

java NBody 50000.0 25000.0 solarSystem.txt (two steps)

5
2.50e+11
5.97400e+24  1.49589e+11  1.48998e+09 -2.96404e+02  2.97993e+04    earth.gif
6.41900e+23  2.27895e+11  1.20500e+09 -1.27720e+02  2.40998e+04     mars.gif
3.30200e+23  5.78258e+10  2.39449e+09 -1.97887e+03  4.78795e+04  mercury.gif
1.98900e+30  9.92618e+01  2.81978e-01  2.64700e-03  1.12791e-05      sun.gif
4.86900e+24  1.08179e+11  1.74994e+09 -5.66597e+02  3.49977e+04    venus.gif

java NBody 75000.0 25000.0 solarSystem.txt (three steps)

5
2.50e+11
5.97400e+24  1.49578e+11  2.23493e+09 -4.44605e+02  2.97978e+04    earth.gif
6.41900e+23  2.27890e+11  1.80748e+09 -1.91579e+02  2.40995e+04     mars.gif
3.30200e+23  5.77516e+10  3.59045e+09 -2.96820e+03  4.78386e+04  mercury.gif
1.98900e+30  1.98523e+02  1.12801e+00  3.97047e-03  3.38411e-05      sun.gif
4.86900e+24  1.08158e+11  2.62477e+09 -8.49891e+02  3.49931e+04    venus.gif

java NBody 31557600.0 25000.0 solarSystem.txt (one year)

5
2.50e+11
5.97400e+24  1.49594e+11 -1.65312e+09  3.29492e+02  2.97978e+04    earth.gif
6.41900e+23 -2.21530e+11 -4.92633e+10  5.18050e+03 -2.36404e+04     mars.gif
3.30200e+23  3.47713e+10  4.57516e+10 -3.82694e+04  2.94146e+04  mercury.gif
1.98900e+30  5.94260e+05  6.23570e+06 -5.85686e-02  1.62847e-01      sun.gif
4.86900e+24 -7.37309e+10 -7.93909e+10  2.54335e+04 -2.39734e+04    venus.gif

You should also test your program against some of the other universes that we provided. Ensure that your program does not misbehave when given more or fewer than five particles.

G. Music

This task is not for credit; it is just for fun. (Depending on the individual configuration of your personal computer, this may not work correctly.) Add the following line somewhere before your time loop: 

StdAudio.play("2001.mid");

Run your program again, and let us know in the readme whether you have sound. The music is the fanfare to Also sprach Zarathustra by Richard Strauss; it was popularised as the key musical motif in Stanley Kubrick’s 1968 film 2001: A Space Odyssey.

A. Extra credit

Submit a .zip file called extra.zip containing an alternate universe (in our input format) along with the necessary image files. If its behavior is sufficiently interesting, we'll award extra credit. Your submission must be in a zip file, even if there are no images, so that our grading scripts can handle it correctly.

B. Technical details on the simulation

This simulation method is called the leapfrog finite difference approximation scheme; it is based on a numeric integration of Newton’s equations, and is the basis for most astrophysical simulations of gravitational systems.

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.

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, (px, py) represents the position at time t whilst (vx, vy) represents the velocity at time t - Δt / 2. Interpreting the position and velocity in this way, the updated position (px + Δt vx, py + Δt vy) 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 in solarSystem.txt. 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.

A. Readme

Complete readme_nbody.txt in the same way that you have done for previous assignments.

B. Submission

Submit NBody.java and readme_nbody.txt on the course website. If you completed the extra credit, also submit extra.zip.