import java.awt.*;
import java.util.LinkedList;

class Sphere extends ImplicitSurface{
  public Sphere(double x, double y, double z, double r){
    o= new double[] {x,y,z};
    c= new double[] {x,y,z};
    this.r  = r;
    v= new double[] {0,0,0};
    f= new double[] {0,0,0};
    m=1;
    contact_normal = null;
  }
}

class Cube extends ImplicitSurface{
  double plane[][] = {
    {-1, 0, 0,-1},
    { 1, 0, 0,-1},
    { 0,-1, 0,-1},
    { 0, 1, 0,-1},
    { 0, 0,-1,-21},
    { 0, 0, 1,-19}
  };
  public Cube(double x, double y, double z, double r){
    o= new double[] {x,y,z};
    c= new double[] {x,y,z};
    this.r  = r;
    v= new double[] {0,0,0};
    f= new double[] {0,0,0};
    m=1;
  }
}

class ImplicitSurface{
  double c[];//transformed center
  double o[];//center
  double v[];
  double f[];
  double r;//radius
  double n = 1.1;//water
  double m;
  int type;
  double ambientColor[] = { 0.2, 0.2, 0.2  };
  double diffuseColor[] = { 0.8, 0.8, 0.8 };
  double specularColor[] = {.1,.1,.1};
  double specularPower = 3.0;
  double contact_normal[];
  double opacity = 1.0;
  
  public ImplicitSurface metal(){
    specularColor[0] = diffuseColor[0];
    specularColor[1] = diffuseColor[1];
    specularColor[2] = diffuseColor[2];
    HW11.normalize(specularColor);
    specularColor[0] *= 0.51;
    specularColor[1] *= 0.51;
    specularColor[2] *= 0.51;
    return this;
  }
  
  public ImplicitSurface blue(){
    ambientColor = new double[] {0.1,0.1,0.2};
    diffuseColor = new double[] {0.1,0.1,0.8};
    return this;
  }
  public ImplicitSurface  green(){
    ambientColor = new double[] {0.1,0.2,0.1};
    diffuseColor = new double[] {0.1,0.8,0.1};
    return this;
  }
  public ImplicitSurface pink(){
    ambientColor = new double[] {0.2,0.1,0.1};
    diffuseColor = new double[] {0.6,0.0,0.3};
    return this;
  }
  public ImplicitSurface purple(){
    ambientColor = new double[] {0.12,0.1,0.18};
    diffuseColor = new double[] {0.3,0.0,0.6};
    return this;
  }
  public ImplicitSurface indigo(){
    ambientColor = new double[] {0.1,0.12,0.18};
    diffuseColor = new double[] {0.0,0.3,0.6};
    return this;
  }
  public ImplicitSurface yellow(){
    ambientColor = new double[] {0.18,0.12,0.1};
    diffuseColor = new double[] {0.4,0.4,0.0};
    return this;
  }
  public ImplicitSurface clear(){
    opacity = 0.0;
    return this;
  }
  public ImplicitSurface semiclear(){
    opacity = 0.5;
    return this;
  }
  
  public ImplicitSurface red(){
    ambientColor = new double[] {0.2,0.1,0.1};
    diffuseColor = new double[] {0.8,0.1,0.1};
    return this;
  }
  
  public void applyForce(double dt){
    v[0] += f[0]/m*dt;
    v[1] += f[1]/m*dt;
    v[2] += f[2]/m*dt;
  }
  public void applyVelocity(double dt){
    c[0] += v[0]*dt;
    c[1] += v[1]*dt;
    c[2] += v[2]*dt;
  }
  public void addToForce(double x, double y, double z){
     f[0]+=x;
     f[1]+=y;
     f[2]+=z;
  }
}

public class HW11 extends MISApplet{
  // MATERIAL PROPERTIES
  double ambientColor[] = { 0.2, 0., 0.0  };
  double diffuseColor[] = { 0.8, 0.0, 0.0 };
  double specularColor[] = {.1,.1,.1};
  double specularPower = 20.0;
  // LIGHTING PROPERTIES
  /**
  double lightDirections[][] = { {1,0,0} };
  double lightColors[][]     = { {1.0,1.0,1.0} };
  */
  double lightDirections[][] = { {1,1,1}, {0,0,1} ,{-1,0,-1}   };
  double lightColors[][]     = { {0.8,0.9,0.8}, {0.5,0.6,0.5}, {0.5,0.6,0.5} };
  boolean useLight[] = new boolean[lightDirections.length];


