/** 4x4 Transform Matrix - designed to do translation, scale, rotation and perspective.
    This matrix assumes right handed coordinate systems, and that it will
    be multiplied on the right by a column vector (or matrix) like : Ax = B */
    
class XFmatrix implements Cloneable {
    float[] m;
    static final double pi = Math.PI;

    /** Create a new identity matrix */
    XFmatrix () {
	    m = new float[16];
		m[0] = 1.0f;
		m[5] = 1.0f;
		m[10] = 1.0f;
		m[15] = 1.0f;
	}
	
	/** create a copy of the matrix */
	public Object clone() // overrides Object.clone()
	{
		XFmatrix r = new XFmatrix();
		System.arraycopy(m,0,r.m,0,m.length);
		return r;
	}
	
	
    /**  Scale by f in all dimensions. 
   	  * Assumes no scaling in w-coord, this is like mult on left by scale matrix */
    void scale(float f) {
		for(int i =0; i<16; i++) {
			m[i] *= f;
		}	
    }	
    
    /** Combine a scale matrix with this. Scale along each axis independently 
    	mult on left by scale matrix */
    void scale(float xf, float yf, float zf) {
	m[0] *= xf;
	m[1] *= xf;
	m[2] *= xf;
	m[3] *= xf;
	
	m[4] *= yf;
	m[5] *= yf;
	m[6] *= yf;
	m[7] *= yf;
	
	m[8] *= zf;
	m[9] *= zf;
	m[10] *= zf;
	m[11] *= zf;
	
    }
    
    /** Combine a translation matrix
      *	Assumes mult on left by translation matrix */
    void translate(float xf, float yf, float zf) {
	m[0] += xf*m[12];
	m[1] += xf*m[13];
	m[2] += xf*m[14];
	m[3] += xf;
	
	m[4] += yf*m[12];
	m[5] += yf*m[13];
	m[6] += yf*m[14];
	m[7] += yf;

	m[8] += zf*m[12];
	m[9] += zf*m[13];
	m[10] += zf*m[14];
	m[11] += zf;
    }
     
    /** Combine a rotation matrix.
      * Rotate theta degrees about the y axis */
	void yrot(double theta) {
	theta *= (pi / 180);
	double ct = Math.cos(theta);
	double st = Math.sin(theta);

	float Nxx = (float) (m[0] * ct + m[8] * st);
	float Nxy = (float) (m[1] * ct + m[9] * st);
	float Nxz = (float) (m[2] * ct + m[10] * st);
	float Nxo = (float) (m[3] * ct + m[11] * st);

	float Nzx = (float) (m[8] * ct - m[0] * st);
	float Nzy = (float) (m[9] * ct - m[1] * st);
	float Nzz = (float) (m[10] * ct - m[2] * st);
	float Nzo = (float) (m[11] * ct - m[3] * st);

	m[3] = Nxo;
	m[0] = Nxx;
	m[1] = Nxy;
	m[2] = Nxz;
	
	m[11] = Nzo;
	m[8] = Nzx;
	m[9] = Nzy;
	m[10] = Nzz;
    }

    /** Combine a rotation matrix.
      * Rotate theta degrees about the x axis */
    void xrot(double theta) {
	theta *= (pi / 180);
	double ct = Math.cos(theta);
	double st = Math.sin(theta);

	float Nyx = (float) (m[4] * ct + m[8] * st);
	float Nyy = (float) (m[5] * ct + m[9] * st);
	float Nyz = (float) (m[6] * ct + m[10] * st);
	float Nyo = (float) (m[7] * ct + m[11] * st);

	float Nzx = (float) ( m[8] * ct - m[4] * st );
	float Nzy = (float) ( m[9] * ct - m[5] * st );
	float Nzz = (float) ( m[10] * ct - m[6] * st );
	float Nzo = (float) ( m[11] * ct - m[7] * st );

	m[4] = Nyx;
	m[5] = Nyy;
	m[6] = Nyz;
	m[7] = Nyo;
	
	m[8] = Nzx;
	m[9] = Nzy;
	m[10] = Nzz;
	m[11] = Nzo;
    }

    /** Combine a rotation matrix.
      * Rotate theta degrees about the z axis */
    void zrot(double theta) {
	theta *= (pi / 180);
	double ct = Math.cos(theta);
	double st = Math.sin(theta);

	float Nyx = (float) ( m[4] * ct + m[0] * st);
	float Nyy = (float) ( m[5] * ct + m[1] * st);
	float Nyz = (float) ( m[6] * ct + m[2] * st);
	float Nyo = (float) ( m[7] * ct + m[3] * st);

	float Nxx = (float) ( m[0] * ct - m[4] * st);
	float Nxy = (float) ( m[1] * ct - m[5] * st);
	float Nxz = (float) ( m[2] * ct - m[6] * st);
	float Nxo = (float) ( m[3] * ct - m[7] * st);

	m[4] = Nyx;
	m[5] = Nyy;
	m[6] = Nyz;
	m[7] = Nyo;
	
	m[0] = Nxx;
	m[1] = Nxy;
	m[2] = Nxz;
 	m[3] = Nxo;
   }

