//W:Mute 2004 - Abstract
//www.wmute.org

//Modified Peter de Jong - folded three-plane attractor
//Based on applet by j.tarbell, complexification.net



//drawing space
int dimx = 800;
int dimy = 800;
int dimz = 800;
BImage bkg;

// 6 parameters for the extended Peter de Jong attractor
float[] xPdJ = new float[6];
float pmin=-2.5;
float pmax= 2.5;

// scaling parameters for display (scale, offset and rotation)
float gsx, gsy, gsz, gx, gy, gz;
float randomscale = 0;
float anglex=0;
float angley=0;
float anglez=0;

// starting point limits
float xneg=-1.0;
float xpos=+1.0;
float yneg=-1.0;
float ypos=+1.0;
float zneg=-1.0;
float zpos=+1.0;

//Parameters for TravelerHenon class
int maxage = 64;
int maxnum = 1000;
int maxcounts=400;
Passenger[] passengers = new Passenger[maxnum]; // based on TravelerHenon class by j.tarbell


//stores explored points
float points[][] = new float[maxcounts*maxnum][3];
color colors[]=new color[maxcounts*maxnum];
int counter=0;

// program flags
boolean resetflag = false; // restart with new parameters
boolean firstflag=true; // toggles clear screen
boolean drag;
ArcBall arcball;
BufferedMouse bufferedMouse;
boolean mouseClick=false;

void setup()
{
  size(640,400);
  bkg=loadImage("bkg.jpg");
  background(bkg);
  fill(0);
  noStroke();
  xPdJ[0] = 2.01;
  xPdJ[1] = -2.53;
  xPdJ[2] = 1.61;
  xPdJ[3] = -0.33 ;
  xPdJ[4] = 0.33;
  xPdJ[5] = 1.55;
  for (int i=0;i<maxnum;i++) {
    passengers[i] = new Passenger(random(xneg,xpos),random(yneg,ypos),random(zneg,zpos),true);
  }
  arcball = new ArcBall(width / 2.0f, height / 2.0f, min(width - 25, height - 25) / 2.0f);
  bufferedMouse=new BufferedMouse(0.5);
  drag=false;
 
}

void loop()
{
  if (counter<(maxcounts*maxnum)) {
    for (int i=0;i<maxnum;i++) {
      passengers[i].update();
    }
    if ((resetflag)&&(!drag)) {
      resetflag=false;
      reset();
    }

  }
  if (counter==(maxcounts*maxnum)){
    display();
    if ((resetflag)&&(!drag)) {
      resetflag=false;
      reset();
    }
  }
}



void keyPressed()
{
   resetflag=true;
}

void mousePressed(){
  arcball.mousePressed();
  mouseClick=true;
}

void mouseDragged()
{
  arcball.mouseDragged();
  
  if (mouseClick) drag=true;
}

void mouseReleased()
{
  mouseClick=false;
  drag=false;
  firstflag=true;

}

void display(){
  translate(width/2,height/2,0);
  translate(gx*dimx,gy*dimy,gz*dimz);
  bufferedMouse.update();
  arcball.run();
  if (drag) {
    background(255);
     for (int i=0;i<maxcounts*maxnum;i+=100){   
        stroke(red(colors[i]),green(colors[i]), blue(colors[i]),200);
        point(points[i][0]/gsx*0.95*dimx,points[i][1]/gsy*0.95*dimy,points[i][2]/gsz*0.95*dimz);
      }
   
  }
  else {
    if (firstflag) {
      probeLimits();
      
      background(255);
      firstflag=false;
    }
    else{
      addpoints();
    }
  }
}

void addpoints(){
  Passenger[] transients = new Passenger[maxnum];
  for (int i=0;i<maxnum;i++) {
    transients[i] = new Passenger(random(xneg,xpos),random(yneg,ypos),random(zneg,zpos),false);
  }
  for (int count=0;count<maxcounts;count++){
    for (int i=0;i<maxnum;i++) {
      transients[i].update();
      stroke(transients[i].r,transients[i].g, transients[i].b,10);
      point(transients[i].x/gsx*dimx*0.95,transients[i].y/gsy*dimy*0.95,transients[i].z/gsz*dimz*0.95);
    }
  }
}

void reset() {
  counter=0;
    background(bkg);
  

  for (int i=0; i<6;i++) xPdJ[i] = random(pmin,pmax);
  while(Lyapunov(xPdJ)<0.005){
    for (int i=0; i<6;i++) xPdJ[i] = random(pmin,pmax);
  }
  probeLimits();
  for (int i=0;i<maxnum;i++) {
    passengers[i] = new Passenger(random(xneg,xpos),random(yneg,ypos),random(zneg,zpos),true);
  }
  firstflag=true;
}