  int res = 10;

  double foc = 5.0;//focal length of pinhole camera

  int type;
  double vec1[] = new double[3];//temporary 3D vector
  double vec2[] = new double[3];//temporary 3D vector

  double oldTime;
  double dt;
  LinkedList<Sphere> sphere_list;
  LinkedList<Cube> cube_list;
  final double EPSILON = 0.0001; 
  final double INFINITY = Double.MAX_VALUE;
  double eye[]; 
  boolean mouseDown = false;
  boolean texture = false;
  int pixel_width_max = 10;
  int pixel_width = pixel_width_max;
  int pass_pixel_width = pixel_width;
  double camRotX;
  double camRotY;
  double camRotZ;
  double rotY;
  double lrotY;
  double rrotY;
  double lrotX;
  double rrotX;
  double xMAX=2;
  double yMAX=2;
  double zMAX=2;
  
  
  /**  Do this once per run */
  public void initialize() {
    cube_list = new LinkedList<Cube>();
    sphere_list = new LinkedList<Sphere>();
    Cube c = new Cube(0,0,2,0.5);
    c.yellow().metal();
    //cube_list.add(c);
    Sphere s = new Sphere(0,0,0,0.5);
    s.indigo().metal();
    s.type = 2;
    sphere_list.add(s);
    Sphere u = new Sphere(0,0,1,0.25);
    u.purple().metal();
    u.type = 1;
    sphere_list.add(u);
    Sphere t = new Sphere(0,1.0,-1.5,1.0);
    t.yellow().semiclear().metal();
    t.type = 0;
    sphere_list.add(t);
  }

  
  /** Overriden method that must define frame in doube pix[] array */
  public void computeImage(double time){
    dt = Math.min(1,time - oldTime);
    oldTime = time;
    
    Matrix4x4 m = new Matrix4x4();
    m.identity();
    m.translate(0,0,-20);
    m.rotateX(camRotX);
    m.rotateY(camRotY);
    m.rotateZ(camRotZ);
    for(Sphere s: sphere_list){
      m.translate(s.o[0],s.o[1],s.o[2]);
      m.transform(0,0,0,s.c);
      m.translate(-s.o[0],-s.o[1],-s.o[2]);
    }
    for(Cube s: cube_list){
      m.translate(s.o[0],s.o[1],s.o[2]);
      m.transform(0,0,0,s.c);
      m.translate(-s.o[0],-s.o[1],-s.o[2]); 
    }

    double rgb[];
    int color[];
    double v  [] = new double[3];
    double w  [] = new double[3];
    if(pass_pixel_width > 1)
      pixel_width = pass_pixel_width;
    pass_pixel_width = -1;
    if(mouseDown)
      pixel_width = pixel_width_max;
    else
      pixel_width = Math.max(1, pixel_width/2);
    
    for (int i = 0 ; i < W ; i+=pixel_width)   // LOOP OVER IMAGE COLUMNS
    for (int j = 0 ; j < H ; j+=pixel_width) {  // LOOP OVER IMAGE ROWS
        //System.out.println("ij "+i+" , "+j);
        v[0] = 0; 
        v[1] = 0;
        v[2] = 0;
        
        w[0] = (double)(i - W/2) / W; // COMPUTE RAY DIRECTION AT EACH PIXEL
        w[1] = (double)(H/2 - j) / W; //
        w[2] = -foc;          // PLACE IMAGE PLANE AT z=-focalLength

        rgb = new double[] {0,0,0};//initialize pixel color to black
        normalize(w);
        rayTrace(v, w, rgb, 0, null , false);          // COMPUTE COLOR AT PIXEL BY RAY TRACING
        color = new int[] {(int)(255*rgb[0]), (int)(255*rgb[1]), (int)(255*rgb[2]) };
        for(int _i = i; _i<W && _i < i+pixel_width; _i++)
        for(int _j = j; _j<H && _j < j+pixel_width; _j++)
          pix[xy2i(_i,_j)] = pack(color[0],color[1],color[2]);
        //setPixel(i, j, color);
    }
  }
  
