public class Matrix3D {
  
  private float entries[][] = new float[4][4];
  
  //
  // Constructors
  //

  public Matrix3D() {                 // initialize with identity transform
    loadIdentity();
  }
  
  public Matrix3D(Matrix3D copy) {    // initialize with copy of source
    int i,j;
    for(i=0; i<4; i++) {
      for(j=0; j<4; j++) {
	entries[i][j] = copy.get(j, i);
      }
    }
  }

  public Matrix3D(Raster r) {         // initialize with a mapping from canonical space to screen space
    
    float width = r.getWidth();
    float height = r.getHeight();
    int right, left, top, bottom, near, far;
    right = bottom = near = 1;
    left = top = far = -1;
    int zmax = 1;
    
    int i,j;
    for(j=0; j<4; j++)
      for(i=0; i<4; i++)
	entries[i][j] = 0;    
    
    entries[0][0] = width/(right-left);
    entries[0][3] = -left*width/(right-left);
    entries[1][1] = height/(bottom-top);
    entries[1][3] = -top*height/(bottom-top);
    entries[2][2] = zmax/(far-near);
    entries[2][3] = -near*zmax/(far-near);
    entries[3][3] = 1;
    
  }

  
  //
  // General interface methods
  //        
  public void set(int i, int j, float value) {         // set element [j][i] to value
    entries[j][i] = value;
  }

  public float get(int i, int j) {                     // return element [j][i]
    return(entries[j][i]);
  }
  
  //
  // Transform points from the in array to the out array
  // using the current matrix. The subset of points transformed
  // begins at the start index and has the specified length
  //
  //    for (i = 0; i < length; i++)
  //        out[start+i] = this * in[start+i]
  //
  public void transform(Point3D in[], Point3D out[], int start, int length) {
    int i;
    float out_length;

    for (i=0; i<length; i++) {
      out_length = ((this.get(0,3) * in[start+i].x) +
		    (this.get(1,3) * in[start+i].y) +
		    (this.get(2,3) * in[start+i].z) + 
		    (this.get(3,3) * 1));
      
      out[start+i].x = ((this.get(0,0) * in[start+i].x) +
			(this.get(1,0) * in[start+i].y) +
			(this.get(2,0) * in[start+i].z) + 
			(this.get(3,0) * 1)) / out_length;
      out[start+i].y = ((this.get(0,1) * in[start+i].x) +
			(this.get(1,1) * in[start+i].y) +
			(this.get(2,1) * in[start+i].z) + 
			(this.get(3,1) * 1)) / out_length;
      out[start+i].z = ((this.get(0,2) * in[start+i].x) +
			(this.get(1,2) * in[start+i].y) +
			(this.get(2,2) * in[start+i].z) + 
			(this.get(3,2) * 1)) / out_length;
    }
  }
  
  public final void compose(Matrix3D src) {                 // this = this * src
    int i,j,k;
    float new_entries[][] = new float[4][4];
    for(i=0; i<4; i++) {
      for(j=0; j<4; j++) {
	new_entries[i][j] = 0;
	for(k=0; k<4; k++) {
	  new_entries[i][j] += entries[i][k] * src.entries[k][j];
	}
      }
    }
    for(i=0; i<4; i++) {
      for(j=0; j<4; j++) {
	entries[i][j] = new_entries[i][j];
      }
    }
  }
  
  public void loadIdentity() {                                  // this = identity
    int i,j;
    for(i=0; i<4; i++) {
      for(j=0; j<4; j++) {
	entries[i][j] = (i==j) ? 1 : 0;
      }
    }
  }

  public void add(Matrix3D m) {                                // this = this + m
    int i,j;
    for(j=0; j<4; j++)
      for(i=0; i<4; i++)
	this.set(i,j,m.get(i,j));
  }
  
  public void translate(float tx, float ty, float tz) {        // this = this * t
    Matrix3D tr_matrix = new Matrix3D();
    tr_matrix.set(3, 0, tx);
    tr_matrix.set(3, 1, ty);
    tr_matrix.set(3, 2, tz);
    this.compose(tr_matrix);
  }

