w: 500,

  h: 500

};

var pivot = {

  x: cnv.w / 2,

  y: cnv.h / 2

};

var fr = 60;

var loops = 60;

​

var tracer;

​

// The parameters of the problem

const g = 9.8;

const l1 = 1;

const l2 = 1;

const m1 = 1;

const m2 = 1;

// First pendulum

const a11_i = 1.57;

const a12_i = -1.5;

const p11_i = 0;

const p12_i = 0;

// Second pendulum

const a21_i = 1.58;

const a22_i = -1.5;

const p21_i = 0;

const p22_i = 0;

​

const dt = 1.0 / fr / loops;

const Lscale = 0.4 * cnv.h / (l1 + l2);

const L1 = Lscale * l1;

const L2 = Lscale * l2;

​

// The variables of the problem 

var t, a1, a2;

var x1p, y1p, x2p, y2p;

​

function setup() {

  // put setup code here

  createCanvas(cnv.w, cnv.h);

  frameRate(fr);

  tracer = createGraphics(cnv.w, cnv.h);

​

  // The initial conditions of the problem

  t = 0;

  a1 = [a11_i, a12_i, p11_i, p12_i];

  a2 = [a21_i, a22_i, p21_i, p22_i];

​

  x1p = pivot.x + L1 * sin(a11_i) + L2 * sin(a12_i);

  y1p = pivot.y + L1 * cos(a11_i) + L2 * cos(a12_i);

  x2p = pivot.x + L1 * sin(a21_i) + L2 * sin(a22_i);

  y2p = pivot.y + L1 * cos(a21_i) + L2 * cos(a22_i);

}

​

function draw() {

  for (let i = 0; i < loops; i++) {

    background(220);

    text("Scale: 1 m is "+Lscale+" px", 5, 15);

​

    let a11 = a1[0];

    let a12 = a1[1];

    let a21 = a2[0];

    let a22 = a2[1];

​

    // Convert solutions to Cartesian

    x11 = pivot.x + L1 * sin(a11);

    y11 = pivot.y + L1 * cos(a11);

    x12 = x11 + L2 * sin(a12);

    y12 = y11 + L2 * cos(a12);

    x21 = pivot.x + L1 * sin(a21);

    y21 = pivot.y + L1 * cos(a21);

    x22 = x21 + L2 * sin(a22);

    y22 = y21 + L2 * cos(a22);

​

    // // tracer

    // let w1 = 1.9 * exp(-(a1[3] ** 2) / 10) + 0.7

    // let w2 = 1.9 * exp(-(a2[3] ** 2) / 10) + 0.7

    // tracer.strokeCap(SQUARE);

    // tracer.strokeWeight(w1);

    // tracer.stroke(150, 0, 0, 100);

    // tracer.line(x1p, y1p, x12, y12);

    // tracer.strokeWeight(w2);

    // tracer.stroke(0, 0, 150, 100);

    // tracer.line(x2p, y2p, x22, y22);

    // // if((frameCount*loops)%8 == 0) tracer.background(220, 1);

    // image(tracer, 0, 0);

​

    // The pivot and the pendulum

    paint(x11, y11, x12, y12);

    paint(x21, y21, x22, y22);

​

    // step the solver

    a1 = step(t, a1);

    a2 = step(t, a2);

    t += dt;

​

    // x1p = x12;

    // y1p = y12;

    // x2p = x22;

    // y2p = y22;

  }

}

​

function evolve(t, a) {

  let a1 = a[0];

  let a2 = a[1];

  let p1 = a[2];

  let p2 = a[3];

​

  let da = a1 - a2;

  let h1 = (p1 * p2 * sin(da)) / (l1 * l2 * (m1 + m2 * (sin(da)) ** 2));

  let h2 = (m2 * l2 * l2 * p1 * p1 + (m1 + m2) * l1 * l1 * p2 * p2 - 2 * m2 * l1 * l2 * p1 * p2 * cos(da)) / (2 * l1 * l1 * l2 * l2 * (m1 + m2 * (sin(da)) ** 2) ** 2);

​

  let da1 = (l2 * p1 - l1 * p2 * cos(da)) / (l1 * l1 * l2 * (m1 + m2 * (sin(da)) ** 2));

  let da2 = (-m2 * l2 * p1 * cos(da) + (m1 + m2) * l1 * p2) / (m2 * l1 * l2 * l2 * (m1 + m2 * (sin(da)) ** 2));

  let dp1 = -((m1 + m2) * g * l1 * sin(a1) + h1 - h2 * sin(2 * da));

  let dp2 = -(m2 * g * l2 * sin(a2) - h1 + h2 * sin(2 * da));

​

  return [da1, da2, dp1, dp2];

}

​

function step(t, a) {

  // Euler method

  let a_ = [];

  let da = evolve(t, a);

  for(let i = 0; i < a.length; i++) {

    a_.push(a[i] + dt*da[i]);

  }

​

  return a_;

}

​

function paint(x1, y1, x2, y2) {

  fill(5);

  ellipse(pivot.x, pivot.y, 6, 6);

  stroke(20);

  fill(5);

  line(pivot.x, pivot.y, x1, y1);

  ellipse(x1, y1, 16);

  fill(80);

  line(x1, y1, x2, y2);

  ellipse(x2, y2, 16);

}