N-body Simulations - Code Part 1

Object-oriented programming

In general, programming languages have built in objects, such as integers, float point numbers, booleans, arrays, etc. that can be used very easily by the user. These allow programmers to code virtually whatever they need, but it is often easier to create objects that better suit the application that is being coded. The following examples are written in Java, so a small background would be useful. Even though the syntax may be slightly different, other programming languages such as C++ and Python can accomplish identical tasks.

It would be nice to declare a body, b, with the following parameters:
Body b = new Body(double rx, double ry, double vx, double vy, double mass, Color color)
This is accomplished, in Java, using the code below. Comments are prefaced by //.

import java.awt.Color;

public class Body {
  private static final double G = 6.673e-11;   // gravitational constant
  private static final double solarmass=1.98892e30;
  
  public double rx, ry;       // holds the cartesian positions
  public double vx, vy;       // velocity components 
  public double fx, fy;       // force components
  public double mass;         // mass
  public Color color;         // color (for fun)
  
  // create and initialize a new Body
  public Body(double rx, double ry, double vx, double vy, double mass, Color color) {
    this.rx    = rx;
    this.ry    = ry;
    this.vx    = vx;
    this.vy    = vy;
    this.mass  = mass;
    this.color = color;
  }
  
  // update the velocity and position using a timestep dt
  public void update(double dt) {
    vx += dt * fx / mass;
    vy += dt * fy / mass;
    rx += dt * vx;
    ry += dt * vy;
  }
  
  // returns the distance between two bodies
  public double distanceTo(Body b) {
    double dx = rx - b.rx;
    double dy = ry - b.ry;
    return Math.sqrt(dx*dx + dy*dy);
  }
  
  // set the force to 0 for the next iteration
  public void resetForce() {
    fx = 0.0;
    fy = 0.0;
  }
  
  // compute the net force acting between the body a and b, and
  // add to the net force acting on a
  public void addForce(Body b) {
    Body a = this;
    double EPS = 3E4;      // softening parameter (just to avoid infinities)
    double dx = b.rx - a.rx;
    double dy = b.ry - a.ry;
    double dist = Math.sqrt(dx*dx + dy*dy);
    double F = (G * a.mass * b.mass) / (dist*dist + EPS*EPS);
    a.fx += F * dx / dist;
    a.fy += F * dy / dist;
  }
  
  // convert to string representation formatted nicely
  public String toString() {
    return "" + rx + ", "+ ry+ ", "+  vx+ ", "+ vy+ ", "+ mass;
  }
}

This defines the Body object that we will use to store information about bodies in our simulation. This file would be saved in a file Body.java and compiled into a class Body.class. The class file should reside in the same directory as all other class files.

On the next page we will create the main class that uses these methods to run the simulation.

N-body Simulations- Code Part 2

Using our Body object

Recall that by creating a new class, we have a Body object that can be called using:

Body b = new Body(double rx, double ry, double vx, double vy, double mass, Color color)
Now, we want to use the brute-force method to update the positions and velocities of the bodies. Here's the code for an applet that we'll see at the end of this tutorial.

BruteForce.java

// tell the compiler where to find the methods you will use.
// required when you create an applet
import java.applet.*;
// required to paint on screen
import java.awt.*;
import java.awt.event.*;


//Start the applet and define a few necessary variables
public class BruteForce extends Applet {
    public int N = 100;
    public Body bodies[]= new Body[10000];
    public TextField t1;
    public Label l1;
    public Button b1;
    public Button b2;
    public boolean shouldrun=false;
    

    
// The first time we call the applet, this function will start
  public void init()
  {
        startthebodies(N);
        t1=new TextField("100",5);
        b2=new Button("Restart");
        b1=new Button("Stop");
        l1=new Label("Number of bodies:");
        ButtonListener myButtonListener = new ButtonListener();
        b1.addActionListener(myButtonListener);
        b2.addActionListener(myButtonListener);
        add(l1);
        add(t1);
        add(b2);
        add(b1);
  }
  
  
// This method gets called when the applet is terminated. It stops the code
  public void stop()
  {
    shouldrun=false;
  }
  
  
//Called by the applet initally. It can be executed again by calling repaint();
  public void paint(Graphics g)
  {
    g.translate(250,250); //Originally the origin is in the top right. Put it in its normal place
    
//check if we stopped the applet, and if not, draw the particles where they are
    if (shouldrun) {
      for(int i=0; i<N; i++) {
        g.setColor(bodies[i].color);
        g.fillOval((int) Math.round(bodies[i].rx*250/1e18),(int) Math.round(bodies[i].ry*250/1e18),8,8);
      }
      //go through the Brute Force algorithm (see the function below)
      addforces(N);
      //go through the same process again until applet is stopped
      repaint();
    }
  }
  
