<!doctype html>
  <title>WASM Collision Demo</title>
  <body style="margin:0; overflow: hidden">
    <canvas id="canvas" style="width: 100vw; height: 100vh; display:block;"></canvas>
  <script type="module">
    import render from 'https://cdn.rawgit.com/guybedford/wasm-demo/dc70a37d/lib/render.js';

    function fetchAndInstantiateWasm (source, importObj) {
      return fetch(source)
      .then(response => {
        if (response.ok)
          return response.arrayBuffer();
        throw new Error(`Unable to fetch wasm ${source}.`)
      })
      .then(WebAssembly.compile)
      .then(module => {
        let instance = new WebAssembly.Instance(module, importObj);
        return instance.exports;
      });
    }

    fetchAndInstantiateWasm('https://cdn.rawgit.com/guybedford/wasm-demo/dc70a37d/lib/dynamics-coll.wasm', {
      env: {
        randomf: Math.random,
        expf: Math.exp
      }
    }).then(m => {
      const circleCount = m.getCircleCount();
      const circleDataOffset = m.getCircleDataOffset();

      const circleData = new Float32Array(m.memory.buffer, circleDataOffset, circleCount * 3);

      render(circleData, circleCount, m.init, m.timeStep);
    });
  </script>
#include <math.h>
#include <stdbool.h>

// (reduced to 1000 for init performance)
#define CIRCLE_COUNT 1000

float randomf ();
float expf (float);

struct Circle {
  float x;
  float y;
  float r;
};

struct CircleV {
  float vx;
  float vy;
};

struct Circle circleData[CIRCLE_COUNT];
struct CircleV circlevData[CIRCLE_COUNT];

int getCircleCount () {
  return CIRCLE_COUNT;
}

struct Circle* getCircleDataOffset () {
  return &circleData[0];
}

bool detectCircleCollision (float x1, float y1, float r1, float x2, float y2, float r2) {
  // before circle intersection, check bounding box intersection
  if (x1 + r1 < x2 - r2 || x1 - r1 > x2 + r2 ||
      y1 + r1 < y2 - r2 || y1 - r1 > y2 + r2)
    return false;
  // circle intersection when distance between centers < radius total
  return sqrtf(powf(x2 - x1, 2) + powf(y2 - y1, 2)) <= r1 + r2;
}

void init (float displayWidth, float displayHeight) {
  for (int i = 0; i < CIRCLE_COUNT; i++) {
    bool collision;
    float x, y, r;
    do {
      collision = false;
      x = displayWidth * randomf();
      y = displayHeight * randomf();
      // size of 0.5 - 128, exponentially distributed (maximum gl_POINT render size)
      r = ceilf(expf(randomf() * 8) / 23.2887);

      // ensure within window bounds
      if (displayWidth - (x + r) < 0 || x - r < 0 ||
          displayHeight - (y + r) < 0 || y - r < 0) {
        collision = true;
      }
      else {
        // ensure no collisions for starting position
        for (int j = 0; j < i; j++) {
          if (detectCircleCollision(x, y, r, circleData[j].x, circleData[j].y, circleData[j].r)) {
            collision = true;
            break;
          }
        }
      }
    }
    while (collision);

    circleData[i].x = x;
    circleData[i].y = y;
    circleData[i].r = r;

    // velocity of -0.1 - +0.1 pixels / iteration
    circlevData[i].vx = (randomf() - 0.5) * 0.1;
    circlevData[i].vy = (randomf() - 0.5) * 0.1;
  }
}

void timeStep (float displayWidth, float displayHeight) {
  for (int i = 0; i < CIRCLE_COUNT; i++) {
    float xi = circleData[i].x;
    float yi = circleData[i].y;
    float ri = circleData[i].r;

    float vxi = circlevData[i].vx;
    float vyi = circlevData[i].vy;

    // gravity
    vyi += 0.0001;

    // window bounds
    if (displayWidth - (xi + ri) < 0 && vxi > 0 || xi - ri < 0 && vxi < 0) {
      vxi = -vxi;
    }
    if (displayHeight - (yi + ri) < 0 && vyi > 0 || yi - ri < 0 && vyi < 0) {
      vyi = -vyi;
    }

    circleData[i].x = xi + vxi;
    circleData[i].y = yi + vyi;
    circlevData[i].vx = vxi;
    circlevData[i].vy = vyi;

    /*
     * Collision detection
     */
    for (int j = i + 1; j < CIRCLE_COUNT; j++) {
      float xj = circleData[j].x;
      float yj = circleData[j].y;
      float rj = circleData[j].r;

      if (detectCircleCollision(xi, yi, ri, xj, yj, rj)) {
        float vxj = circlevData[j].vx;
        float vyj = circlevData[j].vy;

        // calculate collision unit vector
        float collDx = xj - xi;
        float collDy = yj - yi;
        float collLen = sqrtf(collDx * collDx + collDy * collDy);
        collDx = collDx / collLen;
        collDy = collDy / collLen;

        // dot product of unit collision vector with velocity vector gives
        // 1d collision velocities before collision along collisionv vector
        float cui = collDx * vxi + collDy * vyi;
        float cuj = collDx * vxj + collDy * vyj;

        // skip collision if moving away from eachother
        if (cui - cuj <= 0)
          continue;

        // we then use 1d elastic collision equations
        // to get resultant 1d velocities after collision
        // (https://en.wikipedia.org/wiki/Elastic_collision)
        float massSum = ri + rj;
        float massDiff = ri - rj;
        float cvi = (cui * massDiff + 2 * rj * cuj) / massSum;
        float cvj = (2 * ri * cui - cuj * massDiff) / massSum;

        // calculate the collision velocity change
        float dcvi = cvi - cui;
        float dcvj = cvj - cuj;

        // apply that velocity change dotted with the collision unit vector
        // to the original velocities
        circlevData[i].vx = vxi + collDx * dcvi;
        circlevData[i].vy = vyi + collDy * dcvi;
        circlevData[j].vx = vxj + collDx * dcvj;
        circlevData[j].vy = vyj + collDy * dcvj;
      }
    }
  }
}