/*
original code by Staross ( http://openprocessing.org/portal/?userID=6146 )

dt = 0.004;
Nt = 1000;
E = 350;
Nx = 1000;
dx = E/Nx;

A = zeros(Nx,Nx); // B = zeros(m,n) or B = zeros([m n]) returns an m-by-n matrix of zeros.

%letter as initial cond
A = initLetter;

i = sqrt(-1);

alpha = 4;
D = 0.4;
beta = 6;

%diffusion
H = [0 1  0;
     1 -4 1;
     0 1  0;];   
  
 H = D*H / dx^2;
 
for t=1:Nt
        
    D2A = imfilter(A,H); // B = imfilter(A, H) filters the multidimensional array A 
         with the multidimensional filter H. The array A can be logical or a nonsparse 
         numeric array of any class and dimension. The result B has the same size and class as A.
         Each element of the output B is computed using double-precision floating point. 
         If A is an integer or logical array, then output elements that exceed the range of 
         the integer type are truncated, and fractional values are rounded.
    
    A = A + dt*(A +(1+i*alpha).*D2A - (1-i*beta)*( abs(A).^2 ).*A );
    
end

imagesc(real(A)); // The imagesc function scales image data to the full range of the current 
     colormap and displays the image

*/

/*

The Complex Ginzburg-Landau equation
adapted from Staross matlab code

see http://codeinthehole.com/tutorials/cgl/index.html

visual output examples from this code : http://codelab.fr/2027

emoc, dec. 2010
*/

float dt = 0.004;  // à l'origine 0.004
float E = 80; // <-- changer en fonction de la taille de l'image
int Nx = 300; // taille de l'image (carrée)
float dx = E/float(Nx); // à l'origine 0.35 

float D = 0.4;    // à l'origine 0.4
float aalpha = 4; // à l'origine 4
float beta = 6;   // à l'origine 6

PImage base_image;
PFont police;

float[][] H = { { 0.0,  1.0, 0.0 },
                { 1.0, -4.0, 1.0 },
                { 0.0,  1.0, 0.0 } }; 
 
int Hsize = 3;

int t; // temps
int mode = 6; // mode d'affichage
String SKETCH_NAME = "ginzland005";

// matrice A de nombres complexes
float[][] A = new float[Nx][Nx]; // partie réelle
float[][] I = new float[Nx][Nx]; // partie imaginaire

// matrice D2A de nombres complexes
float[][] D2A = new float[Nx][Nx]; // partie réelle
float[][] D2I = new float[Nx][Nx]; // partie imaginaire


void setup() {
  size(400, 400, P2D);
  base_image = loadImage("S.png");
  police = loadFont("OhLaLa-10.vlw");
  
  // initialiser les tableaux
  initialiserTableau(A, 0);
  initialiserTableau(I, 0);
  initialiserTableau(D2A, 0);
  initialiserTableau(D2I, 0);
  
  // modifier la matrice H
  H = matrixMod(H, D, dx, Hsize);

  // Charger les pixels de l'image dans A et I // noir et blanc // normalisé de 0 à 1
  for (int i = 0; i < Nx; i++) {
    for (int j = 0; j < Nx; j++) {
      color cloc = base_image.get(i, j);
      float val = ((red(cloc) + green(cloc) + blue(cloc)) / 3 ) / 255;
      A[i][j] = val;
      I[i][j] = val;
    }
  }
}

