/*
*  reaction-diffusion simulation implementing the Gray-Scott model
 *  main references:
 *  (1) Computational Studies Of Pattern Formation In Turing Systems, Teemu Leppänen, Helsinki University of Technology
 *      http://lib.tkk.fi/Diss/2004/isbn9512273969/
 *  (2) Fronts And Patterns In Reaction-Diffusion Equations, Aric Arild Hagberg, University Of Arizona
 *      http://cnls.lanl.gov/~aric/Simulations/Simulations.html
 *
 *  Frederik Vanhoutte, W:Blut 2009
 */


//connectivity information, stores neighboring cells per index
int[][] connectivity;

//current, previous and next concentration values of substance U
float[] Uc;
float[] Up;
float[] Un;

//current, previous and next concentration values of substance V
float[] Vc;
float[] Vp;
float[] Vn;

//current and previous reaction-diffusion changes of substance U
float[] Fc;
float[] Fp;

//current and previous reaction-diffusion changes of substance V
float[] Gc;
float[] Gp;

// placeholder for array swapping
float[] tmp;

//concentration limits for rendering
float Umin, Umax;
float Vmin, Vmax;

//reaction-diffusion parameters
float F,K; // reaction parameters
float DU,DV; // diffusion coefficients
float dt; // time step

//substrate
int w; // size
PImage substrate; // lo-res image

int steps;
boolean drawing;
int mx,my;


void setup(){
  size(500,500);
  background(0);
  //simulation parameters
  w=120;


  //diffusion coefficients that work well for substrate size and time step
  DU=.082f;
  DV=.041f;


  //reaction parameters that give interesting result
  //found by empirically plotting interesting vs. non-interesting cases in a K vs F phase diagram
  F=random(0.015f,0.06f);
  K=0.05083f+0.26566f*F+random(-0.003f,0.003f);



  tmp = new float[w*w];
  Uc = new float[w*w];
  Up = new float[w*w];
  Un = new float[w*w];
  Fc = new float[w*w];
  Fp = new float[w*w];
  Vc = new float[w*w];
  Vp = new float[w*w];
  Vn = new float[w*w];
  Gc = new float[w*w];
  Gp = new float[w*w];
  connectivity=new int[w*w][];
  substrate=createImage(w,w,ARGB);

  // initial concentrations at trivial steady-state (all V, no U)
  for(int id=0;id<w*w;id++){  
    Uc[id]=0;
    Vc[id]=1;
    Up[id]=Uc[id];      
    Vp[id]=Uc[id];     
  }

  //build connectivity information, first number=number of neighbors, follwoing numbers are the indices of the neighbors.
  for(int i=0;i<w;i++){
    int ip,ipp,im,imm; 

      ip=i+1;
      ipp=i+2;
      im=i-1;
      imm=i-2;



    if(ip>w-1)ip-=w;
    if(ipp>w-1)ipp-=w; 
    if(im<0)im+=w;
    if(imm<0)imm+=w; 
    for(int j=0;j<w;j++){ 


      int id=toIndex(i,j);
      connectivity[id]=new int[9];
      connectivity[id][0]=8;

      int jp,jpp,jm,jmm;

      jp=j+1;
      jpp=j+2;
      jm=j-1;
      jmm=j-2;
      if(jp>w-1)jp-=w;
      if(jpp>w-1)jpp-=w; 
      if(jm<0)jm+=w;
      if(jmm<0)jmm+=w; 

      connectivity[id][1]=toIndex(imm,j);
      connectivity[id][2]=toIndex(im,j);
      connectivity[id][3]=toIndex(ip,j);
      connectivity[id][4]=toIndex(ipp,j);
      connectivity[id][5]=toIndex(i,jmm);
      connectivity[id][6]=toIndex(i,jm);
      connectivity[id][7]=toIndex(i,jp);
      connectivity[id][8]=toIndex(i,jpp);



    }
  }

  // take first step using only current state (Euler)
  Euler(1f);

  steps=10;
  drawing=false;
  colorMode(HSB,1.0f);
}

void draw(){

  if(drawing){
    // add some V at mouse position
    mx=(int)(mouseX*w/(float)width);
    my=(int)(mouseY*w/(float)height);
    if((mx>=0)&&(mx<w)&&(my>=0)&&(my<w)){
      Vc[mx+w*my]=min(Vc[mx+w*my]+.5f,.9f);

    }
  }
  else{
    mx=my=-1; 
  }

  // take 10 steps using previous and current state (Adams-Bashforth)

  //for(int i=0;i<steps;i++){
    //AdamsBashforth(0.75f);
  //}
  
  Euler(1f);

  // find current range of concentrations and adapt color range (but ignore mouse position)
  checkRange(mx,my);
  if(Umin<Umax){
    substrate.loadPixels();
    for(int id=0;id<w*w;id++){
      substrate.pixels[id]=color(1.1f-(Uc[id]-Umin)/(Umax-Umin));
    }
    substrate.updatePixels();
    image(substrate,0,0,width,height);
  }
}


