​

// Documentation: 

// https://wp.nyu.edu/jacky_xu/2021/02/28/simplified-wilberforce-pendulum/

​

// Equations Reference: 

// Misay A. Partnof and Steven C. Richards 

// Basic Coupled Oscillator Theory Applied to the Wilberforce Pendulum

// https://mse.redwoods.edu/darnold/math55/DEproj/sp04/stevemisay/Project1.pdf

​

// Constants Reference: 

// Berg, Ricard E., Marsahll, Todd S..

// Wilberforce Pendulum Oscillations and Normal Modes

// https://aapt.scitation.org/doi/pdf/10.1119/1.16702

​

let bob;

let anchor;

let initHeight = 200;

let boxLength = 64;

let curLength;

let boxHeight = 32;

let pixelPerMeter = 200;

​

// Constants

let k = 2.69;

let TorsionalSpringConstant = 7.86E-4;

let Inertia = 1.45E-4;

let Epsilon = 9.28E-3;

let Mass = 0.5164;

​

// here is the z0 and theta0 I picked

let z0 = 10;

let theta0 = 0;

let t = 0;

let deltaT = 0.05;

​

let wave1 = []; // visualize one side line position of the rect

let wave2 = []; // visualize the current height of the rect

​

function setup() {

  createCanvas(800, 400);

  bob = createVector(width/2, initHeight);

  anchor = createVector(width/2, 0);

  // init the omega1 and omega2 and r via equations

  w1 = pow(TorsionalSpringConstant/Inertia+Epsilon/pow(Mass*Inertia,1/2),1/2);

  w2 = pow(TorsionalSpringConstant/Inertia-Epsilon/pow(Mass*Inertia,1/2),1/2);

  r = pow(Mass/Inertia,1/2);

}

​

function draw() {

  background(112, 50, 126);

  strokeWeight(4);

  stroke(255);

  fill(45, 197, 244);

  // update the time and corresponding shape

  t += deltaT;

  bob.y = initHeight + pixelPerMeter*z(t);

  line(anchor.x, anchor.y, bob.x, bob.y);

  circle(anchor.x, anchor.y, 32);

  curLength = boxLength * cos(theta(t));

  rect(bob.x-curLength/2, bob.y, curLength, boxHeight);

  fill(197, 45, 244);

  ellipse(bob.x+curLength/2, bob.y+boxHeight, curLength/4, boxHeight/2);

  // if you prefer to see the top view, uncomment this block

  // fill(45, 197, 244);

  // ellipse(bob.x-100, height/2+boxHeight/2, boxLength, boxLength);

  // fill(197, 45, 244);

  // ellipse(bob.x-100+boxLength/2*cos(theta(t)), height/2+boxHeight/2+boxLength/2*sin(theta(t)), boxHeight/2, boxHeight/2);

​

  // data visualizations

  wave1.unshift(bob.x-curLength/2+curLength);

  wave2.unshift(bob.y);

  stroke(255, 100);

  strokeWeight(2);

  // visualize one side line position of the rect

  beginShape();

  noFill();

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

    vertex(wave1[i],i+bob.y+boxHeight);

  }

  endShape();

  // visualize the current height of the rect

  beginShape();

  noFill();

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

    vertex(i+bob.x+boxLength/2,wave2[i])

  }

  endShape();

  if (wave1.length > 1500) {

    wave1.pop();

    wave2.pop();

  }

}

​

// Equation for val z

function z(t) {

  let val = z0/2 * 1/r * (cos(w1*t)-cos(w2*t)) + theta0/2 * (cos(w1*t)+cos(w2*t));

  return val

}

​

// Equation for val theta

function theta(t) {

  let val = theta0/2 * r * (cos(w1*t)-cos(w2*t)) + z0/2 * (cos(w1*t)+cos(w2*t));

  return val

}

​