  public void scale(float sx, float sy, float sz) {         // this = this * scale
    Matrix3D sc_matrix = new Matrix3D();
    sc_matrix.set(0, 0, sx);
    sc_matrix.set(1, 1, sy);
    sc_matrix.set(2, 2, sz);    
    this.compose(sc_matrix);
  }

  public void symmetric(float ax, float ay, float az) { // this = this * symmetric
    Matrix3D sym_matrix = new Matrix3D();
    sym_matrix.set(0, 0, ax*ax);
    sym_matrix.set(1, 0, ax*ay);
    sym_matrix.set(2, 0, ax*az);
    sym_matrix.set(0, 1, ax*ay);
    sym_matrix.set(1, 1, ay*ay);
    sym_matrix.set(2, 1, ay*az);
    sym_matrix.set(0, 2, ax*az);
    sym_matrix.set(1, 2, ay*az);
    sym_matrix.set(2, 2, az*az);
    this.compose(sym_matrix);
  }

  public void skew(float ax, float ay, float az) {          // this = this * skew
    Matrix3D sk_matrix = new Matrix3D();
    sk_matrix.set(0, 0, 0);
    sk_matrix.set(1, 0, -az);
    sk_matrix.set(2, 0, ay);
    sk_matrix.set(0, 1, az);
    sk_matrix.set(1, 1, 0);
    sk_matrix.set(2, 1, -ax);
    sk_matrix.set(0, 2, -ay);
    sk_matrix.set(1, 2, ax);
    sk_matrix.set(2, 2, 0);
    this.compose(sk_matrix);
  }
  
  public void shear(float kxy, float kxz, float kyz) {     // this = this * shear
    Matrix3D sh_matrix = new Matrix3D();
    sh_matrix.set(1, 0, kxy);
    sh_matrix.set(2, 0, kxz);
    sh_matrix.set(2, 1, kyz);
    this.compose(sh_matrix);
  }
  
  public void scalar_mult(float a) {                  // this = a * this
    int i,j;
    for(j=0; j<4; j++)
      for(i=0; i<4; i++)
	this.set(i,j, a * this.get(i,j));
  }
  
  public void rotate(float ax, float ay, float az, float angle) { // this = this * rotate
    Matrix3D ro_matrix = new Matrix3D();

    float a_length = (float) (Math.sqrt(ax*ax + ay*ay + az*az));
    ax = ax/a_length;
    ay = ay/a_length;
    az = az/a_length;

    float a, b, c;
    a = (float) (1-Math.cos(angle));
    b = (float) (Math.sin(angle));
    c = (float) (Math.cos(angle));
    
    ro_matrix.set(0,0, a*ax*ax + c);
    ro_matrix.set(1,0, a*ax*ay - b*az);
    ro_matrix.set(2,0, a*ax*az + b*ay);
    ro_matrix.set(0,1, a*ax*ay + b*az);
    ro_matrix.set(1,1, a*ay*ay + c);
    ro_matrix.set(2,1, a*ay*az - b*ax);
    ro_matrix.set(0,2, a*ax*az - b*ay);
    ro_matrix.set(1,2, a*ay*az + b*ax);
    ro_matrix.set(2,2, a*az*az + c);
    

    /*
      ro_matrix.set(0,0,0);
      ro_matrix.set(1,1,0);
      ro_matrix.set(2,2,0);
      ro_matrix.set(3,3,0);
      
      Matrix3D sym_matrix = new Matrix3D();
      Matrix3D sk_matrix = new Matrix3D();
      Matrix3D id_matrix = new Matrix3D();
      
      sym_matrix.symmetric(ax, ay, az);
      sym_matrix.scalar_mult((float) (1-Math.cos(angle)));
      
      sk_matrix.skew(ax, ay, az);
      sk_matrix.scalar_mult((float) (Math.sin(angle)));
      
      id_matrix.scalar_mult((float) (Math.cos(angle)));
      
      ro_matrix.add(sym_matrix);
      ro_matrix.add(sk_matrix);
      ro_matrix.add(id_matrix);
      ro_matrix.set(3,3,1);
      */

    this.compose(ro_matrix);
  }
  
