Particle Lenia simulation In Petri Dish

/*

  Based upon this Tutorial by Alexander Mordvintsev:

  Particle Lenia from scratch

  https://observablehq.com/@znah/particle-lenia-from-scratch

  The simulation shows Particle Lenia: the system of 

  particles that move to minimize local energy of their 

  interaction, leading to the formation of complex and 

  diverse structures. 

  https://arxiv.org/pdf/2305.16706

  This version by Juan Carlos Ponce Campuzano

  16/Apr/2025

  https://www.patreon.com/jcponce

  Modified to constrain particles within a circular boundary (petri dish)

*/

​

let params = {

  mu_k: 4.0, sigma_k: 1.0, w_k: 0.022,

  mu_g: 0.6, sigma_g: 0.15, c_rep: 1.0,

  dt: 0.1

};

​

let point_n = 150;

let points;

let fields;

let steps_per_frame; // Adjust this based on desired speed

let randomCol;

let dishRadius = 10; // Radius of the petri dish

​

function setup() {

  createCanvas(windowWidth, windowHeight);

  points = init(new Float32Array(point_n * 2));

  fields = {

    R_val: new Float32Array(point_n),

    U_val: new Float32Array(point_n),

    R_grad: new Float32Array(point_n * 2),

    U_grad: new Float32Array(point_n * 2)

  };

  colorMode(RGB, 200, 200, 200, 200);

  noStroke();

  steps_per_frame = int(random(3, 12));

  randomCol = random();

}

​

function draw() {

  for (let i = 0; i < steps_per_frame; ++i) step();

  background(250, 50);

  translate(width / 2, height / 2);

  scale(width / 45.0);

  // Draw petri dish outline

  //push();

  //noFill();

  push();

  stroke(10);

  fill(250, 20)

  strokeWeight(0.5);

  ellipse(0, 0, dishRadius * 2, dishRadius * 2);

  pop();

  // Draw particles

  strokeWeight(0.1);

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

    let x = points[i * 2], y = points[i * 2 + 1];

    let r = params.c_rep / (fields.R_val[i] * 5.0);

    //stroke(255)

    push();

    fill(50, 50);

    noStroke();

    ellipse(x, y, r * 2.1, r * 2.1);

    pop();

​

    push();

    stroke(10);

    strokeWeight(0.1);

    point(x, y);

    pop();

    //ellipse(x, y, r * 2, r * 2);

  }

}

​

function init(points) {

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

    // Initialize particles within the dish radius with some margin

    let angle = random(TWO_PI);

    let radius = random(dishRadius * 0.8);

    points[i * 2] = cos(angle) * radius;

    points[i * 2 + 1] = sin(angle) * radius;

  }

  return points;

}

​

function add_xy(a, i, x, y, c) {

  a[i * 2] += x * c;

  a[i * 2 + 1] += y * c;

}

​

function compute_fields() {

  const { R_val, U_val, R_grad, U_grad } = fields;

  const { c_rep, mu_k, sigma_k, w_k } = params;

  R_val.fill(repulsion_f(0.0, c_rep)[0]);

  U_val.fill(peak_f(0.0, mu_k, sigma_k, w_k)[0]);

  R_grad.fill(0);

  U_grad.fill(0);

​

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

    for (let j = i + 1; j < point_n; ++j) {

      let rx = points[i * 2] - points[j * 2];

      let ry = points[i * 2 + 1] - points[j * 2 + 1];

      const r = Math.sqrt(rx * rx + ry * ry) + 1e-20;

      rx /= r;

      ry /= r;

​

      if (r < 1.0) {

        const [R, dR] = repulsion_f(r, c_rep);

        add_xy(R_grad, i, rx, ry, dR);

        add_xy(R_grad, j, rx, ry, -dR);

        R_val[i] += R;

        R_val[j] += R;

      }

      const [K, dK] = peak_f(r, mu_k, sigma_k, w_k);

      add_xy(U_grad, i, rx, ry, dK);

      add_xy(U_grad, j, rx, ry, -dK);

      U_val[i] += K;

      U_val[j] += K;

    }

  }

}

​

function repulsion_f(x, c_rep) {

  const t = Math.max(1.0 - x, 0.0);

  return [0.5 * c_rep * t * t, -c_rep * t];

}

​

function fast_exp(x) {

  let t = 1.0 + x / 32.0;

  t *= t;

  t *= t;

  t *= t;

  t *= t;

  t *= t;

  return t;

}

​

function peak_f(x, mu, sigma, w = 1.0) {

  const t = (x - mu) / sigma;

  const y = w / fast_exp(t * t);

  return [y, -2.0 * t * y / sigma];

}

​

function step() {

  const { R_val, U_val, R_grad, U_grad } = fields;

  const { mu_g, sigma_g, dt } = params;

  compute_fields();

  let total_E = 0.0;

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

    const [G, dG] = peak_f(U_val[i], mu_g, sigma_g);

    const vx = dG * U_grad[i * 2] - R_grad[i * 2];

    const vy = dG * U_grad[i * 2 + 1] - R_grad[i * 2 + 1];

    add_xy(points, i, vx, vy, dt);

    // Boundary constraint - keep particles inside the dish

    const x = points[i * 2];

    const y = points[i * 2 + 1];

    const dist = sqrt(x * x + y * y);

    if (dist > dishRadius) {

      // Move particle back to the edge

      const angle = atan2(y, x);

      points[i * 2] = cos(angle) * dishRadius * 0.99;

      points[i * 2 + 1] = sin(angle) * dishRadius * 0.99;

    }

    total_E += R_val[i] - G;

  }

  return total_E / point_n;

}

​

function mousePressed() {

  params.mu_k = random(1.0, 5.0);

  params.sigma_k = random(0.1, 2.0);

  params.w_k = random(0.01, 0.07);

  params.mu_g = random(0.1, 1.0);

  params.sigma_g = random(0.05, 0.3);

  params.c_rep = random(0.1, 2.0);

}