    /** Combine a rotation matrix.
      * Rotate to a new coord system, with x,y,z being the new axes */
    void rot( float x0, float x1, float x2, 
    		  float y0, float y1, float y2, 
    	 	  float z0, float z1, float z2)
    {
	float Nxx = (float) ( m[0] * x0 + m[4] * x1 + m[8] * x2 );
	float Nxy = (float) ( m[1] * x0 + m[5] * x1 + m[9] * x2 );
	float Nxz = (float) ( m[2] * x0 + m[6] * x1 + m[10] * x2 );
	float Nxo = (float) ( m[3] * x0 + m[7] * x1 + m[11] * x2 );
		
	float Nyx = (float) ( m[0] * y0 + m[4] * y1 + m[8] * y2 );
	float Nyy = (float) ( m[1] * y0 + m[5] * y1 + m[9] * y2 );
	float Nyz = (float) ( m[2] * y0 + m[6] * y1 + m[10] * y2 );
	float Nyo = (float) ( m[3] * y0 + m[7] * y1 + m[11] * y2 );

	float Nzx = (float) ( m[0] * z0 + m[4] * z1 + m[8] * z2 );
	float Nzy = (float) ( m[1] * z0 + m[5] * z1 + m[9] * z2 );
	float Nzz = (float) ( m[2] * z0 + m[6] * z1 + m[10] * z2 );
	float Nzo = (float) ( m[3] * z0 + m[7] * z1 + m[11] * z2 );
	
	m[0] = Nxx;
	m[1] = Nxy;
	m[2] = Nxz;
	m[3] = Nxo;

	m[4] = Nyx;
	m[5] = Nyy;
	m[6] = Nyz;
	m[7] = Nyo;

	m[8] = Nzx;
	m[9] = Nzy;
	m[10] = Nzz;
	m[11] = Nzo;
	}    	

    /** Combine a rotation matrix.
      * Rotate to a new coord system, with Gvectors Vx,Vy,Vz being the new axes */
    void rot( Gvector Vx, Gvector Vy, Gvector Vz)
    {
    	this.rot(Vx.x, Vx.y, Vx.z, 
    			 Vy.x, Vy.y, Vy.z, 
    			 Vz.x, Vz.y, Vz.z);
    } 
    
    /** Combine a perspective matrix.
      * Computes the matrix from the hither and yon planes, assuming
      * hither and yon are specified as distances along the DOV vector */
    void persp ( float hither, float yon)
    {
    	float zmin = -hither/yon;
    	float p22 = 1f/(1f + zmin);
    	float p23 = -zmin/(1f + zmin);
    	float p32 = -1f;

	float Nzx = (float) ( m[8] * p22 + m[12] * p23 );
	float Nzy = (float) ( m[9] * p22 + m[13] * p23 );
	float Nzz = (float) ( m[10] * p22 + m[14] * p23 );
	float Nzo = (float) ( m[11] * p22 + m[15] * p23 );

	m[12] = -m[8];   // w value gets -z
	m[13] = -m[9];
	m[14] = -m[10];
	m[15] = -m[11];

	m[8] = Nzx;
	m[9] = Nzy;
	m[10] = Nzz;
	m[11] = Nzo;
    	
    }
    