  //the bodies are initialized in circular orbits around the central mass.
  //This is just some physics to do that
  public static double circlev(double rx, double ry) {
    double solarmass=1.98892e30;
    double r2=Math.sqrt(rx*rx+ry*ry);
    double numerator=(6.67e-11)*1e6*solarmass;
    return Math.sqrt(numerator/r2);
  }
  
  //Initialize N bodies with random positions and circular velocities
  public void startthebodies(int N) {
    double radius = 1e18;        // radius of universe
    double solarmass=1.98892e30;
    for (int i = 0; i < N; i++) {
      double px = 1e18*exp(-1.8)*(.5-Math.random());
      double py = 1e18*exp(-1.8)*(.5-Math.random());
      double magv = circlev(px,py);
      
      double absangle = Math.atan(Math.abs(py/px));
      double thetav= Math.PI/2-absangle;
      double phiv = Math.random()*Math.PI;
      double vx   = -1*Math.signum(py)*Math.cos(thetav)*magv;
      double vy   = Math.signum(px)*Math.sin(thetav)*magv;
      //Orient a random 2D circular orbit
     
           if (Math.random() <=.5) {
              vx=-vx;
              vy=-vy;
            } 
           
      double mass = Math.random()*solarmass*10+1e20;
      //Color the masses in green gradients by mass
      int red     = (int) Math.floor(mass*254/(solarmass*10+1e20));
      int blue   = (int) Math.floor(mass*254/(solarmass*10+1e20));
      int green    = 255;
      Color color = new Color(red, green, blue);
      bodies[i]   = new Body(px, py, vx, vy, mass, color);
    }
    //Put the central mass in
    bodies[0]= new Body(0,0,0,0,1e6*solarmass,Color.red);//put a heavy body in the center
    
  }
  
  //Use the method in Body to reset the forces, then add all the new forces
  public void addforces(int N) {
    for (int i = 0; i < N; i++) {
      bodies[i].resetForce();
      //Notice-2 loops-->N^2 complexity
      for (int j = 0; j < N; j++) {
        if (i != j) bodies[i].addForce(bodies[j]);
      }
    }
    //Then, loop again and update the bodies using timestep dt
    for (int i = 0; i < N; i++) {
      bodies[i].update(1e11);
    }
  }
   public static double exp(double lambda) {
        return -Math.log(1 - Math.random()) / lambda;
    }
   public boolean action(Event e,Object o)
   {
     N = Integer.parseInt(t1.getText());
     if (N>1000) {
       t1.setText("1000");
       N=1000;
     }
     
       startthebodies(N);
       repaint();
     
     return true;
   }
   public class ButtonListener implements ActionListener{

    public void actionPerformed(ActionEvent evt) 
    {
        // Get label of the button clicked in event passed in
        String arg = evt.getActionCommand();    
        if (arg.equals("Restart"))
        {
          shouldrun=true;
               N = Integer.parseInt(t1.getText());
     if (N>1000) {
       t1.setText("1000");
       N=1000;
     }
     
       startthebodies(N);
       repaint();
        }
        else if (arg.equals("Stop")) 
            stop();
    }
}

}

And that's it! It's pretty simple once we have the Body object defined. This is easily embedded in a webpage using the HTML tag.

<applet code>

You can skip to the Applets page to see this one in action, or head to the Next section of programming, where we will look at the Barnes-Hut algorithm.

N-body Simulations- Code Part 3

Implementing the Barnes-Hut Algorithm: A Quadrant Class

In the Introduction we talked about an object called a Barnes-Hut tree. This is one of two additional parts to the Barnes-Hut algorithm. The other is the concept of a quadrant (in 2D) or an octant (in 3D). For coding simplicity, we will create a quadrant-based Barnes-Hut tree. It is not difficult to add additional nodes if need be.

