final int N = 64;
final int size = (N + 2) * (N + 2);
final int cellSize = 4;

final float fps = 30;
final float dt = 1.0 / fps;
final float diff = 0.01;
final float visc = 1;

float[] u = new float[size];
float[] u_prev = new float[size];

float[] v = new float[size];
float[] v_prev = new float[size];

float[] dens = new float[size];
float[] dens_prev = new float[size];

void setup()
{
  size(N * cellSize, N * cellSize);
  smooth();
  
  frameRate(fps);
}

void draw()
{
  background(0);
  
  getInput();
  
  velStep(u, v, u_prev, v_prev);
  densStep(dens, dens_prev, u, v);
  
  for (int y = 1; y <= N; y++)
  {
    for (int x = 1; x <= N; x++)
    {
      noStroke();
      fill(dens[IX(x, y)] * 255);
      rect((x - 1) * cellSize, (y - 1) * cellSize, cellSize, cellSize);
    }
  }
}

int IX(int i, int j)
{
  return i + (N + 2) * j;
}

void addSource(float[] x, float[] s)
{
  for (int i = 0; i < size; i++)
  {
    x[i] += dt * s[i];
  }
}

void advect(int b, float[] d, float[] d0, float[] u, float[] v)
{
  float dt0 = dt * N;
  for (int i = 1; i <= N; i++)
  {
    for (int j = 1; j <= N; j++)
    {
      float x = i - dt0 * u[IX(i, j)];
      float y = j - dt0 * v[IX(i, j)];
      
      if (x < 0.5)
        x = 0.5;
      if (x > N + 0.5)
        x = N + 0.5;
      int i0 = (int)x;
      int i1 = i0 + 1;
      
      if (y < 0.5)
        y = 0.5;
      if (y > N + 0.5)
        y = N + 0.5;
      int j0 = (int)y;
      int j1 = j0 + 1;
      
      float s1 = x - i0;
      float s0 = 1 - s1;
      
      float t1 = y - j0;
      float t0 = 1 - t1;
      
      d[IX(i, j)] = s0 * (t0 * d0[IX(i0, j0)] + t1 * d0[IX(i0, j1)]) + s1 * (t0 * d0[IX(i1, j0)] + t1 * d0[IX(i1, j1)]);
    }
  }
  
  setBnd(b, d);
}

void densStep(float[] x, float[] x0, float[] u, float[] v)
{
  addSource(x, x0);
  System.arraycopy(x, 0, x0, 0, size);
  diffuse(0, x0, x, diff);
  advect(0, x, x0, u, v);
}

void diffuse(int b, float[] x, float[] x0, float diff)
{
  float a = dt * diff * N * N;
  
  for (int k = 0; k < 20; k++)
  {
    for (int i = 1; i <= N; i++)
    {
      for (int j = 1; j <= N; j++)
      {
        x[IX(i, j)] = (x0[IX(i, j)] + a * (x[IX(i - 1, j)] + x[IX(i + 1, j)] + x[IX(i, j - 1)] + x[IX(i, j + 1)])) / (1 + 4 * a);
      }
    }
  }
  
  setBnd(b, x);
}

