N-Body Simulation

In the spirit of creating more interesting, but ultimately useless simulations, I am working on an n-body simulation. This is just a simulation of some number of bodies floating in space, acting under the force of the gravitational pull between them. It's interesting to watch.

As with the recent Game of Life experiment, I implemented this using D with the ArcLib 2D game library. As before, this means that this should compile under any major operating system that Derelict, SDL, and OpenGL compile on. And again, I recommend compiling with DSSS – see this tutorial for help compiling ArcLib programs.

Source for main.d

This file contains the ArcLib-specific code. The N-Body–specific code is kept in another file, whose source is also included below.

I recommend changing the width and height variables, which are such that the game will scale nicely to fill my widescreen monitor without much distortion.

module main;

import arc.input;
import arc.time;
import arc.window;
import arc.font;
import arc.math.point;
import arc.math.size;
import arc.draw.color;
import arc.draw.shape;
import nbody;
import std.stdio;
import std.math;

const int width = 144, height = 90;
int frames = 0, fps = 30;

float xpos = 0, ypos = 0, zoom = 2e-10;

const Color planet_color = Color.White;

void init_game()
{
  //       xpos      ypos      xvel      yvel      rad     mass
  add_body(3.000e11, 0.000e00, 0.000e00, 0.000e00, 3.0e9, 1.989e30); // sun.gif
  add_body(5.790e10, 0.000e00, 0.000e00, 2.395e04, 1.0e9, 3.302e23); // mercury.gif
  add_body(1.082e11, 0.000e00, 0.000e00, 1.750e04, 1.7e9, 4.869e24); // venus.gif
  add_body(1.496e11, 0.000e00, 0.000e00, 1.490e04, 2.0e9, 5.974e24); // earth.gif
  add_body(2.279e11, 0.000e00, 0.000e00, 1.205e04, 1.3e9, 6.419e23); // mars.gif
}

Point real_to_screen(Point point)
{
  Point ret;

  // Translate according to viewport position
  ret.x = point.x - xpos;
  ret.y = point.y - ypos;

  // Scale by zoom level
  ret.x = real_to_screen(ret.x);
  ret.y = real_to_screen(ret.y);

  // Translate to make center of screen the origin
  ret.x += width/2;
  ret.y += height/2;

  return ret;
}

float real_to_screen(float distance)
{
  return distance * zoom;
}

void draw_bodies()
{
  for(int i = 0; i < bodies.length; i++)
  {
    Body b = bodies[i];
    drawCircle(real_to_screen(Point(b.x, b.y)),
        real_to_screen(b.radius), 30, planet_color, true);
  }
}

void draw_axes()
{
  float xmax = 0, ymax = 0;
  for(int i = 0; i < bodies.length; i++)
  {
    xmax = fmax(xmax, abs(bodies[i].x));
    ymax = fmax(ymax, abs(bodies[i].y));
  }
  drawLine(real_to_screen(Point(-xmax, 0)), real_to_screen(Point(xmax, 0)), Color.White);
  drawLine(real_to_screen(Point(0, ymax)), real_to_screen(Point(0, -ymax)), Color.White);
}

void zoom_to_fit()
{
  float xmax = 0, ymax = 0, xmin = 0, ymin = 0;
  for(int i = 0; i < bodies.length; i++ )
  {
    xmax = fmax(xmax, bodies[i].x);
    ymax = fmax(ymax, bodies[i].y);
    xmin = fmin(xmin, bodies[i].x);
    ymin = fmin(ymin, bodies[i].y);
  }

  float allxmax = fmax(abs(xmax), abs(xmin)),
        allymax = fmax(abs(ymax), abs(ymin));
  if(allxmax / allymax < width / height)
    zoom = height / allymax;
  else
    zoom = width / allxmax;
  //zoom *= 0.6;

  xpos = (xmax + xmin) / 2;
  ypos = (ymax + ymin) / 2;
}

/+ GAME LOOP +/

int main()
{
  arc.window.open("N-Body", width, height, 0 );
  arc.input.open();
  arc.time.open();

  init_game();

  while(!arc.input.keyDown(ARC_QUIT))
  {
    arc.input.process();
    arc.time.process();

    // draw
    arc.window.clear();
    zoom_to_fit();
    draw_axes();
    draw_bodies();

    // paint
    arc.window.swap();

    if(arc.input.keyReleased(ARC_RIGHT))
      fps += 3;
    if(arc.input.keyReleased(ARC_LEFT))
      fps = (fps - 3 > 1) ? fps - 3 : 1;
    if(arc.input.keyReleased(ARC_q) || arc.input.keyReleased(ARC_ESCAPE))
      break;

    // progress simulation
    do_physics(25000);
    frames++;

    limitFPS(fps);
    //if(frames >= 1) break;
  }

  arc.window.close();
  return 0;
}

Source for nbody.d

module nbody;

import std.math;
import std.random;
import std.stdio;

const float G = 6.67e-11;

struct Body
{
  float x, y, xv, yv, radius, mass;
}

struct FVector { float x = 0, y = 0; }

Body[] bodies;

void add_body(float x, float y, float xvel, float yvel, float radius, float mass)
{
  bodies.length = bodies.length + 1;
  bodies[length-1] = Body(x, y, xvel, yvel, radius, mass);
}

// Calculate physics with next body as reference
void do_physics(float dt)
{
  // calculate net forces
  FVector[] forces = new FVector[bodies.length];
  for(int a = 0; a < bodies.length; a++ )
  {
    for(int b = a + 1; b < bodies.length; b++ )
    {
      float distance = sqrt((bodies[a].x - bodies[b].x) * (bodies[a].x - bodies[b].x) + (bodies[a].y - bodies[b].y) * (bodies[a].y - bodies[b].y));
      float force = G * (bodies[a].mass * bodies[b].mass / (distance * distance));
      float d_x = bodies[a].x - bodies[b].x,
            d_y = bodies[a].y - bodies[b].y;
      float f_x = force * d_x / distance,
            f_y = force * d_y / distance;
      forces[a].x += -f_x;
      forces[a].y += -f_y;
      forces[b].x += f_x;
      forces[b].y += f_y;
    }
  }

  // calculate accelerations
  FVector[] accelerations = new FVector[bodies.length];
  for(int i = 0; i < bodies.length; i++ )
  {
    accelerations[i].x = forces[i].x / bodies[i].mass;
    accelerations[i].y = forces[i].y / bodies[i].mass;
  }

  // calculate velocity at time t+dt/2
  for(int i = 0; i < bodies.length; i++ )
  {
    bodies[i].xv += dt * accelerations[i].x;
    bodies[i].yv += dt * accelerations[i].y;
  }

  // calculate position at time t+dt
  for(int i = 0; i < bodies.length; i++ )
  {
    bodies[i].x += dt * bodies[i].xv;
    bodies[i].y += dt * bodies[i].yv;
  }
}

Posted Aug 15, 2007, in the late, late night.

Comments

Click here to view the comments on this post.