Processing: Reaction Diffusion

int N = 150;
 
//System parameters
float diffU;
float diffV;
float paramF;
float paramK;
 
boolean rndInitCondition = false;
 
float[][] U = new float[N][N];
float[][] V = new float[N][N];
float[][] cH = new float[N][N];
float[][] cB = new float[N][N];
 
float[][] dU = new float[N][N];
float[][] dV = new float[N][N];
 
int[][] offset = new int[N][2];

int space = 5;
 
 
void generateInitialState() {
    for (int i = 0; i < N; i++) {
      for (int j = 0; j < N; j++) {
        U[i][j] = 1.0;
        V[i][j] = 0.0;
        cH[i][j] = (float)(1.3-U[i][j]*1.3);
        cB[i][j] = (float)(1-U[i][j]);
      }
    }
}
 
void setup() {
   size(750,750);
   //fullScreen();
   frameRate(25);
   smooth();
   colorMode(HSB,1.0);
    
   //Set default parameters;
   diffU = 0.16;
   diffV = 0.08;
   paramF = 0.035;
   paramK = 0.06;
    
   rndInitCondition = true;
    
   //Populate U and V with initial data
   generateInitialState();
    
    //Set up offsets
    for (int i = 1; i < N-1; i++) {
      offset[i][0] = i-1;
      offset[i][1] = i+1;
    }
     
    offset[0][0] = N-1;
    offset[0][1] = 1;
     
    offset[N-1][0] = N-2;
    offset[N-1][1] = 0;
}
 
void timestep(float F, float K, float diffU, float diffV) {
      for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
          int p = i + j*N;
           
          float u = U[i][j];
          float v = V[i][j];
           
          int left = offset[i][0];
          int right = offset[i][1];
          int up = offset[j][0];
          int down = offset[j][1];
           
          float uvv = u*v*v;    
        
          float lapU = (U[left][j] + U[right][j] + U[i][up] + U[i][down] - 4*u);
          float lapV = (V[left][j] + V[right][j] + V[i][up] + V[i][down] - 4*v);
           
          dU[i][j] = diffU*lapU  - uvv + F*(1 - u);
          dV[i][j] = diffV*lapV + uvv - (K+F)*v;
        }
      }
       
       
    for (int i= 0; i < N; i++) {
      for (int j = 0; j < N; j++){
          U[i][j] += dU[i][j];
          V[i][j] += dV[i][j];
      }
    }
}
 
void draw(){
    background(color(0.0,0.0,0.1));
    for (int k = 0; k < 20; k++) {
       timestep(paramF, paramK, diffU, diffV);
    }
    
    // Draw points
    noStroke();
    for (int i = 0; i < N; i++) {
      for (int j = 0; j < N; j++) {
        cH[i][j] = (float)(1.3-U[i][j]*1.3);
        cB[i][j] = (float)(1-U[i][j]);
        
        color c = color((float)(1.3-U[i][j]*1.3), 0.87, (float)(1-U[i][j]));
        set(i*space + (width-N*5)/2,j*space +(height-N*5)/2,c);
      }
    }
    for (int i = 0; i < N*space-space; i++) {
      for (int j = 0; j < N*space-space; j++) {
        if((i%space != 0) || (j%space !=0)){
          float x = i;
          float x1 = int(i/space)*space;
          float x2 = x1 + space;

          float y = j;
          float y1 = int(j/space)*space;
          float y2 = y1 + space;
     
          float cH1 = (x2 - x)/(x2 - x1)*cH[int(i/space)][int(j/space)] + ((x - x1)/(x2 - x1))*cH[int(i/space)+1][int(j/space)] ;
          float cH2 = ((x2 - x)/(x2 - x1))*cH[int(i/space)][int(j/space)+1] + ((x - x1)/(x2 - x1))*cH[int(i/space)+1][int(j/space)+1];
          
          float cHH = ((y2 - y)/(y2 - y1))*cH1 + ((y - y1)/(y2 - y1))*cH2;
          
          float cB1 = (x2 - x)/(x2 - x1)*cB[int(i/space)][int(j/space)] + ((x - x1)/(x2 - x1))*cB[int(i/space)+1][int(j/space)] ;
          float cB2 = ((x2 - x)/(x2 - x1))*cB[int(i/space)][int(j/space)+1] + ((x - x1)/(x2 - x1))*cB[int(i/space)+1][int(j/space)+1];
          
          float cBB = ((y2 - y)/(y2 - y1))*cB1 + ((y - y1)/(y2 - y1))*cB2;
          color c = color(cHH,0.87,cBB);
          set(i+ (width-N*5)/2,j+(height-N*5)/2,c);
        }
      }
    }
}
 
 
void keyPressed() {
  switch (key) {
    case '1':
          diffU = 0.16;
          diffV = 0.08;
          paramF = 0.035;
          paramK = 0.06;
          generateInitialState();
          break;
    case '2':
          diffU = 0.16;
          diffV = 0.08;
          paramF = 0.042;
          paramK = 0.065;
          generateInitialState();
          break;
    case '3':
          diffU = 0.18;
          diffV = 0.13;
          paramF = 0.025;
          paramK = 0.056;
          generateInitialState();
          break;
    case '4':
          diffU = 0.18;
          diffV = 0.09;
          paramF = 0.02;
          paramK = 0.056;
          generateInitialState();
          break;
    case '5':
          diffU = 0.14;
          diffV = 0.06;
          paramF = 0.035;
          paramK = 0.065;
          generateInitialState();
          break;
    case '6':
          diffU = 0.19;
          diffV = 0.09;
          paramF = 0.062;
          paramK = 0.062;
          generateInitialState();
          break;
     case '7':
          diffU = 0.16;
          diffV = 0.08;
          paramF = 0.05;
          paramK = 0.065;
          generateInitialState();
          break;
    case 'r':
          rndInitCondition = true;
          generateInitialState();
          break;
    case 'n':
          rndInitCondition = false;
          generateInitialState();
  }
}

void mouseDragged() {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
          if(dist(i,j,(mouseX-(width-N*5)/2)/space,(mouseY-(height-N*5)/2)/space) < 4){
             U[i][j] = 0.5*(1 + random(-1, 1));
             V[i][j] = 0.25*(1 + random(-1, 1));
          }
        }
    }
}

void mousePressed() {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
          if(dist(i,j,(mouseX-(width-N*5)/2)/space,(mouseY-(height-N*5)/2)/space) < 4){
             U[i][j] = 0.5*(1 + random(-1, 1));
             V[i][j] = 0.25*(1 + random(-1, 1));
          }
        }
    }
}