  public boolean keyUp(Event e, int key){
    switch(key){
      case '-':
      case '_':
        pixel_width_max= Math.max(1,pixel_width_max/2);
        pass_pixel_width = pixel_width_max;
        break;
      case '+':
      case '=':
        pixel_width_max*=2;
        pass_pixel_width = pixel_width_max;
        break;
      case 't':
      case 'T':
        texture = !texture;
        pass_pixel_width = pixel_width_max;
        break;
    }
    return true;
  }
  
  public boolean mouseDown(Event e, int x, int y){
    mouseDown = true;
    return true;
  }
  public boolean mouseUp(Event e, int x, int y){
    mouseDown = false;
    return true;
  }
  public boolean mouseDrag(Event e, int x, int y) {
    int wx = W/2;
    int wy = H/2;
    y = H - y - 1;
    double radius = Math.min( W, H )/2.;
    double sx = (x-wx)/radius;
    double sy = (y-wy)/radius;
    if(sx*sx+sy*sy>1){
      double denom = Math.sqrt(sx*sx+sy*sy);
      sx /= denom;
      sy /= denom;
    }
    double sz= Math.sqrt(Math.max(0,1-sx*sx-sy*sy));
    double off[] = {sx,sy,sz};
    normalize(off);
    //camRotX = sign(sy)*Math.acos(dot(off,new double[]{1,0,0}));
    //camRotY = sign(sx)*Math.acos(dot(off,new double[]{0,0,1}));
    //camRotZ = Math.acos(dot(off,new double[]{1,0,0}));
    camRotY = Math.PI*(x-W/2)/W;
    camRotX = -Math.PI*(y-H/2)/H;
    //camRotY = sign(x)*Math.PI/2*(x-W/2)/W;
    //camRotX = -sign(y)*Math.PI/2*(y-H/2)/H;
    return true;
  }
  
  private void rayTrace(double v[], double w[], double rgb[], int depth, ImplicitSurface within, boolean skip){
    double t = INFINITY; 
    ImplicitSurface closest_sphere = null;
    //loop over all spheres
    for(Sphere sphere: sphere_list){
      //Try to intersect the ray with Sphn, to get tn;
      //If the ray intersects Sphn at some tn, and tn < t, then this is the
      //nearest sphere, so:
      double tn = raySphereIntersect(v,w,sphere, skip);
      //System.out.println("tn: "+tn );
      if(tn>0&&tn<t){
        t = tn;
        closest_sphere = sphere;
      }
    }
    for(Cube cube: cube_list){
      double tn = rayCubeIntersect(v,w,cube);
      if(tn>0&&tn<t){
        t = tn;
        closest_sphere = cube;
      }
    }
    if(null!=closest_sphere&&depth<20)
      computeShade(v,w,closest_sphere,t,rgb, depth, within, skip);
  }