To create a quadrant object, the code goes as follows (again, comments preceded by //):

//A quadrant object that can self-subdivide. Important for creating a Barnes-Hut tree, since it holds quadrants.

public class Quad {
  
  private double xmid, ymid, length;
  //Create a square quadrant with a given length and midpoints (xmid,ymid)
  public Quad(double xmid, double ymid, double length) {
    this.xmid=xmid;
    this.ymid=ymid; 
    this.length=length;
  }
  
  //How long is this quadrant?
  public double length() {
    return length;
  }
  
  //Check if the current quadrant contains a point
  public boolean contains(double xmid, double ymid) {
    if (xmid<=this.xmid+this.length/2.0 && xmid>=this.xmid-this.length/2.0 && ymid<=this.ymid+this.length/2.0 && ymid>=this.ymid-this.length/2.0) {
      return true;
    }
    else {
      return false;
    }
  }
  //Create subdivisions of the current quadrant
  
  // Subdivision: Northwest quadrant
  public Quad NW() {
    Quad newquad = new Quad(this.xmid-this.length/4.0, this.ymid+this.length/4.0,this.length/2.0);
    return newquad;
  }
  
  // Subdivision: Northeast quadrant
  public Quad NE() {
    Quad newquad = new Quad(this.xmid+this.length/4.0, this.ymid+this.length/4.0,this.length/2.0);
    return newquad;
  }
  
  // Subdivision: Southwest quadrant
  public Quad SW() {
    Quad newquad = new Quad(this.xmid-this.length/4.0, this.ymid-this.length/4.0,this.length/2.0);
    return newquad;
  }
  
  // Subdivision: Southeast quadrant
  public Quad SE() {
    Quad newquad = new Quad(this.xmid+this.length/4.0, this.ymid-this.length/4.0,this.length/2.0);
    return newquad;
  }
  
}
One more thing that is useful to update now is the Body class so we can check if the body is in a certain quadrant. To do this, all we have to do is add the following method, which takes advantage of the methods we just defined in the Quad class.

public boolean in(Quad q) {
    return q.contains(this.rx,this.ry);
  }
}

Now we have a data structure to hold Bodies, but what we really need is the Barnes-Hut tree that hold quadrants and contain information about the bodies they contain. The next section will do just that.

N-body Simulations- Code Part 4
Implementing the Barnes-Hut Algorithm: A Barnes-Hut Tree

Now, we need to create a Barnes-Hut tree that will hold the quadrants and bodies, and allow us to approximate far-away bodies by their total and center of mass.

The code is fairly complex and its methods are non-trivial, so it is necessary to really look through the code to see what is happening. The complexity draws for the recursive methods that are mentioned in the comments.

