// Therm
// CA based equilibrium

float[][] land, future;
float speed = 1;
float decay = .01;

int n = 60;

void setup()
{
  size(300, 300);
  background(0);
  noStroke();
  fill(210, 198, 61);
  land = new float[n][n];
  future = new float[n][n];
  for (int y = 0; y < n; y++)
    for (int x = 0; x < n; x++)
      land[x][y] = 0;
  rectMode(CENTER_DIAMETER);
}

void loop()
{  
  // Calculate velocity towards next cell states
  for (int y = 0; y < n; y++)
    for (int x = 0; x < n; x++)
    {
      // Average orthogonal neighbor 
      float avg = 0;
      avg += land[(x + 1) % n][y];
      avg += land[(x + n - 1) % n][y];
      avg += land[x][(y + 1) % n];
      avg += land[x][(y + n - 1) % n];
      avg = avg / 4;
      float diff = avg - land[x][y];
      
      // Set offset to apply to state
      future[x][y] = diff * speed - decay;
    }

  // User input
  future[max(min(mouseX / (width / n), n-1),0)][max(min(mouseY / (width / n),n-1),0)] = 20;  

  // Update cell states
  for (int y = 0; y < n; y++)
    for (int x = 0; x < n; x++)
    {
      land[x][y] += future[x][y];
      land[x][y] = max(min(land[x][y], width / n), 0);
    }

  // Draw cells
  for (int y = 0; y < n; y++)
    for (int x = 0; x < n; x++)
      rect(x * width/n + width/n/2, y * width/n + width/n/2, land[x][y], land[x][y]);
}