  private void computeShade(
    double v[],
    double w[],
    ImplicitSurface sphere,
    double t,
    double rgb[],
    int depth,
    ImplicitSurface within,
    boolean skip
  ){
    double point[] = new double[]{ 
      (v[0]+w[0]*t),
      (v[1]+w[1]*t),
      (v[2]+w[2]*t)
    };
    double normal[] = new double[] {
      (point[0]-sphere.c[0])/sphere.r,
      (point[1]-sphere.c[1])/sphere.r,
      (point[2]-sphere.c[2])/sphere.r
    };
    //this could be so much cleaner..
    if(null!=sphere.contact_normal){
      normal = sphere.contact_normal;
    }
    for(int k = 0;k<useLight.length;k++)
      useLight[k] = true;
    //first find out which lights are in shadow.
    for(int k=0;k<lightDirections.length;k++){
      if(useLight[k]){
        double L[] = new double[] {
           lightDirections[k][0],
           lightDirections[k][1], 
           lightDirections[k][2]};
        normalize(L);
        double vp[] = new double[] {
          point[0],
          point[1],
          point[2]
        };
        // SHADOWS
        for(Sphere s : sphere_list)
          if(s!=sphere){
            double tn = raySphereIntersect(vp,L,s,skip);
            if(tn< INFINITY&&tn>0){
              useLight[k] = false;
              break;
            }
          }
        for(Cube cube: cube_list){
          double tn = rayCubeIntersect(v,w,cube);
          if(tn>0&&tn<INFINITY){
            useLight[k] = false;
            break;
          }
        }
      }
      
    }
    type = sphere.type;
    ambientColor = sphere.ambientColor;
    diffuseColor = sphere.diffuseColor;
    specularColor = sphere.specularColor;
    specularPower= sphere.specularPower;
    if(!skip && sphere.opacity>0)
      phong(normal, rgb, depth, sphere.opacity);

    // REFLECTION
    double vp[] = new double[] {
      point[0] + EPSILON*normal[0],
      point[1] + EPSILON*normal[1],
      point[2] + EPSILON*normal[2]
    };
    normalize(normal);//superfluous?
    double wdn = dot(w,normal);// w dot normal
    double wp[] = new double[] {
      w[0] - 2*wdn*normal[0], 
      w[1] - 2*wdn*normal[1], 
      w[2] - 2*wdn*normal[2] 
    };
    //recursive call to continue ray
    rayTrace(vp,wp,rgb, depth+1,within, false);  
    
    //REFRACTION
    double c[] = new double[] {
      wdn*normal[0], 
      wdn*normal[1], 
      wdn*normal[2] 
    };
    double s[] = new double[] {
      w[0]-c[0], 
      w[1]-c[1], 
      w[2]-c[2] 
    };
    double n1 = null==within ? 1.0 : within.n;
    double n2 = (null!=within && sphere==within) ? 1.0 : sphere.n;
    
    double theta1 = Math.acos(norm(c));
    double theta2 = Math.asin(Math.sin(theta1)*n1/n2);
    normalize(c);
    normalize(s);
    wp = new double[]{
      c[0]*Math.cos(theta2)+s[0]*Math.sin(theta2),  
      c[1]*Math.cos(theta2)+s[1]*Math.sin(theta2),  
      c[2]*Math.cos(theta2)+s[2]*Math.sin(theta2)  
    };
    //normalize(wp);
    vp = new double[] {
      point[0] + EPSILON*wp[0],
      point[1] + EPSILON*wp[1],
      point[2] + EPSILON*wp[2]
    };
    
    //if(theta1<Math.PI/4)
      rayTrace(vp,wp,rgb, depth+1, sphere, sphere==within);  
  }

  private double rayCubeIntersect(double v[], double w[], Cube c){
    double tl = -1;
    double tr = INFINITY;
    for(int i =0; i<c.plane.length; i++){
      double normal[] = new double[]{c.plane[i][0], c.plane[i][1], c.plane[i][2]}; 
      double A = dot(c.plane[i],new double[] {w[0],w[1],w[2],0});
      double B = dot(c.plane[i],new double[] {v[0],v[1],v[2],1});
      if(A!=0){
        double t = -B/A;
        if(c.plane[i][0]*v[0]+c.plane[i][1]*v[1]+c.plane[i][2]*v[2]+c.plane[i][3]<=0){
          //ENTERING
          if(t>tl){
            tl = t;
            c.contact_normal = normal;
          }
        }else{
          //EXITTING
          if(t<tr){
            tr = t;
          }
        }
      }
    }
    return tl;
  }

  private double raySphereIntersect(double v[], double w[], Sphere sphere, boolean skip){
    //vector vmc = v - c 
    double vmc[] = new double[3];
    for(int i =0;i<vmc.length;i++)
      vmc[i] = v[i] - sphere.c[i];
    double A = dot(w,w);
    double B = 2*dot(vmc,w);
    double C = dot(vmc,vmc)-sphere.r*sphere.r;
    double tn = INFINITY;
    if(B*B-4*A*C>=0){
      tn = (-B-Math.sqrt(B*B-4*A*C))/2*A;
      if(skip && tn<=0){
//        System.out.println("INSIDE");
        tn = (-B+Math.sqrt(B*B-4*A*C))/2*A;
      }
    }
    return tn;
  }
  