To create a Barnes-Hut tree, the code goes as follows (again, comments preceded by //):

import java.awt.Color;

public class BHTree {
  private Body body;     // body or aggregate body stored in this node
  private Quad quad;     // square region that the tree represents
  private BHTree NW;     // tree representing northwest quadrant
  private BHTree NE;     // tree representing northeast quadrant
  private BHTree SW;     // tree representing southwest quadrant
  private BHTree SE;     // tree representing southeast quadrant
  
  //Create and initialize a new bhtree. Initially, all nodes are null and will be filled by recursion
  //Each BHTree represents a quadrant and a body that represents all bodies inside the quadrant
  public BHTree(Quad q) {
    this.quad=q;
    this.body=null;
    this.NW=null;
    this.NE=null;
    this.SW=null;
    this.SE=null;
  }
  //If all nodes of the BHTree are null, then the quadrant represents a single body and it is "external"
  public Boolean isExternal(BHTree t) {
    if (t.NW==null && t.NE==null && t.SW==null && t.SE==null) return true;
    else return false;
  }
  //We have to populate the tree with bodies. We start at the current tree and recursively travel through the branches
  public void insert(Body b) {
    //If there's not a body there already, put the body there.
    if (this.body==null) {
      this.body=b;
    }
    //If there's already a body there, but it's not an external node
    //combine the two bodies and figure out which quadrant of the 
    //tree it should be located in. Then recursively update the nodes below it.
    else if (this.isExternal(this)==false) {
      this.body=b.add(this.body,b);
      
      Quad northwest = this.quad.NW();
      if (b.in(northwest)) {
        if (this.NW==null) {this.NW= new BHTree(northwest);}
        NW.insert(b);    
      }
      else {
        Quad northeast = this.quad.NE();
        if (b.in(northeast)) {
          if (this.NE==null) {this.NE= new BHTree(northeast);}
          NE.insert(b);
        }
        else {
          Quad southeast = this.quad.SE();
          if (b.in(southeast)) {
            if (this.SE==null) {this.SE= new BHTree(southeast);}
            SE.insert(b);
          } 
          else {
            Quad southwest = this.quad.SW();
            if(this.SW==null) {this.SW= new BHTree(southwest);}
            SW.insert(b);
          }
        }
      }
    }
    //If the node is external and contains another body, create BHTrees
    //where the bodies should go, update the node, and end 
    //(do not do anything recursively)
    else if (this.isExternal(this)) {
      Body c = this.body;
      Quad northwest = this.quad.NW();
      if (c.in(northwest)) {
        if (this.NW==null) {this.NW= new BHTree(northwest);}
        NW.insert(c);    
      }
      else {
        Quad northeast = this.quad.NE();
        if (c.in(northeast)) {
          if (this.NE==null) {this.NE= new BHTree(northeast);}
          NE.insert(c);
        }
        else {
          Quad southeast = this.quad.SE();
          if (c.in(southeast)) {
            if (this.SE==null) {this.SE= new BHTree(southeast);}
            SE.insert(c);
          } 
          else {
            Quad southwest = this.quad.SW();
            if(this.SW==null) {this.SW= new BHTree(southwest);}
            SW.insert(c);
          }
        }
      }
      this.insert(b);
    }
  }
  //Start at the main node of the tree. Then, recursively go each branch
  //Until either we reach an external node or we reach a node that is sufficiently
  //far away that the external nodes would not matter much.
  public void updateForce(Body b) {
    if (this.isExternal(this)) {
      if (this.body!=b) b.addForce(this.body);
    }
    else if (this.quad.length()/(this.body.distanceTo(b))<2){
      b.addForce(this.body);
    }
    else {
      if (this.NW!=null) this.NW.updateForce(b);
      if (this.SW!=null) this.SW.updateForce(b);
      if (this.SE!=null) this.SE.updateForce(b);
      if (this.NE!=null) this.NE.updateForce(b);
    }
  }
  // convert to string representation for output
  public String toString() {
    if(NE != null || NW!=null || SW!=null || SE!=null) return "*" + body + "\n" + NW + NE + SW + SE;
    else           return " " + body + "\n";
  }
}

Now these are all the data types we need to run the algorithm. It's obvious that they are much less trivial than the Body data type, but they are essential if we want to cut down run time. Next, we put it all together in our main class that actually runs the simulation.

N-body Simulations- Code Part 5

Implementing the Barnes-Hut Algorithm: Main Method

Now that we have all our objects described, it is fairly straight-forward to run the simulation. Below is the commented code--the function that has been changed is addforces(int N), which has been updated to use the Barnes-Hut method. Other than that difference, our object-oriented approach has allowed this code to look almost identical to the Brute Force case. One big difference is that now, instead of invoking a loop to add all the forces on a body, we now just call the updateForce(Body b) method from our BHTree class, which takes, on average, Log(N) time to run, versus a loop that would take N.

// tell the compiler where to find the methods you will use.
// required when you create an applet
import java.applet.*;
// required to paint on screen
import java.awt.*;
import java.awt.event.*;

//Start the applet and define a few necessary variables
public class BarnesHut extends Applet {
    public int N = 100;
    public Body bodies[]= new Body[10000];
    public TextField t1;
    public Label l1;
    public Button b1;
    public Button b2;
    public boolean shouldrun=false;
    Quad q = new Quad(0,0,2*1e18);

    
// The first time we call the applet, this function will start
  public void init()
  {
        startthebodies(N);
        t1=new TextField("100",5);
        b2=new Button("Restart");
        b1=new Button("Stop");
        l1=new Label("Number of bodies:");
        ButtonListener myButtonListener = new ButtonListener();
        b1.addActionListener(myButtonListener);
        b2.addActionListener(myButtonListener);
        add(l1);
        add(t1);
        add(b2);
        add(b1);
  }
  
  
// This method gets called when the applet is terminated. It stops the code
  public void stop()
  {
    shouldrun=false;
  }
  
  
//Called by the applet initally. It can be executed again by calling repaint();
  public void paint(Graphics g)
  {
    g.translate(250,250); //Originally the origin is in the top right. Put it in its normal place
//check if we stopped the applet, and if not, draw the particles where they are
    if (shouldrun) {
      for(int i=0; i<N; i++) {
        g.setColor(bodies[i].color);
        g.fillOval((int) Math.round(bodies[i].rx*250/1e18),(int) Math.round(bodies[i].ry*250/1e18),8,8);
      }
      //go through the Barnes-Hut algorithm (see the function below)
      addforces(N);
      //repeat
      repaint();
    }
  }
  
  //the bodies are initialized in circular orbits around the central mass.
  //This is just some physics to do that
  public static double circlev(double rx, double ry) {
    double solarmass=1.98892e30;
    double r2=Math.sqrt(rx*rx+ry*ry);
    double numerator=(6.67e-11)*1e6*solarmass;
    return Math.sqrt(numerator/r2);
  }
  
  //Initialize N bodies
  public void startthebodies(int N) {
    double radius = 1e18;        // radius of universe
    double solarmass=1.98892e30;
    for (int i = 0; i < N; i++) {
      double px = 1e18*exp(-1.8)*(.5-Math.random());
      double py = 1e18*exp(-1.8)*(.5-Math.random());
      double magv = circlev(px,py);
      
      double absangle = Math.atan(Math.abs(py/px));
      double thetav= Math.PI/2-absangle;
      double phiv = Math.random()*Math.PI;
      double vx   = -1*Math.signum(py)*Math.cos(thetav)*magv;
      double vy   = Math.signum(px)*Math.sin(thetav)*magv;
      //Orient a random 2D circular orbit
     
           if (Math.random() <=.5) {
              vx=-vx;
              vy=-vy;
            } 
      double mass = Math.random()*solarmass*10+1e20;
      //Color a shade of blue based on mass
      int red     = (int) Math.floor(mass*254/(solarmass*10+1e20));
      int blue   = 255;
      int green    = (int) Math.floor(mass*254/(solarmass*10+1e20));
      Color color = new Color(red, green, blue);
      bodies[i]   = new Body(px, py, vx, vy, mass, color);
    }
    bodies[0]= new Body(0,0,0,0,1e6*solarmass,Color.red);//put a heavy body in the center
    
  }
  //The BH algorithm: calculate the forces
  public void addforces(int N) {
    BHTree thetree = new BHTree(q);
            // If the body is still on the screen, add it to the tree
            for (int i = 0; i < N; i++) {
                if (bodies[i].in(q)) thetree.insert(bodies[i]);
            }
            //Now, use out methods in BHTree to update the forces,
            //traveling recursively through the tree
            for (int i = 0; i < N; i++) {
                bodies[i].resetForce();
                if (bodies[i].in(q)) {
                  thetree.updateForce(bodies[i]);
                  //Calculate the new positions on a time step dt (1e11 here)
                  bodies[i].update(1e11);
                }
            }
  }
  //A function to return an exponential distribution for position
   public static double exp(double lambda) {
        return -Math.log(1 - Math.random()) / lambda;
    }
   public boolean action(Event e,Object o)
   {
     N = Integer.parseInt(t1.getText());
     if (N>10000) {
       t1.setText("10000");
       N=10000;
     }
     
       startthebodies(N);
       repaint();
     
     return true;
   }
   public class ButtonListener implements ActionListener{

    public void actionPerformed(ActionEvent evt) 
    {
        // Get label of the button clicked in event passed in
        String arg = evt.getActionCommand();    
        if (arg.equals("Restart"))
        {
          shouldrun=true;
               N = Integer.parseInt(t1.getText());
     if (N>10000) {
       t1.setText("10000");
       N=10000;
     }
     
       startthebodies(N);
       repaint();
        }
        else if (arg.equals("Stop")) 
            stop();
    }
}

}

And that's it! The objects were more complicated, but the end code is quite simple and elegant. To see both codes in action, go to the Applets page. Here, you can play around with N, the number of bodies, to watch the dynamics and get a sense for the how computationally intensive each method is.

Site Map
N-body Simulations