    /** Replace current values with the product of two XFmatrices. */
    void combine( XFmatrix L, XFmatrix R)
    {
		float X0 = L.m[0]*R.m[0] + L.m[1]*R.m[4] + L.m[2]*R.m[8] +  L.m[3]*R.m[12];
		float X1 = L.m[0]*R.m[1] + L.m[1]*R.m[5] + L.m[2]*R.m[9] +  L.m[3]*R.m[13];
		float X2 = L.m[0]*R.m[2] + L.m[1]*R.m[6] + L.m[2]*R.m[10] + L.m[3]*R.m[14];
		float X3 = L.m[0]*R.m[3] + L.m[1]*R.m[7] + L.m[2]*R.m[11] + L.m[3]*R.m[15];
		
		float X4 = L.m[4]*R.m[0] + L.m[5]*R.m[4] + L.m[6]*R.m[8] +  L.m[7]*R.m[12];
		float X5 = L.m[4]*R.m[1] + L.m[5]*R.m[5] + L.m[6]*R.m[9] +  L.m[7]*R.m[13];
		float X6 = L.m[4]*R.m[2] + L.m[5]*R.m[6] + L.m[6]*R.m[10] + L.m[7]*R.m[14];
		float X7 = L.m[4]*R.m[3] + L.m[5]*R.m[7] + L.m[6]*R.m[11] + L.m[7]*R.m[15];
		
		float X8 = L.m[8]*R.m[0] + L.m[9]*R.m[4] + L.m[10]*R.m[8] +  L.m[11]*R.m[12];
		float X9 = L.m[8]*R.m[1] + L.m[9]*R.m[5] + L.m[10]*R.m[9] +  L.m[11]*R.m[13];
		float X10 = L.m[8]*R.m[2] + L.m[9]*R.m[6] + L.m[10]*R.m[10] + L.m[11]*R.m[14];
		float X11 = L.m[8]*R.m[3] + L.m[9]*R.m[7] + L.m[10]*R.m[11] + L.m[11]*R.m[15];

		float X12 = L.m[12]*R.m[0] + L.m[13]*R.m[4] + L.m[14]*R.m[8] +  L.m[15]*R.m[12];
		float X13 = L.m[12]*R.m[1] + L.m[13]*R.m[5] + L.m[14]*R.m[9] +  L.m[15]*R.m[13];
		float X14 = L.m[12]*R.m[2] + L.m[13]*R.m[6] + L.m[14]*R.m[10] + L.m[15]*R.m[14];
		float X15 = L.m[12]*R.m[3] + L.m[13]*R.m[7] + L.m[14]*R.m[11] + L.m[15]*R.m[15];
		
		m[0] = X0;
		m[1] = X1;
		m[2] = X2;
		m[3] = X3;
		m[4] = X4;
		m[5] = X5;
		m[6] = X6;
		m[7] = X7;
		m[8] = X8;
		m[9] = X9;
		m[10] = X10;
		m[11] = X11;
		m[12] = X12;
		m[13] = X13;
		m[14] = X14;
		m[15] = X15;
	}

    /** Multiply this matrix by a second: M = M*rhs 
    	so rhs must be 4xN matrix. */
/*    float[][] transformFloats(float[][] rhs) {
    	
    	float[][] r = new float[4][rhs[0].length];
		for(int i = 0; i< 4; i++) {
			for (int j = 0; j < rhs[0].length; j++) {
					r[i][j] = m[i][0]*rhs[0][j] + m[i][1]*rhs[1][j] +
					          m[i][2]*rhs[2][j] + m[i][3]*rhs[3][j];
			}
		}
//		System.arraycopy(r,0,m,0,r.length());		
		return r;
    }*/

	/** Transform a list of G-points. Treats the g-points as if they are
	  * arranged in columns in a matrix, and mult like this: M x in = out 
	  * Assumes none of the Gpoints in lists are null!!! */
    public void transform(Gpoint in[], Gpoint out[], int vertices)
    {
        float x, y, z, w;
        for (int i = 0; i < vertices; i++) {
            x = m[0]*in[i].x + m[1]*in[i].y + m[2]*in[i].z + m[3];
            y = m[4]*in[i].x + m[5]*in[i].y + m[6]*in[i].z + m[7];
            z = m[8]*in[i].x + m[9]*in[i].y + m[10]*in[i].z + m[11];
            w = m[12]*in[i].x + m[13]*in[i].y + m[14]*in[i].z + m[15];
            
            out[i].x = x;
            out[i].y = y;
            out[i].z = z;
			out[i].w = w;
			out[i].r = in[i].r;
			out[i].g = in[i].g;
			out[i].b = in[i].b;
			out[i].valid = in[i].valid;
        }
    }

    /** Reinitialize to the unit matrix */
    void unit() {
	m[0] = 1; m[1] = 0; m[2] = 0; m[3] = 0;
	m[4] = 0; m[5] = 1; m[6] = 0; m[7] = 0;
	m[8] = 0; m[9] = 0; m[10] = 1; m[11] = 0;
	m[12] = 0; m[13] = 0; m[14] = 0; m[15] = 1;
   }

	/** toString is helpful for debugging */
    public String toString() {
	return ("[" + m[0] + "," + m[1] + "," + m[2] + "," + m[3] + ";\n"
		  + m[4] + "," + m[5] + "," + m[6] + "," + m[7] +";\n"
		+ m[8] + "," + m[9] + "," + m[10] + "," + m[11] + ";\n"
		+ m[12] + "," + m[13] + "," + m[14] + "," + m[15] + "]\n");
    }
}