gistfile1.txt

// Occasionallly massive slow down due to ...? http://en.wikipedia.org/wiki/Denormal

Grid gri;
int tm;

public void setup() {
  size(800, 800);
  colorMode(HSB, 1.0, 1.0, 1.0);
  // Da, Db, diffusion rates for a, b
  // f = feed, k = kill
  //               width, height    a    b   Da   Db       f       k
  //gri = new Grid(width, height, 1.0, 0.0, 1.0, 0.5, 0.0315, 0.0606);
  //gri = new Grid(width, height, 1.0, 0.0, 1.0, 0.5, 0.01747, 0.0432);
  //gri = new Grid(width, height, 1.0, 0.0, 1.0, 0.3, 0.01747, 0.0432);
  //gri = new Grid(width, height, 1.0, 0.0, 1.0, 0.5, 0.055, 0.062);
  //gri = new Grid(width, height, 1.0, 0.0, 1.0, 0.3, 0.054, 0.062);
  //gri = new Grid(width, height, 1.0, 0.0, 1.0, 0.5, 0.054, 0.062);
  gri = new Grid(width, height, 1.0, 0.0, 1.0, 0.5);//, 0.045,0.0635);
  //gri = new Grid(width, height, 1.0, 0.0, 1.0, 0.5,0.0173, 0.0498);
  
//0.017281079664826393:0.04979756847023964
//0.017390167340636253:0.0498465895652771
  //gri.rand();
  
  //gri.circle(width>>1, height>>1, 100, 0.1, 1.0);
  //gri.circle(width>>1, (height>>1)-25, 70, 1.0, 0.0);
  //gri.rectangle((width>>1)-5, (height>>1)-5, 10,10,0.1,1.0);
}

public void draw() {
  
  gri.edge(0.0, 0.0);                        // clear edges so they act like they aren't there
  for (int i = 0; i < 5; i++) gri.update();

  loadPixels();                              // draw grid  // lazy, making the assumption that the grid is the same size as the processing window
  for (int i = 0; i < gri.sw; i++) {
    double a = gri.a[gri.current][i];
    double b = gri.b[gri.current][i];
    //pixels[i] = color((float)(a+b)*0.5, 1.0, 1.0);
    pixels[i] = color(1.0-(float)a, 1.0-(float)a, 1.0-(float)b);
    //pixels[i]= color(0, 0, b);
  }
  updatePixels();
}

public class Grid {
  boolean kfMap;         
  double corn = 0.05;     // weights for laplace  
  double side = 0.2;   
  double mid = -1;    
  double[][] a, b;   
  int wd, ht, sw;   
  int current = 0;   
  double ts = 1;//.5;        // lower if the diffusion reaction explodes
  double la, lb, Da, Db, f, k;      
  double[] kmap, fmap;
  public Grid(int wd_, int ht_, double aa, double bb, double Da_, double Db_, double f_, double k_) {    // create a grid with single k and f values
    wd = wd_;
    ht = ht_;
    sw = wd*ht;
    a = new double[2][sw];
    b = new double[2][sw];
    Da = Da_;
    Db = Db_;
    f = f_;
    k = k_;  
    init(aa, bb);
    kfMap = false;
  }

  public Grid(int wd_, int ht_, double aa, double bb, double Da_, double Db_) {  // create a grid, with multiple k, f values // set this in initkf()
    this(wd_, ht_, aa, bb, Da_, Db_, 0, 0);
    kfMap = true;
    fmap = new double[sw];
    kmap = new double[sw];
    initkf();
  }
  public void update() {           // update, all drawing should be done before this if it is to show up in this update
    for (int y = 1; y < ht-1; y++) {  // loop from 1 to  1 less than width/height, (makes the laplace easier)
      int ym = y*wd;               
      for (int x = 1; x < wd-1; x++) {
        int index = x+ym;
        la = lap(a[current], index);      // laplace 
        lb = lap(b[current], index);
        double a_ = a[current][index];     // for read-ability below 
        double b_ = b[current][index];
        double abb = a_*b_*b_;
        if (kfMap) {            // if using custom values per index
          f = fmap[index];
          k = kmap[index];
        }
        // www.karlsims.com/rd.html  and   mrob.com/pub/comp/xmorphia/    U = a, V = b
        // Gray-Scott Model
        // Da*la, Db*lb = diffusion term + = increase - = decrease
        // -abb = reaction rate
        // a converts to b, a gets replenished
        //f(1-a) = replenishment term, f = rate
        //(k+f)*b = dimishment term, k + f represents permeability of the membrane to 'a'
        //t(0.5) = time
        a[1-current][index] = a_ + (Da*la - abb + f*(1-a_))*ts;
        b[1-current][index] = b_ + (Db*lb + abb - (k + f)*b_)*ts;
      }
    }
    current = 1-current;
  }