void getInput()
{
  int mx = mouseX / cellSize + 1;
  int my = mouseY / cellSize + 1;
  
  for (int x = 0; x < N + 2; x++)
  {
    dens_prev[IX(x, 0)] = 0;
    u_prev[IX(x, 0)] = 0;
    v_prev[IX(x, 0)] = 0;
    
    dens_prev[IX(x, N + 1)] = 0;
    u_prev[IX(x, N + 1)] = 0;
    v_prev[IX(x, N + 1)] = 0;
  }
  
  for (int y = 0; y < N + 2; y++)
  {
    dens_prev[IX(0, y)] = 0;
    u_prev[IX(0, y)] = 0;
    v_prev[IX(0, y)] = 0;
    
    dens_prev[IX(N + 1, y)] = 0;
    u_prev[IX(N + 1, y)] = 0;
    v_prev[IX(N + 1, y)] = 0;
  }
  
  for (int y = 1; y <= N; y++)
  {
    for (int x = 1; x <= N; x++)
    {
      float d;
      if (mousePressed)
      {
        float r = sqrt(pow(mx - x, 2) + pow(my - y, 2));
        
        if (r > 2)
          d = 0;
        else if (r == 0)
          d = 1;
        else
          d = 1.0 / pow(r, 2);
      }
      else
        d = 0;
      
      dens_prev[IX(x, y)] = d / dt;
    }
  }
  
  for (int y = 1; y <= N; y++)
  {
    for (int x = 1; x <= N; x++)
    {
      float dx;
      float dy;
      
      if (mousePressed)
      {
        float scale;
        
        float r = sqrt(pow(mx - x, 2) + pow(my - y, 2));
        if (r == 0)
        {
          scale = 4.0;
        }
        else
          scale = 4.0 / r;
        
        dx = (mouseX - pmouseX) * scale;
        dy = (mouseY - pmouseY) * scale;
      }
      else
      {
        dx = 0;
        dy = 0;
      }
      
      u_prev[IX(x, y)] = dx;
      v_prev[IX(x, y)] = dy;
    }
  }
}

void project(float[] u, float[] v, float[] p, float[] div)
{
  float h = 1.0 / N;
  
  for (int i = 1; i <= N; i++)
  {
    for (int j = 1; j <= N; j++)
    {
      div[IX(i, j)] = -0.5 * h * (u[IX(i + 1, j)] - u[IX(i - 1, j)] + v[IX(i, j + 1)] - v[IX(i, j - 1)]);
      p[IX(i, j)] = 0;
    }
  }
  setBnd(0, div);
  setBnd(0, p);
  
  for (int k = 0; k < 20; k++)
  {
    for (int i = 1; i <= N; i++)
    {
      for (int j = 1; j <= N; j++)
      {
        p[IX(i, j)] = (div[IX(i, j)] + p[IX(i - 1, j)] + p[IX(i + 1, j)] + p[IX(i, j - 1)] + p[IX(i, j + 1)]) / 4;
      }
    }
    
    setBnd(0, p);
  }
  
  for (int i = 1; i <= N; i++)
  {
    for (int j = 1; j <= N; j++)
    {
      u[IX(i, j)] -= 0.5 * (p[IX(i + 1, j)] - p[IX(i - 1, j)]) / h;
      v[IX(i, j)] -= 0.5 * (p[IX(i, j + 1)] - p[IX(i, j - 1)]) / h;
    }
  }
  setBnd(1, u);
  setBnd(2, v);
}

void setBnd(int b, float[] x)
{
  for (int i = 1; i <= N; i++)
  {
    x[IX(0, i)] = (b == 1) ? -x[IX(1, i)] : x[IX(1, i)];
    x[IX(N + 1, i)] = (b == 1) ? -x[IX(N, i)] : x[IX(N, i)];
    x[IX(i, 0)] = (b == 2) ? -x[IX(i, 1)] : x[IX(i, 1)];
    x[IX(i, N + 1)] = (b == 2) ? -x[IX(i, N)] : x[IX(i, N)];
  }
  
  x[IX(0, 0)] = 0.5 * (x[IX(1, 0)] + x[IX(0, 1)]);
  x[IX(0, N + 1)] = 0.5 * (x[IX(1, N + 1)] + x[IX(0, N)]);
  x[IX(N + 1, 0)] = 0.5 * (x[IX(N, 0)] + x[IX(N + 1, 1)]);
  x[IX(N + 1, N + 1)] = 0.5 * (x[IX(N, N + 1)] + x[IX(N + 1, N)]);
}

void velStep(float[] u, float[] v, float[] u0, float[] v0)
{
  addSource(u, u0);
  addSource(v, v0);
  
  System.arraycopy(u, 0, u0, 0, size);
  System.arraycopy(v, 0, v0, 0, size);
  
  diffuse(1, u0, u, visc);
  diffuse(2, v0, v, visc);
  project(u0, v0, u, v);
  
  advect(1, u, u0, u0, v0);
  advect(2, v, v0, u0, v0);
  project(u, v, u0, v0);
}