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

class Sphere{
  double c[];//transformed center
  double o[];//center
  double v[];
  double f[];
  double r;//radius
  double m;
  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;
  public Sphere(double x, double y, double z, double r){
    o= new double[] {x,z,y};
    c= new double[] {x,z,y};
    this.r  = r;
    v= new double[] {0,0,0};
    f= new double[] {0,0,0};
    m=1;
  }
  public Sphere metal(){
    specularColor[0] = diffuseColor[0];
    specularColor[1] = diffuseColor[1];
    specularColor[2] = diffuseColor[2];
    HW10.normalize(specularColor);
    specularColor[0] *= 0.51;
    specularColor[1] *= 0.51;
    specularColor[2] *= 0.51;
    return this;
  }
  
  public Sphere blue(){
    ambientColor = new double[] {0.1,0.1,0.2};
    diffuseColor = new double[] {0.1,0.1,0.8};
    return this;
  }
  public Sphere green(){
    ambientColor = new double[] {0.1,0.2,0.1};
    diffuseColor = new double[] {0.1,0.8,0.1};
    return this;
  }
  public Sphere pink(){
    ambientColor = new double[] {0.2,0.1,0.1};
    diffuseColor = new double[] {0.6,0.0,0.3};
    return this;
  }
  public Sphere purple(){
    ambientColor = new double[] {0.12,0.1,0.18};
    diffuseColor = new double[] {0.3,0.0,0.6};
    return this;
  }
  public Sphere 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 HW10 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];

  double grav[] = new double[] {0,-.4,0};//gravity force

  double depth;
  
  int res = 10;

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

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

  boolean small = false;
  
  double oldTime;
  double dt;
  LinkedList<Sphere> sphere_list;
  final double EPSILON = 0.0001; 
  final double INFINITY = Double.MAX_VALUE;
  double eye[]; 
  Sphere moon;
  Sphere head;
  Sphere leye;
  Sphere reye;
  Sphere lpupil;
  Sphere rpupil;
  
  double camRotX;
  double camRotY;
  double camRotZ;
  double rotY;
  double lrotY;
  double rrotY;
  double lrotX;
  double rrotX;
  double r;//radius of the spheres
  double xMAX=2;
  double yMAX=2;
  double zMAX=2;
  boolean mouseDown = false;
  boolean repack = false;
  int pixel_width_max = 10;
  int pixel_width = pixel_width_max;

  
  /**  Do this once per run */
  public void initialize() {
    r = 0.56;
    repack  = true;
  }

  public void initFrame(double time){
    if(repack){
      packSpheres();
      repack = false;
      pixel_width = pixel_width_max;
    }
  }
  
  public void packSpheres(){
    double yRow, xRow, zRow;
    sphere_list = new LinkedList<Sphere>();
	boolean isAPlane = true;// used as a switch
	double zAt = EPSILON+r-2*Math.sqrt(6)*2*r/3;//the first ball's z-coord
	double rPrime = (Math.sqrt(6)/2-1)*r;
	while(zAt<=zMAX)//so, when z is less than the z lid on the box...
	{   
	  /*There are two types of planes of balls A and B
	   */   
	  if(isAPlane){
		yRow=0;
		double xAt=-2*r;
		double yAt=r-2*Math.sqrt(3)*r;
		while(yAt<=yMAX){
		  while(xAt<=xMAX){   
			if(xAt<=xMAX-r
			  &&xAt>=r
			  &&yAt<=yMAX-r
			  &&yAt>=r
			  &&zAt<=zMAX-r
			  &&zAt>=r
			){
			  sphere_list.add(new Sphere(xAt, yAt, zAt, r).blue().metal());
			}
			double zSmall = zAt -(rPrime +r);
			if(xAt<=(xMAX-rPrime)
			  &&xAt>=rPrime
			  &&yAt<=(yMAX-rPrime)
			  &&yAt>=rPrime
			  &&zSmall<=(zMAX-rPrime)
			  &&zSmall>=rPrime)
			{
			  if(small) sphere_list.add(new Sphere(xAt, yAt, zSmall, rPrime).green());
			}
			zSmall = zAt +( rPrime +r);
			if(xAt<=(xMAX-rPrime)
			  &&xAt>=rPrime
			  &&yAt<=(yMAX-rPrime)
			  &&yAt>=rPrime
			  &&zSmall<=(zMAX-rPrime)
			  &&zSmall>=rPrime)
			{
			  if(small) sphere_list.add(new Sphere(xAt, yAt, zSmall, rPrime).green());
			}
			xAt = xAt + 2*r;//x coords are always off by 2*r

		  }
		  yRow++;
		  yAt = yAt + Math.sqrt(3)*r;//y coords are always off by this
		  /*
		   *There are two types of rows of balls
		   *they start at different x coords
		   */
		  if(yRow%2==1){
			xAt = r;
		  }
		  if(yRow%2==0){
			xAt = 2*r;
		  }
		}
	  }
	  else
	  {   
		yRow=0;
		double xAt=-3*r;
		double yAt=Math.sqrt(3)/3*r+r-2*Math.sqrt(3)*r;
		while(yAt<=yMAX){
		  while(xAt<=xMAX){
			if(xAt<=xMAX-r
			  &&xAt>=r
			  &&yAt<=yMAX-r
			  &&yAt>=r
			  &&zAt<=zMAX-r
			  &&zAt>=r)
			{
			  sphere_list.add(new Sphere(xAt, yAt, zAt, r).red().metal());
			}
			double zSmall = zAt -(rPrime +r);
			
			if(xAt<=(xMAX-rPrime)
			   &&xAt>=rPrime
			   &&yAt<=(yMAX-rPrime)
			   &&yAt>=rPrime
			   &&zSmall<=(zMAX-rPrime)
			   &&zSmall>=rPrime)
			{
			  if(small)  sphere_list.add(new Sphere(xAt, yAt, zSmall, rPrime).green());
			}
			zSmall = zAt +( rPrime +r);
			if(   xAt<=(xMAX-rPrime)
			  &&xAt>=rPrime
			  &&yAt<=(yMAX-rPrime)
			  &&yAt>=rPrime
			  &&zSmall<=(zMAX-rPrime)
			  &&zSmall>=rPrime)
			{
			  if(small) sphere_list.add(new Sphere(xAt, yAt, zSmall, rPrime).green());
			}
			xAt = xAt + 2*r;
		  }
		  yRow++;
		  yAt = yAt + Math.sqrt(3)*r;
		  if(yRow%2==0)			{
			xAt = r;
		  }
		  if(yRow%2==1){
			xAt = 2*r;
		  }
		}
	  }
	  isAPlane = !isAPlane;
	  zAt = zAt +Math.sqrt(6.)*(2./3.)*r;//z coords are always off by this
	}
  }
        
  /** Overriden method that must define frame in doube pix[] array */
  public void computeImage(double time){
    initFrame(time);
    dt = Math.min(1,time - oldTime);
    oldTime = time;
    
    Matrix4x4 m = new Matrix4x4();
    m.identity();
    m.translate(-1,-1,-21);
    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]);
    }
    // moon.c[2] = -20+1.5*Math.cos(time);
    //moon.c[0] = 1.5*Math.sin(time);


    double rgb[];
    int color[];
    double v  [] = new double[3];
    double w  [] = new double[3];
    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);
        depth = -1;
        rayTrace(v, w, rgb);          // 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);
        //setPixel(i, j, color);
    }
  }
  public boolean keyUp(Event e, int key){
    switch(key){
      case '-':
      case '_':
        r /= 1.1;
        repack=true;
        break;
      case '+':
      case '=':
        r *= 1.1;
        repack=true;
        break;
      case 's':
      case 'S':
        small = !small;
        repack=true;
        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 updatePhysics(){
    /**ball1.addToForce(grav[0]*ball1.m,grav[1]*ball1.m,grav[2]*ball1.m);
    double T[] = {ball1.c[0]-hinge.c[0], ball1.c[1]-hinge.c[1],ball1.c[2]-hinge.c[2] };
    normalize(T);
    double gdt = dot(grav,T);
    double F[] = {-gdt*T[0], -gdt*T[1], -gdt*T[2]};
    ball1.addToForce(F[0]*ball1.m,F[1]*ball1.m,F[2]*ball1.m);
    ball1.applyForce(dt);
    ball1.applyVelocity(dt);
    double h2b[] = {ball1.c[0]-hinge.c[0], ball1.c[1]-hinge.c[1],ball1.c[2]-hinge.c[2] };
    normalize(h2b);
    ball1.c[0] = 5*h2b[0]+hinge.c[0];
    ball1.c[1] = 5*h2b[1]+hinge.c[1];
    ball1.c[2] = 5*h2b[2]+hinge.c[2];
    */
  }
  
  private void rayTrace(double v[], double w[], double rgb[]){
    depth++;
    double t = INFINITY; 
    Sphere 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);
      //System.out.println("tn: "+tn );
      if(tn>0&&tn<t){
        t = tn;
        closest_sphere = sphere;
      }
    }
    if(null!=closest_sphere)
      computeShade(v,w,closest_sphere,t,rgb);
  }

  private void computeShade(
    double v[],
    double w[],
    Sphere sphere,
    double t,
    double rgb[]
  ){
    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
    };
    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);
        v = new double[] {
          point[0],
          point[1],
          point[2]
        };
        for(Sphere s : sphere_list)
          if(s!=sphere){
            double tn = raySphereIntersect(v,L,s);
            if(tn< INFINITY&&tn>0){
              useLight[k] = false;
              break;
            }
          }
      }
      
    }
    ambientColor = sphere.ambientColor;
    diffuseColor = sphere.diffuseColor;
    specularColor = sphere.specularColor;
    specularPower= sphere.specularPower;
    phong(normal, rgb);
    v = new double[] {
      point[0] + EPSILON*normal[0],
      point[1] + EPSILON*normal[1],
      point[2] + EPSILON*normal[2]
    };
    double wdn = dot(w,normal);// w dot normal
    w = 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(v,w,rgb);  
  }

  private double raySphereIntersect(double v[], double w[], Sphere sphere){
    //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;
    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[]){
    // 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] += 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;
}