float Lyapunov(float[] par){

  float x1=random(-0.01,0.01);
  float y1=random(-0.01,0.01);
  float z1=random(-0.01,0.01);
  float x2=x1+1e-6;
  float y2=y1;
  float z2=z1;
  float xt, yt, zt, factor;
  float tmpL=0.0;
  float d=1e-6;
  for(int i=0;i<10000;i++){
    xt = sin(par[0]*y1) - cos(par[3]*z1);
    yt = sin(par[1]*x1) - cos(par[4]*z1);
    zt = sin(par[2]*x1) - cos(par[5]*y1);
    x1=xt;
    y1=yt;
    z1=zt;
    xt = sin(par[0]*y2) - cos(par[3]*z2);
    yt = sin(par[1]*x2) - cos(par[4]*z2);
    zt = sin(par[2]*x2) - cos(par[5]*y2);
    x2=xt;
    y2=yt;
    z2=zt;
    tmpL+=(float)Math.log(1e6*dist(x1,y1,z1,x2,y2,z2));

    if (d<1e-12){
      x2=x1+1e-6;
      y2=y1;
      z2=z1;
    }
    else{
      x2=x1+(x2-x1)/d*1e-6;
      y2=y1+(y2-y1)/d*1e-6;
      z2=z1+(z2-z1)/d*1e-6;
    }
  }
  return(tmpL/69314.718055994530941723212145818);
}



void probeLimits() {
  float xmin=10.0;
  float xmax=-10.0;
  float ymin=10.0;
  float ymax=-10.0;
  float zmin=10.0;
  float zmax=-10.0;
  float xavg=0.0;
  float yavg=0.0;
  float zavg=0.0;
  for(int i = 0; i<maxcounts*maxnum; i++){
    xmin=min(xmin,points[i][0]);
    ymin=min(ymin,points[i][1]);
    xmax=max(xmax,points[i][0]);
    ymax=max(ymax,points[i][1]);
    zmin=min(zmin,points[i][2]);
    zmax=max(zmax,points[i][2]);
    xavg+=points[i][0];
    yavg+=points[i][1];
    zavg+=points[i][2];
  }
  xavg=xavg/(maxcounts*maxnum);
  yavg=yavg/(maxcounts*maxnum);
  zavg=zavg/(maxcounts*maxnum);
  gsx=2*max(xmax-xavg,xavg-xmin);
  gsy=2*max(ymax-yavg,yavg-ymin);
  gsz=2*max(zmax-zavg,zavg-zmin);
  gx=-xavg/gsx;
  gy=-yavg/gsy;
  gz=-zavg/gsz;

  
}

class Passenger {
  float x, y, z;
  float xnext, ynext, znext;
  int r,g,b;
  int age;
  boolean remember;
  Passenger() {
    x=random(xneg,xpos);
    y=random(yneg,ypos);
    z=random(zneg,zpos);
    remember=false;
  }

  Passenger(float xx, float yy, float zz, boolean rr) {
    x=xnext=xx;
    y=ynext=yy;
    z=znext=zz;
    remember=rr;
    age=0;
    r=g=b=0;
    if (remember) {
      points[counter][0]=x;
      points[counter][1]=y;
      points[counter][2]=z;
      colors[counter]=color(r,g,b);
      counter++;
    }
  }

  void update() {
    // trace path
    xnext = sin(xPdJ[0]*y) - cos(xPdJ[1]*z);
    ynext = sin(xPdJ[2]*x) - cos(xPdJ[3]*z);
    znext = sin(xPdJ[4]*x) - cos(xPdJ[5]*y); 
    int rold=r;
    r=int(dist(x,y,z,xnext,ynext,znext)/1.404*125);
   
    b=50*int(float(rold)/r);
    g=(r+b)/2;
    x=xnext;
    y=ynext;
    z=znext;
    
    if (remember) {
      points[counter][0]=x;
      points[counter][1]=y;
      points[counter][2]=z;
      colors[counter]=color(r,g,b);
      counter++;
    }
    age++;
    if (age>maxage) {
      restart();
    }
  }

  void restart() {
    x=xnext=random(xneg,xpos);
    y=ynext=random(yneg,ypos);
    z=znext=random(zneg,zpos);
    r=g=b=0;
    if (remember) {
      points[counter][0]=x;
      points[counter][1]=y;
      points[counter][2]=z;
      colors[counter]=color(r,g,b);
      counter++;
    }
    age = 0;
  }

}

// Arcball, related classes and functions by arielm - June 23, 2003
// http://www.chronotext.org
class ArcBall
{
  float center_x, center_y, radius;
  Vec3 v_down, v_drag;
  Quat q_now, q_down, q_drag;
  Vec3[] axisSet;
  int axis;

  ArcBall(float center_x, float center_y, float radius)
  {
    this.center_x = center_x;
    this.center_y = center_y;
    this.radius = radius;

    v_down = new Vec3();
    v_drag = new Vec3();

    q_now = new Quat();
    q_down = new Quat();
    q_drag = new Quat();

  axisSet = new Vec3[] {new Vec3(1.0f, 0.0f, 0.0f), new Vec3(0.0f, 1.0f, 0.0f), new Vec3(0.0f, 0.0f, 1.0f)};
    axis = -1;  // no constraints...
  }

