Orbit simulator on an Arduino Uno

// Orbit simulator: a planet and its moon, pulling on each other on a 128x64
// OLED, with a knob that turns gravity up and down.
//
//   OLED SCL -> A5   SDA -> A4   VCC -> 5 V   GND -> GND
//   potentiometer: one end 5 V, the other GND, wiper -> A0   (gravity)
//   LAUNCH button -> pin 2 and GND: put the moon back in a round orbit
//
// # The one law that does all the work
//
// Newton found that any two masses pull on each other along the line between
// them, and that the pull gets weaker with the square of the distance:
//
//   force = G * m1 * m2 / r^2
//
// G is how strong gravity is. The knob changes G. Everything else here is
// just working out what that pull does to the two bodies, over and over.
//
// # How the sketch moves them
//
// Time goes by in tiny steps of 1/60 of a second. For each step:
//
//   1. Work out the pull on each body from the distance between them.
//   2. Change each body's velocity by pull * step (a push makes it speed up).
//   3. Change each body's position by velocity * step.
//
// That is called Euler's method. It is the simplest way there is, and doing
// step 2 before step 3 (rather than the other way round) keeps the orbit from
// slowly spiralling out, because it saves energy well.
//
// The planet is nine times heavier than the moon, so the moon does most of the
// moving. But the planet is pulled too, and wobbles in a tiny circle of its
// own: both bodies really circle their shared balance point. Watch the planet
// closely and you will see it.
//
// # What to try
//
// With the knob in the middle, G is 1.00 and the moon runs round in a neat
// circle. Press LAUNCH after turning the knob and it starts a fresh round
// orbit at the new gravity. But if you turn the knob WHILE it is flying, the
// moon keeps its speed and meets a different pull:
//
//   stronger gravity   the moon is yanked inward and swings in a long oval,
//                      and if you push it too far it crashes into the planet
//   weaker gravity     the moon is too fast for the pull and drifts outward
//                      in a wider oval; below about half, it escapes for good
//
// Every time the moon completes a lap the ORB counter goes up.
//
// # Numbers on this board
//
// The Uno has no hardware for fractions, so the maths library does them in
// software: a handful of float calculations a step is easy, thousands would
// not be. Distances are in screen pixels and time in seconds, so a speed is
// pixels per second.

#include "mokxi_i2c.h"
#include "mokxi_ssd1306.h"
#include <math.h>

const uint8_t SCL_PIN = A5;
const uint8_t SDA_PIN = A4;
const uint8_t KNOB_PIN = A0;
const uint8_t LAUNCH_PIN = 2;

SoftI2c bus(SCL_PIN, SDA_PIN, 400);
Ssd1306 oled(bus, SSD1306_ADDRESS);

// The middle of the screen, which is also the balance point of the pair.
const float CX = 64.0f;
const float CY = 34.0f;

const float M_PLANET = 9.0f;
const float M_MOON = 1.0f;
/** G when the knob is at 1.00, picked so a round orbit at 24 pixels takes six seconds. */
const float G0 = 1500.0f;
/** Starting distance between the two bodies. */
const float R0 = 24.0f;
/** One step of simulated time, and how many steps go into each frame drawn. */
const float DT = 1.0f / 60.0f;
const uint8_t STEPS = 2;

float px, py, pvx, pvy; // the planet
float mx, my, mvx, mvy; // the moon
float gain = 1.0f;      // the knob: G / G0
uint16_t orbits = 0;
bool prevBelow = false;

// The trail: the last 36 places the moon was, drawn as dots.
const uint8_t TRAIL = 36;
int8_t trailX[TRAIL];
int8_t trailY[TRAIL];
uint8_t trailAt = 0;
uint8_t trailCount = 0;
uint8_t frameNo = 0;

int16_t shownMx = -999, shownMy = -999, shownPx = -999, shownPy = -999;
int shownGain = -1;
int shownOrbits = -1;
uint32_t lastFrame = 0;
uint32_t lastReport = 0;
bool wasDown = false;

/** The knob, turned into a gain between 0.25 and 3.25 (1.00 in the middle). */
static float readGain() {
  float x = (float)analogRead(KNOB_PIN) / 1023.0f;
  return 0.25f + 3.0f * x * x;
}

static void drawBodies(bool on) {
  if (shownPx > -900) {
    oled.fillCircle(shownPx, shownPy, 4, on);
    oled.fillCircle(shownMx, shownMy, 2, on);
  }
}

/** Put the moon into a round orbit for the gravity the knob is set to now. */
static void launch() {
  gain = readGain();
  float g = G0 * gain;
  // Speed of one body past the other for a circle: v^2 = G * (M + m) / r.
  float v = sqrtf(g * (M_PLANET + M_MOON) / R0);
  float total = M_PLANET + M_MOON;
  // Put the balance point at the middle of the screen and give the pair no
  // overall drift: the lighter one is further out and moves faster.
  mx = CX + R0 * M_PLANET / total;
  my = CY;
  px = CX - R0 * M_MOON / total;
  py = CY;
  mvx = 0.0f;
  mvy = v * M_PLANET / total;
  pvx = 0.0f;
  pvy = -v * M_MOON / total;
  orbits = 0;
  prevBelow = true;
  oled.clear();
  trailCount = 0;
  trailAt = 0;
  shownPx = shownPy = shownMx = shownMy = -999;
  shownGain = -1;
  shownOrbits = -1;
  // A few stars for a sky.
  uint16_t seed = 12345u;
  for (uint8_t i = 0; i < 16; i++) {
    seed = (uint16_t)(seed * 25173u + 13849u);
    int16_t sx = (int16_t)((seed >> 4) % 128);
    seed = (uint16_t)(seed * 25173u + 13849u);
    int16_t sy = (int16_t)(10 + (seed >> 4) % 54);
    oled.setPixel(sx, sy, true);
  }
  Serial.println("launch");
}