  /** Phong algorithm for shading vertices
    * takes a normal and an initial color then adds phong shade to color given
    */
  private void phong(double[] normal, double rgb[], int depth, double opacity){
        // SAMPLE THE FUNCTION FOUR TIMES TO GET GRADIENT INFO

         if(texture){
         double f0 = f(normal[0]      ,normal[1]      ,normal[2]      ),
                fx = f(normal[0]+.0001,normal[1]      ,normal[2]      ),
                fy = f(normal[0]      ,normal[1]+.0001,normal[2]      ),
                fz = f(normal[0]      ,normal[1]      ,normal[2]+.0001);

         // SUBTRACT THE FUNCTION'S GRADIENT FROM THE SURFACE NORMAL

         
         normal[0] -= (fx - f0) / .0001;
         normal[1] -= (fy - f0) / .0001;
         normal[2] -= (fz - f0) / .0001;
         normalize(normal);
         }
    // ADD AMBIENT
    for (int i = 0 ; i < 3 ; i++){
      color[i] = 0;
      if(depth==0) color[i] = ambientColor[i];
    }
    // LOOP THROUGH ALL LIGHT SOURCES
    for (int l = 0 ; l < lightDirections.length ; l++) {
      if(useLight[l]){
        // COMPUTE DIFFUSE
        double n_dot_l = dot(lightDirections[l], normal);
        if (n_dot_l > 0)
          for (int i = 0 ; i < 3 ; i++)
            diffuse[i] = diffuseColor[i] * n_dot_l;
        else
          for (int i = 0 ; i < 3 ; i++)
            diffuse[i] = 0;

        // COMPUTE REFLECTION DIRECTION
        for (int i = 0 ; i < 3 ; i++)
          reflection[i] = 2 * n_dot_l * normal[i] - lightDirections[l][i];

        double r_dot_e = dot(reflection, eyeDirection);

        // COMPUTE SPECULAR
        if (r_dot_e > 0) {
          double specularFactor = Math.pow(r_dot_e, specularPower);
          for (int i = 0 ; i < 3 ; i++)
          specular[i] = specularColor[i] * specularFactor;
        }
        else
          for (int i = 0 ; i < 3 ; i++)
            specular[i] = 0;
        
        // ADD DIFFUSE AND SPECULAR
        for (int i = 0 ; i < 3 ; i++)
          color[i] +=  lightColors[l][i] * (diffuse[i] + specular[i]);
      }
    }
      
      
    // SEND COLOR TO THE FRAME BUFFER
    for (int i = 0 ; i < 3 ; i++) {
      rgb[i] += opacity*Math.pow(0.25,depth)* Math.max(0, Math.min(1.0, color[i]));
      rgb[i] = Math.min(rgb[i],1.0);
    }
  }
  
  
  
  /** Universal Merge sort using dist function and point P */
  private double[][] mergeSort(double[][] m){

    //sub-problem arrays and result
    double[][] left, right, result;
    
    //base cases
    if(m.length <= 1)
      return m;

    //find middle and create and fille sub-problem arrays
    int middle = m.length/2;
    left = new double[middle][6];
    right = new double[m.length-middle][6];
    for(int j=0;j<left.length;j++)
      left[j] = m[j];
    for(int j=0;j<right.length;j++)
      right[j] = m[middle+j];

    //recursive calls
    left = mergeSort(left);
    right = mergeSort(right);

    //merge sub-problem results into result
    result = merge(left, right);
    return result;
  }

  private double[][] merge(double[][] left, double[][] right){
    //place for merged result 
    double[][] result = new double[left.length+right.length][6];
    int li = 0;
    int ri = 0;
    //merge using dist to compare
    while(li < left.length && ri < right.length){
      if(left[li][1] <= right[ri][1]){
        result[li+ri] = left[li];
        li++;
      }else{
        result[li+ri] = right[ri];
        ri++;
      }
    }

    //clean up left overs
    while(li<left.length){
      result[li+ri] = left[li];
      li++;
    }    
    while(ri<right.length){
      result[li+ri] = right[ri];
      ri++;
    }
    
    return result;
  }

  // USEFUL MATH
  public static void lerp(double f , double y0[], double y1[], double dst[]){
    for(int i=0;i<y0.length;i++)
      dst[i] = y0[i] + (y1[i]-y0[i])*f;
  }
  public static void lerp(double f , int y0[], int y1[], int dst[]){
    for(int i=0;i<y0.length;i++)
      dst[i] = (int) (y0[i] + (y1[i]-y0[i])*f);
  }
  public static double lerp(double f , double y0, double y1){
    return y0 + (y1-y0)*f;
  }