  void mousePressed()
  {
    v_down = mouse_to_sphere(bufferedMouse.currentX, bufferedMouse.currentY);
    q_down.set(q_now);
    q_drag.reset();
  }

  void mouseDragged()
  {
    v_drag = mouse_to_sphere(bufferedMouse.currentX, bufferedMouse.currentY);
    q_drag.set(Vec3.dot(v_down, v_drag), Vec3.cross(v_down, v_drag));
  }

  void run()
  {
    q_now = Quat.mul(q_drag, q_down);
    applyQuat2Matrix(q_now);
  }

  Vec3 mouse_to_sphere(float x, float y)
  {
    Vec3 v = new Vec3();
    v.x = (x - center_x) / radius;
    v.y = (y - center_y) / radius;

    float mag = v.x * v.x + v.y * v.y;
    if (mag > 1.0f)
    {
      v.normalize();
    }
    else
    {
      v.z = sqrt(1.0f - mag);
    }

    return (axis == -1) ? v : constrain_vector(v, axisSet[axis]);
  }

  Vec3 constrain_vector(Vec3 vector, Vec3 axis)
  {
    Vec3 res = new Vec3();
    res.sub(vector, Vec3.mul(axis, Vec3.dot(axis, vector)));
    res.normalize();
    return res;
  }

  void applyQuat2Matrix(Quat q)
  {
    // instead of transforming q into a matrix and applying it...

    float[] aa = q.getValue();
    rotate(aa[0], aa[1], aa[2], aa[3]);
  }
}

static class Vec3
{
  float x, y, z;

  Vec3()
  {
  }

  Vec3(float x, float y, float z)
  {
    this.x = x;
    this.y = y;
    this.z = z;
  }

  void normalize()
  {
    float length = length();
    x /= length;
    y /= length;
    z /= length;
  }

  float length()
  {
    return (float) Math.sqrt(x * x + y * y + z * z);
  }

  static Vec3 cross(Vec3 v1, Vec3 v2)
  {
    Vec3 res = new Vec3();
    res.x = v1.y * v2.z - v1.z * v2.y;
    res.y = v1.z * v2.x - v1.x * v2.z;
    res.z = v1.x * v2.y - v1.y * v2.x;
    return res;
  }

  static float dot(Vec3 v1, Vec3 v2)
  {
    return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
  }

  static Vec3 mul(Vec3 v, float d)
  {
    Vec3 res = new Vec3();
    res.x = v.x * d;
    res.y = v.y * d;
    res.z = v.z * d;
    return res;
  }

  void sub(Vec3 v1, Vec3 v2)
  {
    x = v1.x - v2.x;
    y = v1.y - v2.y;
    z = v1.z - v2.z;
  }

  void add(Vec3 v1, Vec3 v2)
  {
    x = v1.x + v2.x;
    y = v1.y + v2.y;
    z = v1.z + v2.z;
  }
}

static class Quat
{
  float w, x, y, z;

  Quat()
  {
    reset();
  }

  Quat(float w, float x, float y, float z)
  {
    this.w = w;
    this.x = x;
    this.y = y;
    this.z = z;
  }

  void reset()
  {
    w = 1.0f;
    x = 0.0f;
    y = 0.0f;
    z = 0.0f;
  }

  void set(float w, Vec3 v)
  {
    this.w = w;
    x = v.x;
    y = v.y;
    z = v.z;
  }

  void set(Quat q)
  {
    w = q.w;
    x = q.x;
    y = q.y;
    z = q.z;
  }

  static Quat mul(Quat q1, Quat q2)
  {
    Quat res = new Quat();
    res.w = q1.w * q2.w - q1.x * q2.x - q1.y * q2.y - q1.z * q2.z;
    res.x = q1.w * q2.x + q1.x * q2.w + q1.y * q2.z - q1.z * q2.y;
    res.y = q1.w * q2.y + q1.y * q2.w + q1.z * q2.x - q1.x * q2.z;
    res.z = q1.w * q2.z + q1.z * q2.w + q1.x * q2.y - q1.y * q2.x;
    return res;
  }

  float[] getValue()
  {
    // transforming this quat into an angle and an axis vector...

    float[] res = new float[4];

    float sa = (float) Math.sqrt(1.0f - w * w);
    if (sa < EPSILON)
    {
      sa = 1.0f;
    }

    res[0] = (float) Math.acos(w) * 2.0f;
    res[1] = x / sa;
    res[2] = y / sa;
    res[3] = z / sa;

    return res;
  }
}

class BufferedMouse{
  float stiffness;
  float currentX, currentY;

  BufferedMouse(float ss){
    if (ss<0.0) ss=0.0;
    if (ss>1.0) ss=1.0;
    stiffness=ss;
    currentX=mouseX;
    currentY=mouseY;
  }

  void update(){
    currentX=currentX*stiffness+(1.0-stiffness)*mouseX;
    currentY=currentY*stiffness+(1.0-stiffness)*mouseY;
  }

  void clear(){
    currentX=mouseX;
    currentY=mouseY;
  }

}