import java.util.*;

public class Matrix3D {
  private float A[][];

  //
  // Constructors
  //
  public Matrix3D() {                 // initialize with identity transform
    A = new float[4][4];
    loadIdentity();
  }
  public Matrix3D(Matrix3D copy) {    // initialize with copy of source
    A = new float[4][4];
    for(int i=0; i<4; ++i)
      for(int j=0; j<4; ++j)
	A[i][j] = copy.A[i][j];
  }
  public Matrix3D(Raster r) {         // initialize with a mapping from
                                      // canonical space to screen space
    A = new float[4][4];

    Point3D p[] = new Point3D[2];
    Point3D p2[] = new Point3D[2];

    loadIdentity();
    scale(r.getWidth()/2, r.getHeight()/2, 0);
    translate(1, 1, 0);
  }

  public Matrix3D(float value[]) {
    int i, j;

    if (value.length != 16)
      throw new InternalError("value.length != 16");

    A = new float[4][4];

    for(i=0; i<4; ++i)
      for(j=0; j<4; ++j)
	A[i][j] = value[i*4 + j];
  }
  
  //
  // General interface methods
  //        
  public void set(int i, int j, float value) {
    A[i][j] = value;
  }
  public float get(int i, int j) {
    return(A[i][j]);
  }
  
  //
  // 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 transform1(Point3D in, Point3D out) {
    int j;
    float v[] = new float[4];

