Brachistochrone

function setup() {

  createCanvas(windowWidth, windowHeight);

  start = createVector(windowWidth/3, windowHeight/3);

  end = createVector(2*windowWidth/3, 2*windowHeight/3);

  endpoint_rad = 5;

  start_selected = false;

  end_selected = false;

  play = false;

  dt = 0.1;

  g = createVector(0, 1);

  ball = new Ball(p5.Vector.copy(start), createVector(0,0), 20);

  /*

  let {org, nor} = compute_ramp(start, end, start.dist(end)*1);

  ramp = new LinearRamp(org, nor);

  //*/

  //*

  //let {ctr, rad} = compute_circle_rad(start, end, start.dist(end)*1);

  let {ctr, rad} = compute_circle_vert(start, end);

  ramp = new CircularRamp(ctr, rad);

  //*/

  /*

  let {ctr, k} = compute_parabola_k(start, end, -0.001);

  ramp = new ParabolicRamp(ctr, k);

  //*/

}

​

class Ball {

  constructor(pos=createVector(0,0), vel=createVector(0,0), rad=10) {

    this.pos = pos;

    this.vel = vel;

    this.rad = rad;

  }

  draw() {

    push();

    noStroke();

    fill('grey');

    circle(this.pos.x, this.pos.y, this.rad);

    pop();

    //drawArrow(this.pos, p5.Vector.mult(this.vel,5), 'black');

  }

}

​

class LinearRamp {

  constructor(org=createVector(0,0), dir=createVector(1,0)) {

    this.org = org;

    this.dir = dir;

  }

  get_potential(pos) {

    return p5.Vector.sub(this.org, pos).dot(this.dir);

  }

  get_force(pos) {

    return this.dir;

  }

  draw() {

    let T = p5.Vector.rotate(this.dir, HALF_PI).mult(windowWidth);

    let A = p5.Vector.sub(this.org, T);

    let B = p5.Vector.add(this.org, T);

    line(A.x, A.y, B.x, B.y);

  }

}

​

​

class ParabolicRamp {

  constructor(org=createVector(0,0), k=1) {

    this.org = org;

    this.k = k;

  }

  proj(pos) {

    //*

    // ref: https://www.shadertoy.com/view/ws3GD7

    let pos_ = createVector(abs(pos.x - this.org.x), pos.y-this.org.y);

    let ik = 1.0/this.k;

    let p = ik*(pos_.y - 0.5*ik)/3.0;

    let q = 0.25*ik*ik*pos_.x;

    let h = q*q - p*p*p;

    let r = sqrt(abs(h));

​

    let x = (h>0.0) ? 

        // 1 root

        pow(q+r,1.0/3.0) + pow(abs(q-r),1.0/3.0)*Math.sign(p) :

        // 3 roots

        2.0*cos(atan2(r,q)/3.0)*sqrt(p);

    return createVector(this.org.x + x*Math.sign(pos.x - this.org.x), this.k*x*x + this.org.y);

  }

  get_potential(pos) {

    //*

    let pos_proj = this.proj(pos);

    let d = pos.dist(pos_proj);

    return (pos.y>pos_proj.y)? -d : d;

    /*/

    return this.proj(pos).y - pos.y

    //*/

  }

  get_force(pos) {

    //*

    let pos_proj = this.proj(pos);

    return p5.Vector.sub(pos, pos_proj);

    /*/

    return p5.Vector.sub(pos, this.proj(pos));

    //*/

  }

  draw() {

    push();

    for(let x=0; x<windowWidth; ++x) {

      let y = (x - this.org.x)*(x - this.org.x)*this.k + this.org.y;

      let x_ = x+1;

      let y_ = (x_ - this.org.x)*(x_ - this.org.x)*this.k + this.org.y;

      line(x, y, x_, y_);

    }

    pop();

  }

}

​

​

​

class CircularRamp {

  constructor(ctr=createVector(0,0), rad=1) {

    this.ctr = ctr;

    this.rad = rad;

  }

  get_potential(pos) {

    return this.rad - pos.dist(this.ctr);

  }

  get_force(pos) {

    return p5.Vector.sub(pos, this.ctr).normalize();

  }

  draw() {

    push();

    noFill();

    circle(this.ctr.x, this.ctr.y, this.rad*2);

    pop();

    // drawX(this.ctr); // Debug

  }

}

​

function animate() {

  /// PBD style integration

  // start with standard integration

  ball.vel.add(p5.Vector.mult(g,dt));

  let prev_pos = p5.Vector.copy(ball.pos);

  ball.pos.add(p5.Vector.mult(ball.vel,dt));

  // solve constraints

  let p = -abs(ramp.get_potential(ball.pos)); // using -abs to get a bilateral constraint

  let dp = ramp.get_force(ball.pos);

  let iter = 0, max_iter = 10;

  while(p < 0 && iter < max_iter) {

    let alpha = p / dp.magSq();

    ball.pos.add(p5.Vector.mult(dp, alpha));

    p = -abs(ramp.get_potential(ball.pos));

    dp = ramp.get_force(ball.pos);

    ++ iter;

  }

  // back-propagate position constraints onto velocity

  ball.vel = p5.Vector.sub(ball.pos, prev_pos).div(dt);

}

​

function draw() {

  // Update

  if(play) {

    animate();

  }

  // Clear

  background(255);

  // Ramp

  stroke('black');

  strokeWeight(3);

  ramp.draw();

  // EndPoints

  noStroke();

  fill('red');

  circle(start.x, start.y, endpoint_rad*2);

  fill('green');

  circle(end.x, end.y, endpoint_rad*2);

  // Ball

  ball.draw();

  return;

  /*

  // Debug

  let m = createVector(mouseX, mouseY);

  let pos_proj = ramp.proj(m);

  fill('black')

  text(ramp.get_potential(ball.pos), 10,10)

  circle(pos_proj.x, pos_proj.y, 20);

  text(ramp.get_potential(createVector(mouseX, mouseY)), 10,30)

  line(mouseX, mouseY, pos_proj.x, pos_proj.y);

  //*/

}

​

function keyPressed() {

  if (keyCode === 32) {

    play = 1-play;

  }

  if (keyCode === 82) {

    ball.pos = p5.Vector.copy(start);

    ball.vel = createVector(0,0);

  }

}

​

​

function mousePressed() {

  let m = createVector(mouseX, mouseY);

  if(start.dist(m) < endpoint_rad) {

    start_selected = true;

  }

  if(end.dist(m) < endpoint_rad) {

    end_selected = true;

  }

}

​

function mouseReleased() {

  start_selected = false;

  end_selected = false;

}

​

function mouseDragged() {

  let update = false;

  if(start_selected) {

    start.add(createVector(mouseX-pmouseX, mouseY-pmouseY));

    update = true;

  }

  else if(end_selected) {

    end.add(createVector(mouseX-pmouseX, mouseY-pmouseY));

    update = true;

  }

  if(update) {

    let {ctr, rad} = compute_circle_vert(start, end);

    ramp = new CircularRamp(ctr, rad);

  }

}