static void showMessage(const char *text) {
  oled.fillRect(20, 26, 88, 14, false);
  oled.rect(20, 26, 88, 14, true);
  int16_t w = 0;
  for (const char *p = text; *p != '\0'; p++) w += 6;
  oled.text((int16_t)((128 - w) / 2), 30, text, 1, true, false);
  oled.display();
  delay(1200);
}

/** Move both bodies forward by one step: the pull, then the velocities, then the places. */
static void step() {
  float dx = mx - px;
  float dy = my - py;
  // A little softening (0.5) stops the pull blowing up if the bodies touch.
  float r2 = dx * dx + dy * dy + 0.5f;
  float r = sqrtf(r2);
  float inv3 = 1.0f / (r2 * r);
  float g = G0 * gain;
  // The pull on each body is G times the OTHER body's mass over r squared,
  // along the line between them. Dividing by r once more turns (dx, dy) into
  // a direction, which is where the 1/r^3 comes from.
  float aMoon = -g * M_PLANET * inv3;
  float aPlanet = g * M_MOON * inv3;
  mvx += aMoon * dx * DT;
  mvy += aMoon * dy * DT;
  pvx += aPlanet * dx * DT;
  pvy += aPlanet * dy * DT;
  mx += mvx * DT;
  my += mvy * DT;
  px += pvx * DT;
  py += pvy * DT;
  // One lap is the moon crossing the line to the right of the planet,
  // going downwards on the screen (increasing y).
  bool below = (my - py) >= 0.0f;
  if (below && !prevBelow && (mx - px) > 0.0f) orbits++;
  prevBelow = below;
}

static void drawHud() {
  int g100 = (int)(gain * 100.0f + 0.5f);
  if (g100 != shownGain) {
    shownGain = g100;
    oled.fillRect(0, 0, 50, 8, false);
    oled.text(0, 0, "G", 1, true, false);
    int16_t x = 8;
    x = oled.number(x, 0, g100 / 100, 1, true, false);
    x = oled.text(x, 0, ".", 1, true, false);
    if (g100 % 100 < 10) x = oled.text(x, 0, "0", 1, true, false);
    oled.number(x, 0, g100 % 100, 1, true, false);
  }
  if ((int)orbits != shownOrbits) {
    shownOrbits = orbits;
    oled.fillRect(84, 0, 44, 8, false);
    oled.text(84, 0, "ORB", 1, true, false);
    oled.number(108, 0, orbits, 1, true, false);
  }
}

static void frame() {
  gain = readGain();
  for (uint8_t i = 0; i < STEPS; i++) step();

  float dx = mx - px;
  float dy = my - py;
  float r = sqrtf(dx * dx + dy * dy);
  float away = sqrtf((mx - CX) * (mx - CX) + (my - CY) * (my - CY));
  if (r < 6.0f) {
    drawBodies(false);
    showMessage("IMPACT");
    Serial.println("impact");
    launch();
    return;
  }
  if (away > 110.0f) {
    showMessage("ESCAPED");
    Serial.println("escaped");
    launch();
    return;
  }

  // Rub out the old bodies, add a dot to the trail, draw the bodies where they are now.
  drawBodies(false);
  frameNo++;
  if ((frameNo & 3u) == 0u) {
    if (trailCount == TRAIL) {
      oled.setPixel(trailX[trailAt], trailY[trailAt], false);
    } else {
      trailCount++;
    }
    int16_t tx = (int16_t)(mx + 0.5f);
    int16_t ty = (int16_t)(my + 0.5f);
    trailX[trailAt] = (int8_t)(tx > 127 ? 127 : (tx < 0 ? 0 : tx));
    trailY[trailAt] = (int8_t)(ty > 63 ? 63 : (ty < 0 ? 0 : ty));
    trailAt = (uint8_t)((trailAt + 1) % TRAIL);
    oled.setPixel(tx, ty, true);
  }
  shownPx = (int16_t)(px + 0.5f);
  shownPy = (int16_t)(py + 0.5f);
  shownMx = (int16_t)(mx + 0.5f);
  shownMy = (int16_t)(my + 0.5f);
  drawBodies(true);
  drawHud();
  oled.display();
}

void setup() {
  Serial.begin(115200);
  Serial.println("Mokxi Uno: orbit simulator");
  pinMode(LAUNCH_PIN, INPUT_PULLUP);
  oled.begin();
  launch();
  lastFrame = millis();
}

void loop() {
  bool down = digitalRead(LAUNCH_PIN) == LOW;
  if (down && !wasDown) launch();
  wasDown = down;

  uint32_t now = millis();
  if ((uint32_t)(now - lastFrame) < 33u) {
    delay(1);
    return;
  }
  lastFrame += 33u;
  frame();

  if ((uint32_t)(now - lastReport) >= 500u) {
    lastReport = now;
    float dx = mx - px;
    float dy = my - py;
    Serial.print("g ");
    Serial.print((int)(gain * 100.0f + 0.5f));
    Serial.print(" r ");
    Serial.print((int)(sqrtf(dx * dx + dy * dy) + 0.5f));
    Serial.print(" orbits ");
    Serial.println(orbits);
  }
}