Particles moving in a vector field defined as the gradient of a two-variable function

sketch_241005a_fields.pde

/* Particles following vector fields
 * ---------------------------------
 *  
 * Copyright (C) 2025 Waldeck Schutzer (@infinitymathart) 
 * Based on code by Étienne Jacob (@etinjcb)
 *  
 * This program is free software: you can redistribute it and/or modify  
 * it under the terms of the GNU General Public License as published by  
 * the Free Software Foundation, either version 3 of the License, or  
 * (at your option) any later version.  
 *  
 * This program is distributed in the hope that it will be useful,  
 * but WITHOUT ANY WARRANTY; without even the implied warranty of  
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the  
 * GNU General Public License for more details.  
 *  
 * You should have received a copy of the GNU General Public License  
 * along with this program. If not, see <https://www.gnu.org/licenses/>.  
 */  
 
final boolean recording = true;
final boolean chromatic = false;    // Enable chromatic aberration effect
final boolean oscillation = true;  // Enable particla oscillation
final float dt = (oscillation ? 0.2 : 1) * 0.025; // Integration time step

final int duration = 10; // seconds
final float scl = 1.35;
float fac = 1;

// Requires the PostFX library for effects like chromatic aberrration (optional)
// Otherwise comment out the following 4 lines and remove references therein.
import ch.bildspur.postfx.builder.*;
import ch.bildspur.postfx.pass.*;
import ch.bildspur.postfx.*;
PostFX fx; 

///////////////////////////////////////////////////
/// various parameters to control the aesthetic
 
// This one is quite explicit
boolean use_white_rectangle = true;
// Border margin
int border = int(50*fac);
// Inverting colors or not
boolean invert_colors = false;
// Maximum point size
float maximimum_point_size = 4*fac;
 
/////////////////////////////////////////////////
/// FLOW FIELD ANIMATION ALGORITHM
 
// Number of steps
int nsteps = 2000;
// Number of particles per path
int number_of_particles_per_path = 20;
// Number of paths
int NPath = 2000;
// The total number of particles will be NPath*number_of_particles_per_path

// Number of drawings used to render each final frame with motion blur
int samplesPerFrame = 7;
// Total number of frames in the gif
int numFrames = duration*60;
// Kind of the time interval used for each frame in the motion blur
float shutterAngle = 1.5;
float bht = 1.0; // Brightness adjustment
PFont courier;

int[][] result;
float t, c;
 
float ease(float p) {
  return 3*p*p - 2*p*p*p;
}
 
float ease(float p, float g) {
  if (p < 0.5) 
    return 0.5 * pow(2*p, g);
  else
    return 1 - 0.5 * pow(2*(1 - p), g);
}
 
float mn = .5*sqrt(3), ia = atan(sqrt(.5));
 
void push() {
  pushMatrix();
  pushStyle();
}
 
void pop() {
  popStyle();
  popMatrix();
}
 
void draw() {
 
  if (!recording) {
     t = 1-map(frameCount-1, 0, numFrames, 0, 1);
     //t = mouseX*1.0/width;
    c = mouseY*1.0/height;
    if (mousePressed)
      println(c);
    pushStyle();
    pushMatrix();
    draw_();
    popMatrix();
    popStyle();
    postDraw();
  } else {
    for (int i=0; i<width*height; i++)
      for (int a=0; a<3; a++)
        result[i][a] = 0;
 
    c = 0;
    for (int sa=0; sa<samplesPerFrame; sa++) {
      t = 1-map(frameCount-1 + sa*shutterAngle/samplesPerFrame, 0, numFrames, 0, 1);
      draw_();
      loadPixels();
      for (int i=0; i<pixels.length; i++) {
        result[i][0] += pixels[i] >> 16 & 0xff;
        result[i][1] += pixels[i] >> 8 & 0xff;
        result[i][2] += pixels[i] & 0xff;
      }
    }
 
    loadPixels();
    for (int i=0; i<pixels.length; i++)
      pixels[i] = 0xff << 24 | 
        int(constrain(result[i][0]*bht/samplesPerFrame,0,255)) << 16 | 
        int(constrain(result[i][1]*bht/samplesPerFrame,0,255)) << 8 | 
        int(constrain(result[i][2]*bht/samplesPerFrame,0,255));
    updatePixels();
     
    postDraw();

    if(invert_colors){
      filter(INVERT);
    }
 
    saveFrame("/tmp/r/frame_####.png");
    println(frameCount,"/",numFrames);
     
    if (frameCount==numFrames)
      exit();
  }
}

/// A class to define paths particles take
class Path
{
  float x = random(width);
  float y = random(height);
  float s = 0.038; //random(0.03, 0.05);
  float r = 20*fac; //random(40,200); // speed
   
  ArrayList<PVector> positions = new ArrayList<PVector>();
   
  // point size
  float sz = random(fac,maximimum_point_size);
   
  // Nunmber of particles per path
  int npart = number_of_particles_per_path;
   
  // offset so that particles don't appear at the same time for each path
  float t_off = random(1);
   
  Path()
  {
    positions.add(new PVector(x,y));
  }
   
  void update()
  {
    // Enhanced Euler Integrator
    PVector k1 = field(x,y,s).mult(r);
    PVector k2 = field(x + dt*k1.x, y + dt*k1.y, s).mult(r);
    PVector res = PVector.add(k1, k2).mult(0.5);
    x += dt*res.x;
    y += dt*res.y;
    positions.add(new PVector(x,y));
  }
   