void draw() {
  colorMode(RGB, 1);
  background(1, 1, 1);
  textFont(police, 10);
  fill(1, 0, 0);
  //text("mode : " + mode + " - génération : " + t, 20, 380);
  text("mode : " + mode + " - génération : " + t + " - durée " + (millis() / 1000) + " s", 20, 365);
  text("D : " + D + " aalpha : " + aalpha + " beta : " + beta + " dx : " + dx, 20, 380); 
  text("dt : " + dt, 20, 395);
  
  
  // convolution de A par H (en évitant les bords...) partie réelle et imaginaire
  for (int i = 1; i < Nx-1; i++) {
    for (int j = 1; j < Nx-1; j++) {
      D2A[i][j]  = (A[i-1][j-1] * H[0][0]) + (A[i][j-1] * H[1][0]) + (A[i+1][j-1] * H[2][0]); 
      D2A[i][j] += (A[i-1][j]   * H[0][1]) + (A[i][j]   * H[1][1]) + (A[i+1][j]   * H[2][1]);
      D2A[i][j] += (A[i-1][j+1] * H[0][2]) + (A[i][j+1] * H[1][2]) + (A[i+1][j+1] * H[2][2]);
    }
  }
  for (int i = 1; i < Nx-1; i++) {
    for (int j = 1; j < Nx-1; j++) {
      D2I[i][j]  = (I[i-1][j-1] * H[0][0]) + (I[i][j-1] * H[1][0]) + (I[i+1][j-1] * H[2][0]); 
      D2I[i][j] += (I[i-1][j]   * H[0][1]) + (I[i][j]   * H[1][1]) + (I[i+1][j]   * H[2][1]);
      D2I[i][j] += (I[i-1][j+1] * H[0][2]) + (I[i][j+1] * H[1][2]) + (I[i+1][j+1] * H[2][2]);
    }
  }
  
  // modifier A
  for (int i = 0; i < Nx; i++) {
    for (int j = 0; j < Nx; j++) {
      A[i][j] = A[i][j] + dt * ( A[i][j] + D2A[i][j] - aalpha * D2I[i][j] - (pow(A[i][j], 2) + pow(I[i][j], 2)) * (A[i][j] + beta * I[i][j]));
      I[i][j] = I[i][j] + dt * ( I[i][j] + D2I[i][j] + aalpha * D2A[i][j] - (pow(A[i][j], 2) + pow(I[i][j], 2)) * (I[i][j] - beta * A[i][j])); 
    }
  }
  
  
  // afficher l'image à partir de la partie réelle de A
  translate(50, 50);
  for (int i = 0; i < Nx; i++) {
    for (int j = 0; j < Nx; j++) {
      /* DEBUG : afficher quelques valeurs brutes de A
      if (t == 100) {
        if ((i > 150) && (i < 160)) println(A[i][j]);
      }*/
      float c = modeAffichage(A[i][j], mode);
      stroke(c, c, c);
      point(i, j);
    }
  }
  println("t : " + t + " millis " + millis());
  t++;
}

float[][] matrixMod(float[][] m, float D, float dx, int matrixsize) {
  float[][] mnew = new float[matrixsize][matrixsize];
  for (int i = 0; i < matrixsize; i++){
    for (int j= 0; j < matrixsize; j++){
      mnew[i][j] = (m[i][j] * D) / (pow(dx , 2));
      println(mnew[i][j]);
    }
  }
  return mnew;
}

float[][] initialiserTableau(float[][] t, float valeur) {
  for (int i = 0; i < t.length; i++){
    for (int j= 0; j < t[0].length; j++){
      t[i][j] = valeur;
    }
  }
  return t;
}

float modeAffichage(float v, int mode) {
  switch(mode) {
    case 0: v = v;            break;
    case 1: v = v * 0.2 + 1;  break;
    case 2: v = v * 10;       break;
    case 3: v = v * 100;      break;  
    case 4: v = v * 1000;     break;  
    
    case 5: 
      if (v < 0) v *= -1;
      v = v % 1;
      break;
    
    case 6: 
      if (v < 0) v *= -1;
      v = (v % 4) * 0.25;
      break;

    case 7: v = (v + 5) * 0.2; break;
    case 8: v = v * v * v;     break;
    case 9: v = map(v, -2, 2, 1, 0); break;
  }
  return v;
}

void keyPressed() {
  if (key == 's') {
    saveFrame(SKETCH_NAME+"_"+minute()+second()+millis()+"mode"+mode+".png");
  }
  if (key == '0') { mode = 0; println("mode : " + mode); }
  if (key == '1') { mode = 1; println("mode : " + mode); }
  if (key == '2') { mode = 2; println("mode : " + mode); }
  if (key == '3') { mode = 3; println("mode : " + mode); }
  if (key == '4') { mode = 4; println("mode : " + mode); }
  if (key == '5') { mode = 5; println("mode : " + mode); }
  if (key == '6') { mode = 6; println("mode : " + mode); }
  if (key == '7') { mode = 7; println("mode : " + mode); }
  if (key == '8') { mode = 8; println("mode : " + mode); }
  if (key == '9') { mode = 9; println("mode : " + mode); }
}

/* nombres complexes et java

http://en.literateprograms.org/Complex_numbers_%28Java%29
http://commons.apache.org/math/apidocs/org/apache/commons/math/complex/Complex.html

syntaxe matlab :

http://www.cs.princeton.edu/introcs/11matlab/
getting started with matlab : http://www.mathworks.com/help/pdf_doc/matlab/getstart.pdf
talab tutorial(s) : http://www.mathworks.com/academia/student_center/tutorials/launchpad.html

open-processing.org 
diffusion-reaction (Staross) : 10133
Belousov-Zhabotinsky (BZ) reaction : 9813
reaction-diffusion : 6195
mandelbrot : 6343
*/