// 1st order Euler time step
void Euler(float dt){ 
  int id=0;
  for(int i=0;i<w;i++){
    for(int j=0;j<w;j++){ 
      float[] Fres=F(Uc,Vc,j,i); // calculate reaction-diffusion
      Fc[id]=Fres[0];
      Gc[id]=Fres[1];
      Un[id]=Uc[id]+dt*Fc[id]; 
      Vn[id]=Vc[id]+dt*Gc[id];  
      id++;
    }
  }
  // swap current and next arrays
  tmp=Uc;
  Uc=Un;
  Un=tmp;
  tmp=Vc;
  Vc=Vn;
  Vn=tmp;
  tmp=Fp;
  Fp=Fc;
  Fc=tmp;
  tmp=Gp;
  Gp=Gc;
  Gc=tmp;
}

// 2nd order Adams-Bashforth time step, less sensitive to numerical errors thane Euler
void AdamsBashforth(float dt){
  int id=0;
  for(int i=0;i<w;i++){
    for(int j=0;j<w;j++){ 

      float[] Fres=F(Uc,Vc,j,i); // calculate reaction-diffusion
      Fc[id]=Fres[0];
      Gc[id]=Fres[1];
      Un[id]=Uc[id]+1.5f*dt*Fc[id]-0.5f*dt*Fp[id]; 
      Vn[id]=Vc[id]+1.5f*dt*Gc[id]-0.5f*dt*Gp[id]; 
      id++;
    }
  }
  // swap previous, current and next arrays
  tmp=Up;
  Up=Uc;
  Uc=Un;
  Un=tmp;
  tmp=Vp;
  Vp=Vc;
  Vc=Vn;
  Vn=tmp;
  tmp=Fp;
  Fp=Fc;
  Fc=tmp;
  tmp=Gp;
  Gp=Gc;
  Gc=tmp;
}


//reaction-diffusion: change of concentration = change by diffusion + change by reaction
// model: Gray-Scott
float[] F(final float U[],final float V[], int i, int j){
  int id=i+w*j;
  float[] result = new float[2];  
  float product=U[id]*V[id]*V[id];
  result[0]=DU*L4(U,id)-product+F*(1f-U[id]);
  result[1]=DV*L4(V,id)+product-(F+K)*V[id];
  return result;
}


// Diffusion is driven by concentration differences. In the equations these are expressed by the Laplacian operator

//4th order approximation of Laplacian 
float L4(final float U[], int id){
  return 1/12f*(-U[connectivity[id][1]]+16*U[connectivity[id][2]]-30*U[id]+16*U[connectivity[id][3]]-U[connectivity[id][4]]-U[connectivity[id][5]]+16*U[connectivity[id][6]]-30*U[id]+16*U[connectivity[id][7]]-U[connectivity[id][8]]);
}

//2nd order approximation of Laplacian, slightly faster but more prone to numerical errors
float L2(final float U[], int id){
  return U[connectivity[id][2]]+U[connectivity[id][3]]+U[connectivity[id][6]]+U[connectivity[id][7]]-4*U[id]; 
}



// find minimum and maximum concentration of substances U and V
void checkRange (int i, int j)
{
  Umin = 100;
  Umax = -1;
  Vmin = 100;
  Vmax = -1;
  for ( int id=0; id < w*w; id++ ){
    if(id!=i+w*j){
      Umin = min (Umin, Uc[id]);
      Umax = max (Umax, Uc[id]);
      Vmin = min (Vmin, Vc[id]);
      Vmax = max (Vmax, Vc[id]);
    }
  }
} 



void keyPressed(){
  //reinitiliaze everything (identical to setup())
  background(0);

  F=random(0.015f,0.06f);
  K=0.05083f+0.26566f*F+random(-0.003f,0.003f);
  tmp = new float[w*w];
  Uc = new float[w*w];
  Up = new float[w*w];
  Un = new float[w*w];
  Fc = new float[w*w];
  Fp = new float[w*w];
  Vc = new float[w*w];
  Vp = new float[w*w];
  Vn = new float[w*w];
  Gc = new float[w*w];
  Gp = new float[w*w];

  substrate=createImage(w,w,ARGB);
  for(int id=0;id<w*w;id++){
    Uc[id]=0;
    Vc[id]=1;
    Up[id]=Uc[id];      
    Vp[id]=Uc[id];      
  }

  Euler(1f);

}




// slow down update while drawing
void mousePressed(){
  steps=1;
  drawing=true; 

}

void mouseReleased(){
  drawing=false;
  steps=10; 
}

int toIndex(int i, int j){
  return i+w*j; 
}