  void show()
  {
     
    strokeWeight(sz);
     
    float tt = oscillation ? ((sin(3*TAU*(t+t_off)))%1) : ((TAU*(t+t_off)))%1; // No oscillation

    int len = positions.size();
     
    for(int i=0;i<npart;i++){
      // Particle location calculated by linear interpolation from the computed positions
      float loc = constrain(map(i+tt,0,npart,0,len-1),0,len-1-0.001);
      int i1 = floor(loc);
      int i2 = i1+1;
      float interp = loc - floor(loc);
      float xx = lerp(positions.get(i1).x,positions.get(i2).x,interp);
      float yy = lerp(positions.get(i1).y,positions.get(i2).y,interp);
       
      float fact = 1;
      if(use_white_rectangle && (xx<border||xx>width-border||yy<border||yy>height-border)) fact = 0;
       
      // This is to make the particles appear and disappear gradually
      float alpha = fact*255*pow(sin(PI*(len-1-loc)/(len-1)),0.25);
       
      stroke(255,alpha);
       
      point(xx,yy);
    }
  }
}
 
Path[] array2 = new Path[NPath];
 
void path_step(){
  for(int i=0;i<NPath;i++){
    array2[i].update();
  }
}
 
//////////////////////////////////////
/// Definition of the flow field
 
// A differentiable scalar function f on the variables x and y

// This gives two cells within a bigger cell surrounded with fast moving points
//
float f0(float x, float y)
{
  x = scl*(x - width / 2);
  y = scl*(y - height / 2);
  float u = x / (50*fac);
  float v = y / (50*fac);
  float su = sin(u);
  float sv = sin(v);
  float s = (su-sv); //*(sqrt(x*x+y*y)/(50*fac))*exp(-(u*u+v*v)/(100*fac));
  return s*s*exp(-u*u-v*v)+pow((u*u+v*v),5.0)/100000000;
}

// This gives particles moving in random directions that seem to be alive
//
float f1(float x, float y)
{
  float nscl = 0.05;
  x = scl*(x - width / 2);
  y = scl*(y - height / 2);
  float u = x / (50*fac);
  float v = y / (50*fac);
  float su = sin(u);
  float sv = sin(v);
  float s = su*sv*map(noise(nscl*x,nscl*y),0,1,-0.1,1.0); //*(sqrt(x*x+y*y)/(50*fac))*exp(-(u*u+v*v)/(100*fac));
  return s;
}

// This gives (what?)
//
float f(float x, float y)
{
  x = scl*(x - width / 2);
  y = scl*(y - height / 2);
  float u = x / (50*fac);
  float v = y / (50*fac);
  float su = sin(u);
  float sv = sin(v);
  float s = su*sv*(sqrt(x*x+y*y)/(50*fac))*exp(-(u*u+v*v)/(100*fac));
  return s;
}
// The vector field is actually the gradient field of the scalar function
//
PVector field(float x, float y, float s)
{
  float h = 0.01;
  float dx = 0.5*(f(x+h,y)-f(x-h,y))/h;
  float dy = 0.5*(f(x,y+h)-f(x,y-h))/h;
  return new PVector(dy, -dx).mult(50*fac);
}

////////////////////
/// SETUP AND DRAW_
 
void settings()
{
  if (recording)
    size(2160,2160,P2D);
  else
    size(800,800,P2D);
  smooth(8);
}

void setup()
{
  /// drawing size
  fac = width/(chromatic ? 1.8*800.0 : 800.0);
  noiseSeed(1337);
  fx = new PostFX(this);  

  courier = createFont("Courier New",16*width/800.0,true);
  border = int(border*fac);
   
  /// Initialization of the array used to render frames
  result = new int[width*height][3];
   
  /// Initilization of Paths
  for(int i=0;i<NPath;i++){
    array2[i] = new Path();
  }
   
  /// Computation of Paths
  for(int i=0;i<nsteps;i++){
    // println(i+1,"/",nsteps);
    path_step();
  }
}
 
void draw_()
{
  background(0);
  //image(f_img, 0, 0);
 
  for(int i=0;i<NPath;i++)
  {
    array2[i].show();
  }
  
  // add bloom filter
  if (chromatic)
    fx.render()
      //.sobel()
      .bloom(0.5, int(20*width/800.0), 15*width/800.0)
      .chromaticAberration()
      .compose();
    else
    fx.render()
      //.sobel()
      .bloom(0.5, int(20*width/800.0), 15*width/800.0)
      // .chromaticAberration()
      .compose();
   
  if(use_white_rectangle)
  {
    fill(0);
    noStroke();
    rect(0, 0, width, border);
    rect(0, 0, border, height);
    rect(0, height-border, width, border);
    rect(width-border, 0, border, height);
    noFill();
    stroke(255);
    strokeWeight(1);
    rect(border,border,width-2*border,height-2*border);
  }

}

void postDraw()
{
  // texts 
  pushMatrix();
  translate(width/2,height/2);
  fill(255);stroke(255);
  textAlign(CENTER, CENTER);
  textFont(courier);
  text("@infinitymathart • 2025",0.0*width,0.40*height);
  popMatrix();
}