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