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

//inspiration http://www.chaoscope.org

tracer[] dots;
int numdots;
float[][] img;
boolean drag;
ArcBall arcball;
BufferedMouse bufferedMouse;
boolean mouseClick=false;

float a1,a2,a3,b1,b2,b3;
float r=1.7;

void setup(){
  size(250,250);
  background(255);
  noFill();
  numdots=10000;
  img = new float[width*height][3];
  for(int i=0;i<width*height;i++){
    img[i][0]=img[i][1]=img[i][2]=0.0;
  }
  dots=new tracer[numdots];
  
  newParameters();
  
  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(){
  translate(width/2,height/2);
  bufferedMouse.update();
  arcball.run();
  for(int i=0;i<numdots;i++){
    dots[i].update();
    
  }
  drawAll();
}

void clearAll(){
  background(0);
  for(int i=0;i<width*height;i++){
    img[i][0]=img[i][1]=img[i][2]=0.0;
  }
}

void drawAll(){
  for(int i=0;i<numdots;i++){
    if (dots[i].count>1){
      float xi=screenX(60.0*dots[i].xc, 60.0*dots[i].yc, 60.0*dots[i].zc);
      float yi=screenY(60.0*dots[i].xc, 60.0*dots[i].yc, 60.0*dots[i].zc);
      int u=int(xi);
      int v=int(yi);
      float dx=xi-u;
      float dy=yi-v;
      if((u>0)&&(v>0)&&(u<width-2)&&(v<height-2)){
        img[u+v*width][0]+=(1.0-dx)*(1.0-dy)*dots[i].r;
        img[u+1+v*width][0]+= dx*(1.0-dy)*dots[i].r;
        img[u+(v+1)*width][0]+= (1.0-dx)*dy*dots[i].r;
        img[u+1+(v+1)*width][0]+= dx*dy*dots[i].r;
        img[u+v*width][1]+=(1.0-dx)*(1.0-dy)*dots[i].g;
        img[u+1+v*width][1]+= dx*(1.0-dy)*dots[i].g;
        img[u+(v+1)*width][1]+= (1.0-dx)*dy*dots[i].g;
        img[u+1+(v+1)*width][1]+= dx*dy*dots[i].g;
        img[u+v*width][2]+=(1.0-dx)*(1.0-dy)*dots[i].b;
        img[u+1+v*width][2]+= dx*(1.0-dy)*dots[i].b;
        img[u+(v+1)*width][2]+= (1.0-dx)*dy*dots[i].b;
        img[u+1+(v+1)*width][2]+= dx*dy*dots[i].b;
      }
    }
  }
  for(int i=0;i<width*height;i++){
    int r=int(255.0*(1.0-Math.exp(-img[i][0]/140.0)));
    int g=int(255.0*(1.0-Math.exp(-img[i][1]/140.0)));
    int b=int(255.0*(1.0-Math.exp(-img[i][2]/140.0)));
    pixels[i]=color(r,g,b);
  }
  if(drag){
  for(int i=0;i<numdots;i+=10){
    dots[i].draw();
  }
  }
}

float Lyapunov(float a1,float a2,float a3,float b1,float b2,float b3,float r){
  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=b1*(z1+a1);
    yt=b2*(x1+a2);
    zt=b3*(y1+a3);
    d=sqrt(xt*xt+yt*yt+zt*zt);
    factor=1.0;
    if (d>abs(r)){
      factor=(2*abs(r)-d)/d;
      xt*=factor;
      yt*=factor;
      zt*=factor;
    }
    x1=xt;
    y1=yt;
    z1=zt;
    xt=b1*(z2+a1);
    yt=b2*(x2+a2);
    zt=b3*(y2+a3);
    d=sqrt(xt*xt+yt*yt+zt*zt);
    factor=1.0;
    if (d>abs(r)){
      factor=(2*abs(r)-d)/d;
      xt*=factor;
      yt*=factor;
      zt*=factor;
    }
    x2=xt;
    y2=yt;
    z2=zt;
    d=dist(x1,y1,z1,x2,y2,z2);
    if(d>0.0) 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 newParameters(){
  clearAll();
  a1=random(-1.0,1.0);
  a2=random(-1.0,1.0);
  a3=random(-1.0,1.0);
  b1=random(1.0,3.0);
  b2=random(1.0,3.0);
  b3=random(1.5,3.0);
  //r=random(0.7,2.0);
  Lyapunov(a1,a2,a3,b1,b2,b3,r);
  while(Lyapunov(a1,a2,a3,b1,b2,b3,r)<0.0005){
    a1=random(-1.0,1.0);
    a2=random(-1.0,1.0);
    a3=random(-1.0,1.0);
    b1=random(1.0,3.0);
    b2=random(1.0,3.0);
    b3=random(1.5,3.0);
    //r=random(0.7,2.0);
  }
  for(int i=0;i<numdots;i++){
    dots[i]=new tracer(a1,a2,a3,b1,b2,b3,r);
  }
}

void keyPressed(){
  newParameters();
}

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

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

void mouseReleased()
{
  mouseClick=false;
  drag=false;
}

class tracer{
  float xc,yc,zc,xn,yn,zn,xp,yp,zp,vel,acc;
  float a1,a2,a3,b1,b2,b3,rad;
  int count;
  float r,g,b;
  tracer(float aa1,float aa2,float aa3,float bb1,float bb2,float bb3,float rr){
    a1=aa1;
    a2=aa2;
    a3=aa3;
    b1=bb1;
    b2=bb2;
    b3=bb3;
    rad=rr;
    r=g=b=1.0;
    xc=random(-1.0,1.0);
    yc=random(-1.0,1.0);
    zc=random(-1.0,1.0);
    float d=sqrt(xc*xc+yc*yc+zc*zc);
    float factor=1.0;
    if(d>abs(rad)){
      factor=abs(rad)/d;
      xc*=factor;
      yc*=factor;
      zc*=factor;
    }
    xp=xn=0;
    yp=yn=0;
    zp=zn=0;
    count=int(random(100));
    vel=acc=1.0;
  }

  void update(){
    xn=b1*(zc+a1);
    yn=b2*(xc+a2);
    zn=b3*(yc+a3);
    float d=sqrt(xn*xn+yn*yn+zn*zn);
    float factor=1.0;
    if (d>abs(rad)){
      factor=(2*abs(rad)-d)/d;
      xn*=factor;
      yn*=factor;
      zn*=factor;
    }
    count++;
    vel=dist(xc,yc,zc,xn,yn,zn)/(2.0*r);
    acc=0.5*((xn-xc)*(xc-xp)+(yn-yc)*(yc-yp)+(zn-zc)*(zc-zp))/(dist(xc,yc,zc,xn,yn,zn)*dist(xc,yc,zc,xp,yp,zp));
    g=constrain(vel,0.0,1.0);
    r=constrain(0.5+acc,0.0,1.0);
    xp=xc;
    yp=yc;
    zp=zc;
    xc=xn;
    yc=yn;
    zc=zn;
    if ((count>1000)&&(vel<0.0001)) reset();
  }

  void reset(){
    xc=random(-1.0,1.0);
    yc=random(-1.0,1.0);
    zc=random(-1.0,1.0);
    float d=sqrt(xc*xc+yc*yc+zc*zc);
    float factor=1.0;
    if(d>abs(rad)){
      factor=(2*abs(rad)-d)/d;
      xc*=factor;
      yc*=factor;
      zc*=factor;
    }
    xp=xn=0;
    yp=yn=0;
    zp=zn=0;
    count=0;    vel=acc=1.0;
  }

  void draw(){
    stroke(255,120);
    point(60.0*xc, 60.0*yc, 60.0*zc);
  }
}


// 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;
  }

}