// Inspired by advanced character physics paper

// https://www.cs.cmu.edu/afs/cs/academic/class/15462-s13/www/lec_slides/Jakobsen.pdf

​

let particles = [];

let temp_particles = []

let old_particles = [];

let forces = []

let locks = []

let w, w2;

const G = 3

const N = 20

const dt = 1;

const max_iters = 10;

let slider;

​

function setup() {

  slider = createSlider(0, 10, 0);

  slider.position(10, 10);

  createCanvas(400, 400);

  w = (0.5*width)/N;

  w2 = w*w;

  let offset = 0.25 * width

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

    for(let j = 0; j < N; ++j) {

      let pos = [i*w + offset, j*w*0.8+ 10]

      particles.push(pos)

      temp_particles.push(pos)

      old_particles.push(pos)

      forces.push([random()*5 - 2.5, 0])

    }

  }

  // Lock top row

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

    locks[i*N] = true;

  }

}

​

function draw() {

  background(220);

  strokeWeight(0.5)

  // Show all particles

  for(let i = 0; i < N-1; ++i) {

      beginShape(TRIANGLE_STRIP);

      for(let j = 0; j < N; ++j) {

        let p_1 = particles[i*N + j];

        let p_2 = particles[(i+1)*N + j];

        vertex(p_1[0], p_1[1]);

        vertex(p_2[0], p_2[1]);

      }

      endShape();

   }

  // Compute Forces

  for(let i = 0; i < N*N; ++i) {

    let r = particles[i]

    forces[i][0] = randomGaussian()*slider.value()

    forces[i][1] = G + randomGaussian()*slider.value()

  }

  // Verlet integration

  // 1. Copy current positions in temp buffer

  for(let i = 0; i < N*N; ++i) {

    temp_particles[i] = particles[i];

  }

  // 2. Compute new positions

  for(let i = 0; i < N*N; ++i) {

    if(locks[i]) continue;

    particles[i][0] = 2*particles[i][0] - old_particles[i][0] + forces[i][0]*dt*dt;

    particles[i][1] = 2*particles[i][1] - old_particles[i][1] + forces[i][1]*dt*dt; 

  }

  // 3. Store old positions

  for(let i = 0; i < N*N; ++i) {

    old_particles[i] = temp_particles[i];

  }

  // Circle

  let c_x = mouseX;

  let c_y = mouseY;

  let R = 30;

  let R2 = R*R;

  circle(c_x, c_y, 2*R);

  // Satisfy constraints using single relaxed jacobi iteration

  for(let iter = 0; iter < max_iters; ++iter) {

    // 1. Boundary conditions

    for(let i = 0; i < N*N; ++i) {

      // Box

      particles[i][0] = max(particles[i][0], 0)

      particles[i][0] = min(particles[i][0], width)

      particles[i][1] = max(particles[i][1], 0)

      particles[i][1] = min(particles[i][1], height)

      // Circle

      let r = [particles[i][0] - c_x, particles[i][1] - c_y]

      let r2 = r[0]*r[0] + r[1]*r[1]

      if(r2 <= R2) {

        particles[i][0] = r[0] / sqrt(r2) * R + c_x;

        particles[i][1] = r[1] / sqrt(r2) * R + c_y;

      }

    }

​

    // 2. Spring constraints

    // 2.1 Horizontal constraints

    for(let i = 1; i < N; ++i) {

      for(let j = 0; j < N; ++j) {

        // horizontal constraint

        let p1 = particles[(i-1)*N + j];

        let p2 = particles[i*N + j];

        // delta = x2 - x1

        let delta = [p2[0] - p1[0], p2[1] - p1[1]];

        // if(delta[0] != 0) print(delta)

        // delta2 = dot(delta, delta)

        let delta2 = delta[0]*delta[0] + delta[1]*delta[1];

        let s = w2 / (delta2 + w2) - 0.5;

        delta[0] *= s;

        delta[1] *= s;

        if(!locks[(i-1)*N + j]){

          p1[0] -= delta[0];

          p1[1] -= delta[1];

        }

        if(!locks[i*N + j]){

          p2[0] += delta[0];

          p2[1] += delta[1];  

        }

      }

    }

​

    // 2.2 Vertical constraint

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

      for(let j = 1; j < N; ++j) {

        // horizontal constraint

        let p1 = particles[i*N + (j-1)];

        let p2 = particles[i*N + j];

        // delta = x2 - x1

        let delta = [p2[0] - p1[0], p2[1] - p1[1]];

        // if(delta[0] != 0) print(delta)

        // delta2 = dot(delta, delta)

        let delta2 = delta[0]*delta[0] + delta[1]*delta[1];

        let s = w2 / (delta2 + w2) - 0.5;

        delta[0] *= s;

        delta[1] *= s;

        if(!locks[i*N + (j-1)]){

          p1[0] -= delta[0];

          p1[1] -= delta[1];

        }

        if(!locks[i*N + j]){

          p2[0] += delta[0];

          p2[1] += delta[1];  

        }

      }

    }

  }

  // 

  // 3. Lock constraints

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

    for(let j = 0; j < N; ++j) {

      if(locks[i*N + j]) {

        particles[i*N + j] = old_particles[i*N + j];

      }

    }

  }

}