  public void initkf() {                // set a grid of k and f values, used if missing the last two grid constructor variables
    for (int x = 0; x < wd; x++) {
      for (int y = 0; y < ht; y++) {
        int index = x + y*wd;
        //fmap[index] = map(x, 0, wd, 0, 0.26);
        //kmap[index] = map(y, 0, ht, 0, 0.100);
        //fmap[index] = map(x, 0, wd, 0.016385, 0.018333);
        //kmap[index] = map(y,0.0,ht,0.049368,0.050013);
        //fmap[index] = sin(x*0.02)*0.01 + 0.05;
        //kmap[index] = sin(y*0.02)*0.005 + 0.06;
        float d = map(dist(x,y,wd*0.5, ht*0.5),0,sqrt(wd*wd+ht*ht)*0.5,0.0,1.0);
        d = (1.0/(d+1.0)-0.5)*2;//(d < 0.5) ? 0.0 : 1.0;
        fmap[index] = map(d, 0, 1, 0.047314748, 0.05962701);
        kmap[index] = map(d, 0, 1, 0.061306357, 0.06479236);
      }
    }
  }

  public double lap(double[] c, int index) {      // part of the diffusion term
    //       top left           top right          bottom left          bottom right                 top             left             right        bottom                    middle      
    return (c[index - 1 -wd] + c[index + 1 -wd] + c[index - 1 + wd] + c[index + 1 + wd]) * corn + (c[index  -wd] + c[index - 1 ] + c[index + 1 ] + c[index  + wd]) * side + c[index]*mid;
    
  }

  public void rand() {                        // seed grid with random a and b values 
    for (int i = 0; i < sw; i++) {
      
      b[current][i] = random(0.25)+0.75;
      a[current][i] = 1.0 - b[current][i];
    }
  }

  public void init(double aa, double bb) {      // seed the grid with values of a and b
    for (int i = 0; i < sw; i++) {
      a[current][i] = aa;
      b[current][i] = bb;
    }
  }
  public void rectangle(int x_, int y_, int w_, int h_, double aa, double bb) {   // set all elements in a rectangle to aa, bb
    for (int y = y_; y < y_+h_; y++) {
      if (y < 0 || y >= ht) continue;
      for (int x = x_; x < x_+w_; x++) {
        if (x < 0 || x >= wd) continue;
        a[current][x+y*wd] = aa;
        b[current][x+y*wd] = bb;
      }
    }
  }
  public void circle(int x_, int y_, int r_, double aa, double bb) {   // set all elements in a circle to aa, bb
    int rr = r_*r_;
    for (int y = -r_; y <= +r_; y++) {
      int ym = y+y_;
      if (ym < 0 || ym >= ht) continue;
      int xr = (int)sqrt(rr-y*y);
      for (int x = -xr; x<= +xr; x++) {
        int xm = x+x_;
        if (xm < 0 || xm >= wd) continue;
        a[current][xm+ym*wd] = aa;
        b[current][xm+ym*wd] = bb;
      }
    }
  }
  
  public void edge(double aa, double bb) {  // set the edges of the grid to aa, bb
    for (int x = 0; x < wd; x++) {
      a[current][x] = aa;
      a[current][x+(ht-1)*wd] = aa;
      b[current][x] = bb;
      b[current][x+(ht-1)*wd] = bb;
    }
    for (int y = 0; y < wd; y++) {
      a[current][y*wd] = aa;
      a[current][wd-1 + y*wd] = aa;
      b[current][y*wd] = bb;
      b[current][wd-1 + y*wd] = bb;
    }
  }

  public void cellset(int x, int y, double aa, double bb) {             // set a single element to aa, bb
    int index = constrain(x, 0, wd-1) + constrain(y, 0, ht-1)*wd;
    a[0][index] = aa;
    a[1][index] = aa;
    b[0][index] = bb;
    b[1][index] = bb;
  }
}

public void mousePressed() {
  //println((float)mouseX/gri.wd+":"+(float)mouseY/gri.ht);
  int index = mouseX + mouseY*gri.wd;
  //if (gri.kfMap) println(gri.fmap[index]+":"+gri.kmap[index]);
}

public void mouseDragged() {
  if (mouseButton == LEFT) {
    gri.cellset(mouseX, mouseY, 0.1, 1.0);
  } else {
    gri.cellset(mouseX, mouseY, 1.0, 0.1);
  }
}