class ControlPoint {
  float x, y, z, ox, oy, oz;
  ControlPoint() {
  }
  void create() {
    ox = x = random(-1.0,1.0);
    oy = y = random(-1.0,1.0);
    oz = z = random(-1.0,1.0);
  }
  void create(Extent ex) {
    ox = x = random(ex.minx, ex.maxx);
    oy = y = random(ex.miny, ex.maxy);
    oz = z = random(ex.minz, ex.maxz);
  }
  void draw(Extent ex) {
    float scrx = x - ex.ctrx;
    float scry = y - ex.ctry;
    float scrz = z - ex.ctrz;
    float h = 128.0 + (atan2(y-oy, x-ox)/PI) * 128.0;
    float e = sphereSize * ex.width/width;
    pushMatrix();
    translate(scrx,scry,scrz);
    switch(renderMode) {
      case rmPOINTS:
        stroke(color(h,160,160));
        point(0,0);
        break;
      case rmRECTS:
        noStroke();
        fill(color(h,160,160));
        rect(-e,-e, e,e);
        break;
      case rmBOXES:
        noStroke();
        fill(color(h,160,160));
        box(e*2.0);
        break;
      case rmSPHERES:
        noStroke();
        fill(color(h,160,160));
        sphere(e);
        break;
    }
    popMatrix();
  }
} 

class Extent {
  float minx,miny,maxx,maxy,minz,maxz;
  float ctrx,ctry,ctrz;
  float width,height,depth;
  Extent() {
  }
  void calc() {
    width = maxx - minx;
    height = maxy - miny;
    depth = maxz - minz;
    ctrx = (minx+maxx)*0.5;
    ctry = (miny+maxy)*0.5;
    ctrz = (minz+maxz)*0.5;
  }
  void untrack() {
    minx = miny = minz = 99;
    maxx = maxy = maxz = -99;
  }
  void track(float x, float y, float z, float trackrate) {
    if (x < minx) minx += (x-minx) * trackrate;
    if (x > maxx) maxx += (x-maxx) * trackrate;
    if (y < miny) miny += (y-miny) * trackrate;
    if (y > maxy) maxy += (y-maxy) * trackrate;
    if (z < minz) minz += (z-minz) * trackrate;
    if (z > maxz) maxz += (z-maxz) * trackrate;
    calc();
  }
  boolean valid() {
    calc();
    if ((minx <= -1E9) || (miny <= -1E9) || (minz <= -1E9)) return false;
    if ((maxx >= 1E9) || (maxy >= 1E9) || (maxz >= 1E9)) return false;
    return true;
  }
}

class IfsFunction {
  static final int ncoefs = 12;
  float [] coefs;
  float prob;
  IfsFunction() {
    coefs = new float[ncoefs]; 
  }
  void create() {
    for (int i=0; i<ncoefs; i++)
      coefs[i] = random(-PI,PI);
    prob = random(1.0); 
  }
  void apply(ControlPoint cp) {
    // retrieve current values
    float x = cp.x;
    float y = cp.y;
    float z = cp.z;
    // calculate new position
    float xp =  cos(coefs[0]*x) + sin(coefs[1]*y) - sin(coefs[2]*z) - coefs[3];
    float yp =  sin(coefs[4]*x) + cos(coefs[5]*y) + sin(coefs[6]*z) - coefs[7];
    float zp =  cos(coefs[8]*x) - cos(coefs[9]*y) + cos(coefs[10]*z) - coefs[11];
    // update historical values
    cp.ox = x;
    cp.oy = y;
    cp.oz = z;
    // update current values
    cp.x = xp;
    cp.y = yp;
    cp.z = zp;
  }
}

class IfsSystem {
  IfsFunction [] funcs;
  int nfuncs;
  ControlPoint cp;
  Extent ex;
  IfsSystem(int n) {
    nfuncs = n;
    funcs = new IfsFunction[nfuncs];
    for (int i=0; i<nfuncs; i++)
      funcs[i] = new IfsFunction();
    cp = new ControlPoint();
    ex = new Extent();
  }
  void create() {
    float sumprob = 0.0;
    for (int i=0; i<nfuncs; i++) {
      funcs[i].create();
      sumprob += funcs[i].prob;
    }
    for (int i=0; i<nfuncs; i++)
      funcs[i].prob /= sumprob; 
    estimate(100,200,100);
  }
  boolean valid() {
    float exmin = 0.001;
    float exmax = 100;
    if ((ex.width < exmin) || (ex.height < exmin) || (ex.depth < exmin)) return false;
    if ((ex.width > exmax) || (ex.height > exmax) || (ex.depth > exmax)) return false;
    return true;
  }
  IfsFunction pickfunc() {
    float sumprob = 0.0;
    int which = 0;
    float r = random(1.0);
    for (int i=0; i<nfuncs; i++) {
      sumprob += funcs[i].prob;
      if (r < sumprob)
        which = i;
    }
    return funcs[which];    
  }
  void estimate(int ncp, int nprime, int niters) {
    ex.untrack();
    for (int i=0; i<ncp; i++) {
      cp.create();
      for (int j=0; j<nprime; j++) {
        IfsFunction f = pickfunc();
        f.apply(cp);
      }
      for (int j=0; j<niters; j++) {
        IfsFunction f = pickfunc();
        f.apply(cp);
        ex.track(cp.x,cp.y,cp.z,1.0);
      }
    }
  }    
  void iterate(int ncp, int nprime, int niters) {
    for (int i=0; i<ncp; i++) {
      cp.create();
      for (int j=0; j<nprime; j++) {
        IfsFunction f = pickfunc();
        f.apply(cp);
      }
      for (int j=0; j<niters; j++) {
        IfsFunction f = pickfunc();
        f.apply(cp);
        cp.draw(ex);
      }
    }
  }   
}