Gaussian Free Field

gff.pde

int gffSize = 30; // field density
int rays = 150;
float rsteps = TWO_PI / rays;
float[][] gff;

void setup() {
  size(800,800);
  smooth(8);
  noStroke();
  fill(235);
  background(21,20,19);
  colorMode(HSB,360,100,100,100);
  
  // initialize field
  gff = new float[gffSize][gffSize];
  for(int x=0;x<gffSize;x++) {
    for(int y=0;y<gffSize;y++) {
      gff[x][y] = randomGaussian();
    }
  }
}

void draw() {
  
  fill(0,0,10,200);
  rect(0,0,width,height);
  
  float kappa = frameCount * 0.0031;
  float hf = sqrt(8.0 * (kappa / PI)) / (4.0 - kappa);
  for(int t=0; t<rays; t++) {
    float theta = t * rsteps;
    float x = width * 0.5;
    float y = height * 0.5;
    
    fill(degrees(theta),90,100,95);
    for(int i=0; i<350;i++) {
      float v = hf * getFieldValue(x,y);
      // uncomment below to use Perlin noise instead of Gaussian field
      //float v = hf * 2.0 * (noise(x/(float)gffSize,y/(float)gffSize)-0.5);
      float sx = sin(v + theta);
      float sy = cos(v + theta);
      
      rect(x,y,0.8,0.8);
      
      x+=sx;
      y+=sy;
    }
  }
}

// bilinear interpolation
float getFieldValue(float x, float y) {
  float px = map(x,0,width,0,gffSize);
  float py = map(y,0,height,0,gffSize);
  int ix = floor(px);
  int iy = floor(py);
  float rx = px - ix;
  float ry = py - iy;
  
  ix = ix % gffSize;
  iy = iy % gffSize;
  int ixp = (ix+1) % gffSize;
  int iyp = (iy+1) % gffSize;
  
  float vy1 = lerp(gff[ix][iy], gff[ixp][iy], rx);
  float vy2 = lerp(gff[ix][iyp], gff[ixp][iyp], rx);
  
  return lerp(vy1,vy2,ry);
}

void keyPressed() {
  if(keyCode==32) {
    saveFrame("######.jpg");
  }
}