Movers with Mutual Gravity

At the end of class, on Thursday, Feb. 10, we almost had the code below running well.

It illustrates mutual gravity among movers. I had to make four tweaks to the code to make the simulation look nice. Search the code for the word “TWEAK” and you will find the four tweaks with comments explaining what I did.

Implementation

// This code illustrates *The Nature of Code* Sections 2.9 and 2.10.

class Mover {
  PVector location;
  PVector velocity;
  PVector acceleration;
  float mass;

  Mover(float m, float x, float y, float vx, float vy) {
    mass = m;
    location = new PVector(x, y);
    velocity = new PVector(vx, vy);
    acceleration = new PVector(0, 0);
  }

  void applyForce(PVector force) {
    PVector f = PVector.div(force, mass);
    acceleration.add(f);
  }

  void update() {
    velocity.add(acceleration);
    location.add(velocity);
    acceleration.mult(0);
  }

  void display() {
    stroke(0);
    fill(175);
    ellipseMode(CENTER);
    ellipse(location.x, location.y, mass*16, mass*16);
  }

  void checkEdges() {
    if (location.x > width) {
      location.x = width;
      velocity.x *= -1;
    } else if (location.x < 0) {
      velocity.x *= -1;
      location.x = 0;
    }
    if (location.y > height) {
      location.y = height;
      velocity.y *= -1;
    } else if (location.y < 0) {
      location.y = 0;
      velocity.y *= -1;
    }
  }
}

Mover[] movers = new Mover[4];

void setup() {
  size(1000, 600, P2D);
  for (int i = 0; i < movers.length; ++i) {
    // make our movers randomly
    float mass = random(1, 5);
    float x = random(0, width);
    float y = random(0, height);
    // TWEAK 1: Calmed the initial velocities down a little.
    float vx = random(-3, 3);
    float vy = random(-3, 3);
    movers[i] = new Mover(mass, x, y, vx, vy);
  }
}

void draw() {
  background(255);
  for (int i = 0; i < movers.length; ++i) {
    for (int j = 0; j < movers.length; ++j) {
      // force of mover j on mover i
      if (j != i) {
        PVector gravity_of_j_on_i = PVector.sub(movers[j].location, movers[i].location);
        float separation = gravity_of_j_on_i.mag();
        // TWEAK 2: We were getting close to division by zero when movers got very close.
        // Make sure the denominator doesn't get too small (causing a huge acceleration).
        if (separation < 10) {
          separation = 10;
        }
        gravity_of_j_on_i.mult(200 * movers[i].mass * movers[j].mass / pow(separation, 3));
        movers[i].applyForce(gravity_of_j_on_i);
      }
    }
  }
  // TWEAK 3: Added a wee bit of friction to further calm things down. It's not realistic. Space
  // has no friction. But it makes the simulation nicer.
  for (int i = 0; i < movers.length; ++i) {
    PVector friction = movers[i].velocity.copy();
    friction.normalize();
    friction.mult(-0.012);
    movers[i].applyForce(friction);
  }
  // TWEAK 4: Only update, display and check edges after ALL the forces have been applied.
  // Otherwise can have violations of Newton's 3rd Law.
  for (int i = 0; i < movers.length; ++i) {
    movers[i].update();
    movers[i].display();
    movers[i].checkEdges();
  }
}