//bz reaction by Alasdair Turner 
//http://www.openprocessing.org/visuals/?visualID=1263
class BZ{
float [][][] a;
float [][][] b;
float [][][] c;
int sz;
int p = 0, q = 1;
int time=0;
float ch=0.05,ad=2; 
BZ()
{
sz=30; 
  a = new float [sz][sz][2];
  b = new float [sz][sz][2];
  c = new float [sz][sz][2];
  for (int x = 0; x < sz; x++) {
    for (int y = 0; y < sz; y++) {
      a[x][y][p] = random(0.0,0.7);
      b[x][y][p] = random(0.0,0.7);
      c[x][y][p] = random(0.0,0.7);
    }
  }
} 
void run()
{
  time++;
  for (int x = 0; x < sz; x++) {
    for (int y = 0; y < sz; y++) {
      float c_a = 0.0;
      float c_b = 0.0;
      float c_c = 0.0;
      for (int i = x - 1; i <= x+1; i++) {
        for (int j = y - 1; j <= y+1; j++) {
          c_a += a[(i+sz)%sz][(j+sz)%sz][p];
          c_b += b[(i+sz)%sz][(j+sz)%sz][p];
          c_c += c[(i+sz)%sz][(j+sz)%sz][p];
        }
      }
      c_a /= 9.0*(noise((frameCount)*ch,x*ch,y*ch)*ad);
      c_b /= 9.0*(noise((frameCount)*ch,x*ch,y*ch)*ad);
      c_c /= 9.0*(noise((frameCount)*ch,x*ch,y*ch)*ad);
      // adjust these values to alter behaviour
      a[x][y][q] = constrain(c_a + c_a * (c_b - c_c), 0, 1);
      b[x][y][q] = constrain(c_b + c_b * (c_c - c_a), 0, 1);
      c[x][y][q] = constrain(c_c + c_c * (c_a - c_b), 0, 1);
    }
  }
              if (p == 0) {
    p = 1; q = 0;
  }
  else {
    p = 0; q = 1;
  }


}
}

//compilation of some codes by respected gentlemen
//most of the stuff by toxi
//echoechonoisenoise 2010
//open architecture open design : use it for any purpose, just link the results of your work

import peasy.org.apache.commons.math.*;
import peasy.*;
import peasy.org.apache.commons.math.geometry.*;
import toxi.geom.*;
import toxi.volume.*;
import toxi.math.noise.*;

PVector l[];//lights
float str;

BZ bz;
boolean runReaction=true;

int DIMX=30;//for explanation of isosurface creation, properties, methods etc, check toxiclibs examples
int DIMY=30;
int DIMZ=30;
int sc=400;
float t=0.007;
float ISO_THRESHOLD = 0.3;
float NS=0.3;
Vec3D SCALE=new Vec3D(1,1,1).scaleSelf(sc);
IsoSurface surface;
float[][][] PastVolumeData=new float[DIMX][DIMY][DIMZ];
VolumetricSpace volume=new VolumetricSpace(SCALE,DIMX,DIMY,DIMZ);

PeasyCam cam;
void setup() {
  size(600,600,P3D);
  colorMode(HSB,255);
  boolean w=true;
  bz=new BZ();
  initLight();
  stroke(0,20);
  cam=new PeasyCam(this,800);
}
void draw() {
  background(55);
  for(int i=0; i<l.length; i++) {
    directionalLight(str,str,str, l[i].x,l[i].y,l[i].z);
  }
  int index=0;
  if(runReaction){
    bz.run();//progress bz reaction simulation 
    float[] volumeData=volume.getData();
    for(int z=0;z<DIMZ; z++) {
      for(int y=0; y<DIMY; y++) {
        for(int x=0; x<DIMX; x++) {
          if(z<DIMZ-1){
            volumeData[index++]=PastVolumeData[x][y][z];
            PastVolumeData[x][y][z]=PastVolumeData[x][y][z+1];
          }
          else{
            volumeData[index++]=bz.a[x][y][bz.q];
            PastVolumeData[x][y][z]=bz.a[x][y][bz.q];
          }
        } 
      }
    }
    volume.closeSides();
    surface=new IsoSurface(volume);
    surface.computeSurface(ISO_THRESHOLD);
  }
  beginShape(TRIANGLES);
  int num=surface.getNumFaces(); 
  Vec3D[] v=null;
  for(int i=0; i<num; i++) {
    v=surface.getVerticesForFace(i,v);
    fill(PastVolumeData[(int)map(v[0].x,-sc,sc,0,DIMX-1)][(int)map(v[0].y,-sc,sc,0,DIMY-1)][(int)map(v[0].z,-sc,sc,0,DIMZ-1)]*100,100,150);//didn't figure out a better way yet...
    vertex(v[0].x,v[0].y,v[0].z);
    fill(PastVolumeData[(int)map(v[1].x,-sc,sc,0,DIMX-1)][(int)map(v[1].y,-sc,sc,0,DIMY-1)][(int)map(v[1].z,-sc,sc,0,DIMZ-1)]*100,100,150);
    vertex(v[1].x,v[1].y,v[1].z);
    fill(PastVolumeData[(int)map(v[2].x,-sc,sc,0,DIMX-1)][(int)map(v[2].y,-sc,sc,0,DIMY-1)][(int)map(v[2].z,-sc,sc,0,DIMZ-1)]*100,100,150);
    vertex(v[2].x,v[2].y,v[2].z);
  }
  endShape();
}
void keyTyped(){
  runReaction=!runReaction;  
}
void initLight() {//by Marius Watz
  randomSeed(0);
  l=new PVector[7];
  for(int i=0; i<l.length; i++) {
    str=random(TWO_PI);
    l[i]=new PVector(cos(str)*10,0.3,sin(str)*10);
  }
  str=random(120,180);
}