  public static void normalize(double vec[]) {
     double len = norm(vec);
     for (int i = 0 ; i < 3 ; i++)
        vec[i] /= len; 
  }
  
  /**
   * return z value of surface normal
   * Uses first 3 values
   */
  public static double normalZ(double[][] triangle){
    double[] p1 = triangle[0];
    double[] p2 = triangle[1];
    double[] p3 = triangle[2];
    
    //if rect then I can do a little err check and get a true triangle
    if (triangle.length>3){
      if(equals(p1,p2)){
        p1 = triangle[3];
      }else if(equals(p1,p3)){
        p1 = triangle[3];
      }else if(equals(p2,p3)){
        p3 = triangle[3];
      }
    }
    
    
    double[] a = new double[] {p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]};
    double[] b = new double[] {p3[0]-p1[0],p3[1]-p1[1],p3[2]-p1[2]};
    
    /* Surface normal is (p2-p1)x(p3-p1)
     * of which we only need the z-coordinate
     * so if n = axb then the z-coordinate is 
     * a1b2 - a2b1
     * which in this case is
     * (p21 - p11)(p32-p12)-(p32-p12)(p31-p11)
     */
    //let's try that again
    //(a2b3 - a3b2, a3b1 - a1b3, a1b2 - a2b1)
    double[] cross = new double[] {a[1]*b[2] - a[2]*b[1],
      a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]};
    return cross[2];  
  }
  
  public static boolean equals(double p1[], double p2[]){
    for(int i =0; i<p1.length; i++)
      if(p1[i]!=p2[i]) return false;
    return true;
  }
  
  public static double norm(double vec[]) {
     return Math.sqrt(dot(vec, vec));
  }

  public static double dot(double a[], double b[]) {
     return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
  }
  
  //0 is positive
  public int sign(double d){
    if(d==0) return 0;
    else return (int)(d/Math.abs(d));
  }
  //FINAL CONSTANTS
  final int NUM_SIDES_ON_TRIANGLE = 3;
  final double dampening = 0.99;

  // WORKING VARIABLES
  double color[] = new double[3];
  double diffuse[] = new double[3];
  double specular[] = new double[3];
  double point[] = {0,0,0};
  double normal[] = {0,0,0};
  double eyeDirection[] = {0,0,1};
  double reflection[] = {0,0,0};
  double radius = 0.2;
 // STRIPES TEXTURE (GOOD FOR MAKING MARBLE)

   double stripes(double x, double f) {
      double t = .5 + .5 * Math.sin(f * 2*Math.PI * x);
      return t * t - .5;
   }

   // TURBULENCE TEXTURE

   double turbulence(double x, double y, double z, double freq) {
      double t = -.5;
      for ( ; freq <= W/12 ; freq *= 2)
         t += Math.abs(noise(x,y,z,freq) / freq);
      return t;
   }
        
        // CHOOSE A TYPE OF SPACE FILLING TEXTURE

   double f(double x,double y,double z) {
     switch(type){

      case 0:  return .008 * stripes(x+ 2.2*turbulence(x,y,0,1), 1.6); 
      case 1:  return .02 * stripes(x+ 2.2*turbulence(x,y,0,2), 0.4); 
      default: return .01 * stripes(x+ 2.2*turbulence(x,y,0,1), 0.5); 
       //case 2: return 0.11*turbulence(x+stripes(x,0.5),y,z,0.5); 
       //case 1: return .01 * noise(x,y,z, 10); 
       //default: return 0.12*turbulence(x+stripes(x,0.3),y+stripes(y,0.3),z,1);
     }
   }
   
   // NOISE TEXTURE

   double noise(double x, double y, double z, double freq) {
      double x1, y1, z1;
      x1 = .707*x-.707*z;
      z1 = .707*x+.707*z;
      y1 = .707*x1+.707*y;
      x1 = .707*x1-.707*y;
      return ImprovedNoise.noise(freq*x1 + 100, freq*y1, freq*z1);
   }

}