/*
 * Copyright (c) 1996 CONNIE K. PENG  All Rights Reserved.
 *
 * Permission to use, copy, modify, and distribute this software
 * and its documentation for NON-COMMERCIAL purposes.
 *
 */
public class Matrix 
{
  float[][] element;

  public Matrix() 
  {
    element = new float[4][4];

    this.reset();
  }

  public Matrix(Matrix m) 
  {
    element = new float[4][4];

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

  public void reset ()
  {
    for(int i = 0; i < 4; ++i)
      for(int j = 0; j < 4; ++j)
        if (i == j)
          element[i][j] = 1;
        else
          element[i][j] = 0;
  }

  public void multiply(Matrix m) 
  {
    float[][] e = new float[4][4];
        
    for(int i = 0; i < 4; ++i)
      for(int j = 0; j < 4; ++j)
        for(int k = 0; k < 4; ++k)
          e[i][j] += element[k][j] * m.element[i][k];

    for(int i = 0; i < 4; ++i)
      for(int j = 0; j < 4; ++j)
          element[i][j] = e[i][j];
  }

  public Matrix inverse()
  {
    Matrix m = new Matrix();

    // from row 1 to row 4

    for (int j = 0; j < 4; j++)
    {
      // normalize the diagonal element

      float e = element[j][j];
      for (int i = 0; i < 4; i++)
      {
        element[i][j] = element[i][j] / e;
        m.element[i][j] = m.element[i][j] /e;
      }
      
      // make other element to be 0

      for (int k = 0; k < 4; k++)
      {
        if (k != j)
        {
          e = -element[j][k];
          for (int i = 0; i < 4; i++)
          {
            element[i][k] = element[i][k] + e * element[i][j];
            m.element[i][k] = m.element[i][k] + e * m.element[i][j];
          }
        }
      }
    }
    return m;
  }

  public void scale(float f) 
  {
    Matrix m = new Matrix();

    for(int i = 0; i < 4; ++i)
        element[i][i] = f;

    this.multiply(m);
  }

  public void translate(Vector vec)
  {
    translate(vec.element[0], vec.element[1], vec.element[2]);
  }

  public void translate (float x, float y, float z)
  {
    Matrix m = new Matrix();
    m.element[0][3] = x;
    m.element[1][3] = y;
    m.element[2][3] = z;

    this.multiply(m);
  }

  // rotate along x axis
  public void rotatex(int degree)
  {
    float d;
    float cos, sin;

    d = (float) (degree * Math.PI/180.);
    sin = (float) Math.sin(d);
    cos = (float) Math.cos(d);

    rotatex(sin, cos);
  }

  // rotate along x axis
  public void rotatex(float sin, float cos)
  {
    Matrix m = new Matrix();

    m.element[1][1] =  cos;
    m.element[2][2] =  cos;
    m.element[1][2] = -sin;
    m.element[2][1] =  sin;

    this.multiply(m);
  }

  // rotate along y axis
  public void rotatey(int degree)
  {
    float d;
    float cos, sin;

    d = (float) (degree * Math.PI/180.);
    sin = (float) Math.sin(d);
    cos = (float) Math.cos(d);

    rotatey(sin, cos);
  }

  // rotate along y axis
  public void rotatey(float sin, float cos)
  {
    Matrix m = new Matrix();

    m.element[0][0] =  cos;
    m.element[2][2] =  cos;
    m.element[0][2] =  sin;
    m.element[2][0] = -sin;

    this.multiply(m);
  }

  // rotate along z axis
  public void rotatez(int degree)
  {
    float d;
    float cos, sin;

    d = (float) (degree * Math.PI/180.);
    sin = (float) Math.sin(d);
    cos = (float) Math.cos(d);

    rotatez(sin, cos);
  }

  // rotate along z axis
  public void rotatez(float sin, float cos)
  {
    Matrix m = new Matrix();

    m.element[0][0] =  cos;
    m.element[1][1] =  cos;
    m.element[0][1] = -sin;
    m.element[1][0] =  sin;

    this.multiply(m);
  }

  // rotate along a given axis and point 
  public void rotate(Vector point, Vector axis, int angle)
  {
    float sinx, cosx;   // angle rotate along x axis
    float siny, cosy;	// angle rotate along y axis
    float yz_m;		// the magnitude of the yz component of axis
    float m;		// total magnitude of the axis vector

    Matrix matrix = new Matrix();

    yz_m = (float)Math.sqrt(axis.element[1] * axis.element[1] + 
                            axis.element[2] * axis.element[2]);
    if (yz_m == 0)
    {
      sinx = 0;
      cosx = 1;
    }
    else
    {
      sinx = axis.element[1]/yz_m;
      cosx = axis.element[2]/yz_m;
    }

    m = axis.magnitude();
    if (m == 0)
    {
      siny = 0;
      cosy = 1;
    }
    else
    {
      siny = -axis.element[0]/m;
      cosy =  yz_m/m;
    }

    // first map the origin to the point, then rotate z axis around x axis
    // (up/down), then rotate z axix around y axis (left/right), then
    // this map the z axix to be the axix vector given, then rotate around
    // this z axix to the angle given.

    Vector origin = new Vector(point);
    matrix.translate(origin.negative());
    matrix.rotatex(sinx, cosx);
    matrix.rotatey(siny, cosy);
    matrix.rotatez(angle);

    // reverse the above process

    matrix.rotatey(-siny, cosy);
    matrix.rotatex(-sinx, cosx);
    matrix.translate(point);

    this.multiply(matrix);
  }

  public void print()
  {
    for (int j = 0; j < 4; j++)
      System.out.println(element[0][j] + " " + element[1][j] +
                         " " + element[2][j] + " " + element[3][j]);
  }
}