    out.x = 0;
    out.y = 0;
    out.z = 0;
    float w = 0;
    v[0] = in.x;
    v[1] = in.y;
    v[2] = in.z;
    v[3] = 1;
    for (j=0; j<4; ++j) {
      out.x += A[0][j] * v[j];
      out.y += A[1][j] * v[j];
      out.z += A[2][j] * v[j];
      w += A[3][j] * v[j];
    }
    if (w != 0) {
      out.x /= w;
      out.y /= w;
      out.z /= w;
    }
  }

  public void transform(Point3D in[], Point3D out[], int start, int length) {
    int i;

    for (i = start; i < start+length; ++i)
      transform1(in[i], out[i]);
  }

  
  public final void compose(Matrix3D src) {  // this = this * src
    Matrix3D result = new Matrix3D();

    for(int i=0; i<4; ++i) 
      for(int j=0; j<4; ++j) {
	result.A[i][j] = 0;
	for(int k=0; k<4; ++k)
	  result.A[i][j] += A[i][k] * src.A[k][j];
      }

    for(int i=0; i<4; ++i)
      for(int j=0; j<4; ++j)
	A[i][j] = result.A[i][j];
  }

  public final Matrix3D add(Matrix3D src) {  // this = this + src
    for(int i=0; i<4; ++i) 
      for(int j=0; j<4; ++j)
	A[i][j] += src.A[i][j];
    return(this);
  }

  public final Matrix3D scalarMultiply(double d) {  // this = d * this
    for(int i=0; i<4; ++i) 
      for(int j=0; j<4; ++j)
	A[i][j] = (float)(A[i][j] * d);
    return(this);
  }


  public void loadIdentity() {              // this = identity
    for(int i=0; i<4; ++i)
      for(int j=0; j<4; ++j)
	A[i][j] = (i==j)?1.0f:0.0f;
  }

  public void translate(float tx, float ty, float tz) {  // this = this * t
    float t[] = {1, 0, 0, tx,
                 0, 1, 0, ty,
                 0, 0, 1, tz,
                 0, 0, 0, 1};

    compose(new Matrix3D(t));
  }

  public void scale(float sx, float sy, float sz) {    // this = this * scale
    float s[] = {sx,  0,  0, 0,
		  0, sy,  0, 0,
		  0,  0, sz, 0,
		  0,  0,  0, 1};

    compose(new Matrix3D(s));
  }

  public void skew(float kxy, float kxz, float kyz) {  // this = this * skew
    float s[] = { 1, kxy, kxz, 0,
		  0,   1, kyz, 0,
		  0,   0,   1, 0,
		  0,   0,   0, 1};

    compose(new Matrix3D(s));
  }

  public void rotate(float ax, float ay, float az, float angle) {  
                                                       // this = this * rotate

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

    float Sym[] = { ax*ax, ax*ay, ax*az, 0,
		    ax*ay, ay*ay, ay*az, 0,
		    ax*az, ay*az, az*az, 0,
		    0,     0,     0,     1};

    float Skew[] = { 0,   -az,  ay, 0,
		     az,  0,   -ax, 0,
		     -ay, ax,  0,   0,
		     0,   0,   0,   1};

    
    Matrix3D M = (new Matrix3D(Sym)).scalarMultiply(1 - Math.cos(angle));
    M.add((new Matrix3D(Skew)).scalarMultiply(Math.sin(angle)));
    M.add((new Matrix3D()).scalarMultiply(Math.cos(angle)));
    M.set(3,3,1);
    compose(M);
  }

  private Point3D unitCrossProduct(Point3D v1, Point3D v2) {
    if (v1 == null) throw new InternalError();
    float m[] = {    0, -v1.z,  v1.y, 0,
	          v1.z,     0, -v1.x, 0,
                 -v1.y,  v1.x,     0, 0, 
                     0,     0,     0, 0};

    Point3D o = new Point3D();

    (new Matrix3D(m)).transform1(v2, o);

    float len = (float) Math.sqrt(o.x*o.x + o.y*o.y + o.z*o.z);
    o.x /= len;
    o.y /= len;
    o.z /= len;
    return(o);
  }

  public void lookAt(float eyex, float eyey, float eyez,
		     float atx,  float aty,  float atz,
		     float upx,  float upy,  float upz) { 
                                                       // this = this * lookat
    float lx = atx - eyex;
    float ly = aty - eyey;
    float lz = atz - eyez;

    float len = (float) Math.sqrt(lx*lx + ly*ly + lz*lz);
    lx /= len;
    ly /= len;
    lz /= len;

    Point3D p1 = new Point3D(lx, ly, lz);
    Point3D p2 = new Point3D(upx, upy, upz);

    Point3D r = unitCrossProduct(p1, p2);

    Point3D u = unitCrossProduct(r, new Point3D(lx, ly, lz));

    float RdotI = r.x * eyex + r.y * eyey + r.z * eyez;
    float UdotI = u.x * eyex + u.y * eyey + u.z * eyez;
    float LdotI =  lx * eyex +  ly * eyey +  lz * eyez;

    float v[] = {r.x, r.y, r.z, -RdotI,
                 u.x, u.y, u.z, -UdotI,
                 -lx, -ly, -lz,  LdotI,
                   0,   0,   0,      1};

    compose(new Matrix3D(v));
  }
  //
  // 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) {
    float p[] =
    { 2*near/(right-left), 0, -(right+left)/(right-left), 0,
      0, 2*near/(bottom-top), -(bottom+top)/(bottom-top), 0,
      0,                   0,  (far+near)/(far-near), -2*(far*near)/(far-near),
      0,                   0,              1,            0};
    
    compose(new Matrix3D(p));
  }
  public void orthographic(float left, float right,    // this = this * ortho
			   float bottom, float top,
			   float near, float far) {
    /*
    System.out.println("left: " + left);
    System.out.println("right: " + right);
    System.out.println("bottom: " + bottom);
    System.out.println("top: " + top);
    System.out.println("near: " + near);
    System.out.println("far: " + far);
    */
    
    float o[] =
    { 2/(right-left), 0,              0,            -(right+left)/(right-left),
      0,              2/(bottom-top), 0,            -(bottom+top)/(bottom-top),
      0,              0,              2/(far-near), -(far+near)/(far-near),
      0,              0,              0,            1};

    compose(new Matrix3D(o));
  }

  public String toString() {
    String s = "[";

    for(int i=0; i<4; ++i) {
      for(int j=0; j<4; ++j)
	s += A[i][j] + ((j==3)?"":", ");
      if (i!=3) s += ",\n ";
    }
    s += "]";
    return(s);
  }
}