  public void lookAt(float eyex, float eyey, float eyez,
		     float atx,  float aty,  float atz,
		     float upx,  float upy,  float upz) {   // this = this * lookat
    
    Matrix3D look_matrix = new Matrix3D();
    float lx,ly,lz,l_length,lxh,lyh,lzh;
    float rx,ry,rz,r_length,rxh,ryh,rzh;
    float ux,uy,uz,u_length,uxh,uyh,uzh;
      
    lx = atx - eyex;
    ly = aty - eyey;
    lz = atz - eyez;
    l_length = (float) Math.sqrt(lx*lx + ly*ly + lz*lz);
    
    lxh = lx/l_length;
    lyh = ly/l_length;
    lzh = lz/l_length;
    
    // r = l x up
    rx = -lz*upy + ly*upz;
    ry = lz*upx - lx*upz;
    rz = -ly*upx + lx*upy;
    r_length = (float) Math.sqrt(rx*rx + ry*ry + rz*rz);
    
    rxh = rx/r_length;
    ryh = ry/r_length;
    rzh = rz/r_length;
    
    // u = r x l
    ux = -rz*ly + ry*lz;
    uy = rz*lx - rx*lz;
    uz = -ry*lx + rx*ly;
    u_length = (float) Math.sqrt(ux*ux + uy*uy + uz*uz);
    
    uxh = ux/u_length;
    uyh = uy/u_length;
    uzh = uz/u_length;
    
    look_matrix.set(0,0,rxh);
    look_matrix.set(1,0,ryh);
    look_matrix.set(2,0,rzh);
    look_matrix.set(3,0,-(rxh*eyex + ryh*eyey + rzh*eyez));
    look_matrix.set(0,1,uxh);
    look_matrix.set(1,1,uyh);
    look_matrix.set(2,1,uzh);
    look_matrix.set(3,1,-(uxh*eyex + uyh*eyey + uzh*eyez));
    look_matrix.set(0,2,-lxh);
    look_matrix.set(1,2,-lyh);
    look_matrix.set(2,2,-lzh);
    look_matrix.set(3,2,(lxh*eyex + lyh*eyey + lzh*eyez));
    look_matrix.set(0,3,0);
    look_matrix.set(1,3,0);
    look_matrix.set(2,3,0);
    look_matrix.set(3,3,1);
    
    this.compose(look_matrix);
  }

  
  //
  // Assume the following projection transformations
  // transform points into the canonical viewing space
  //
  
  public void perspective(float left, float right,           // this = this * persp
			  float bottom, float top,
			  float near, float far) {
    Matrix3D per_matrix = new Matrix3D();
    per_matrix.set(0, 0, 2*near/(right-left));
    per_matrix.set(2, 0, -(right+left)/(right-left));
    per_matrix.set(1, 1, 2*near/(bottom-top));
    per_matrix.set(2, 1, -(bottom+top)/(bottom-top));
    per_matrix.set(2, 2, (far+near)/(far-near));
    per_matrix.set(3, 2, -2*far*near/(far-near));
    per_matrix.set(2, 3, 1);
    per_matrix.set(3, 3, 0);
    this.compose(per_matrix);
  }
  
  public void orthographic(float left, float right,         // this = this * ortho
			   float bottom, float top,
			   float near, float far) {
    Matrix3D orth_matrix = new Matrix3D();
    orth_matrix.set(0, 0, 2/(right-left));
    orth_matrix.set(3, 0, -(right+left)/(right-left));
    orth_matrix.set(1, 1, 2/(bottom-top));
    orth_matrix.set(3, 1, -(bottom+top)/(bottom-top));
    orth_matrix.set(2, 2, 2/(far-near));
    orth_matrix.set(3, 2, -(far+near)/(far-near));
    this.compose(orth_matrix);
  }

}