public abstract class prematrix
{
public double M[][];

//--------- optimization ---------
public static prematrix temp[] = null;
public static vector   tempv[] = null;
//--------------------------------
    
public prematrix()
    {
	M = new double[4][4];
	M[0][0] = 0;
	M[0][1] = 0;
	M[0][2] = 0;
	M[0][3] = 0;
	M[1][0] = 0;
	M[1][1] = 0;
	M[1][2] = 0;
	M[1][3] = 0;
	M[2][0] = 0;
	M[2][1] = 0;
	M[2][2] = 0;
	M[2][3] = 0;
	M[3][0] = 0;
	M[3][1] = 0;
	M[3][2] = 0;
	M[3][3] = 0;
    }

public static void identity(prematrix m)
    {
	m.M[0][0] = 1;
	m.M[0][1] = 0;
	m.M[0][2] = 0;
	m.M[0][3] = 0;
	m.M[1][0] = 0;
	m.M[1][1] = 1;
	m.M[1][2] = 0;
	m.M[1][3] = 0;
	m.M[2][0] = 0;
	m.M[2][1] = 0;
	m.M[2][2] = 1;
	m.M[2][3] = 0;
	m.M[3][0] = 0;
	m.M[3][1] = 0;
	m.M[3][2] = 0;
	m.M[3][3] = 1;	
    }

public static void translate(vector v, prematrix m)
    {
	m.M[0][0] = 1;
	m.M[0][1] = 0;
	m.M[0][2] = 0;
	m.M[0][3] = v.xyz[0];
	m.M[1][0] = 0;
	m.M[1][1] = 1;
	m.M[1][2] = 0;
	m.M[1][3] = v.xyz[1];
	m.M[2][0] = 0;
	m.M[2][1] = 0;
	m.M[2][2] = 1;
	m.M[2][3] = v.xyz[2];
	m.M[3][0] = 0;
	m.M[3][1] = 0;
	m.M[3][2] = 0;
	m.M[3][3] = 1;	
    }

public static void translate(double x, double y, double z, prematrix m)
    {
	m.M[0][0] = 1;
	m.M[0][1] = 0;
	m.M[0][2] = 0;
	m.M[0][3] = x;
	m.M[1][0] = 0;
	m.M[1][1] = 1;
	m.M[1][2] = 0;
	m.M[1][3] = y;
	m.M[2][0] = 0;
	m.M[2][1] = 0;
	m.M[2][2] = 1;
	m.M[2][3] = z;
	m.M[3][0] = 0;
	m.M[3][1] = 0;
	m.M[3][2] = 0;
	m.M[3][3] = 1;	
    }

public static void scale(double x, double y, double z, prematrix m)
    {
	m.M[0][0] = x;
	m.M[0][1] = 0;
	m.M[0][2] = 0;
	m.M[0][3] = 0;
	m.M[1][0] = 0;
	m.M[1][1] = y;
	m.M[1][2] = 0;
	m.M[1][3] = 0;
	m.M[2][0] = 0;
	m.M[2][1] = 0;
	m.M[2][2] = z;
	m.M[2][3] = 0;
	m.M[3][0] = 0;
	m.M[3][1] = 0;
	m.M[3][2] = 0;
	m.M[3][3] = 1;	
    }

public static void rotate(vector r, prematrix m)
    {
	double cx = (double)Math.cos(r.xyz[0]);
	double cy = (double)Math.cos(r.xyz[1]);
	double cz = (double)Math.cos(r.xyz[2]);
	double sx = (double)Math.sin(r.xyz[0]);
	double sy = (double)Math.sin(r.xyz[1]);
	double sz = (double)Math.sin(r.xyz[2]);
	m.M[0][0] = cy*cz;
	m.M[0][1] = -cy*sz;
	m.M[0][2] = sy;
	m.M[0][3] = 0;
	m.M[1][0] = cx*sz + sx*sy*cz;
	m.M[1][1] = cx*cz - sx*sy*sz;
	m.M[1][2] = -sx*cy;
	m.M[1][3] = 0;
	m.M[2][0] = sx*sz - cx*sy*cz;
	m.M[2][1] = sx*cz + cx*sy*sz;
	m.M[2][2] = cx*cy;
	m.M[2][3] = 0;
	m.M[3][0] = 0;
	m.M[3][1] = 0;
	m.M[3][2] = 0;
	m.M[3][3] = 1;
    }

// don't worry, this function is machine generated ;-)
public static void mult(prematrix a, prematrix b, prematrix r)
    {
	r.M[0][0] = a.M[0][0] * b.M[0][0] + a.M[0][1] * b.M[1][0] + a.M[0][2] * b.M[2][0] + a.M[0][3] * b.M[3][0];
	r.M[0][1] = a.M[0][0] * b.M[0][1] + a.M[0][1] * b.M[1][1] + a.M[0][2] * b.M[2][1] + a.M[0][3] * b.M[3][1];
	r.M[0][2] = a.M[0][0] * b.M[0][2] + a.M[0][1] * b.M[1][2] + a.M[0][2] * b.M[2][2] + a.M[0][3] * b.M[3][2];
	r.M[0][3] = a.M[0][0] * b.M[0][3] + a.M[0][1] * b.M[1][3] + a.M[0][2] * b.M[2][3] + a.M[0][3] * b.M[3][3];
	r.M[1][0] = a.M[1][0] * b.M[0][0] + a.M[1][1] * b.M[1][0] + a.M[1][2] * b.M[2][0] + a.M[1][3] * b.M[3][0];
	r.M[1][1] = a.M[1][0] * b.M[0][1] + a.M[1][1] * b.M[1][1] + a.M[1][2] * b.M[2][1] + a.M[1][3] * b.M[3][1];
	r.M[1][2] = a.M[1][0] * b.M[0][2] + a.M[1][1] * b.M[1][2] + a.M[1][2] * b.M[2][2] + a.M[1][3] * b.M[3][2];
	r.M[1][3] = a.M[1][0] * b.M[0][3] + a.M[1][1] * b.M[1][3] + a.M[1][2] * b.M[2][3] + a.M[1][3] * b.M[3][3];
	r.M[2][0] = a.M[2][0] * b.M[0][0] + a.M[2][1] * b.M[1][0] + a.M[2][2] * b.M[2][0] + a.M[2][3] * b.M[3][0];
	r.M[2][1] = a.M[2][0] * b.M[0][1] + a.M[2][1] * b.M[1][1] + a.M[2][2] * b.M[2][1] + a.M[2][3] * b.M[3][1];
	r.M[2][2] = a.M[2][0] * b.M[0][2] + a.M[2][1] * b.M[1][2] + a.M[2][2] * b.M[2][2] + a.M[2][3] * b.M[3][2];
	r.M[2][3] = a.M[2][0] * b.M[0][3] + a.M[2][1] * b.M[1][3] + a.M[2][2] * b.M[2][3] + a.M[2][3] * b.M[3][3];
	r.M[3][0] = a.M[3][0] * b.M[0][0] + a.M[3][1] * b.M[1][0] + a.M[3][2] * b.M[2][0] + a.M[3][3] * b.M[3][0];
	r.M[3][1] = a.M[3][0] * b.M[0][1] + a.M[3][1] * b.M[1][1] + a.M[3][2] * b.M[2][1] + a.M[3][3] * b.M[3][1];
	r.M[3][2] = a.M[3][0] * b.M[0][2] + a.M[3][1] * b.M[1][2] + a.M[3][2] * b.M[2][2] + a.M[3][3] * b.M[3][2];
	r.M[3][3] = a.M[3][0] * b.M[0][3] + a.M[3][1] * b.M[1][3] + a.M[3][2] * b.M[2][3] + a.M[3][3] * b.M[3][3];
    }
    
public static void screen(double nx, double ny, double xmin, double ymin, double xmax, double ymax, prematrix m)
    {
	translate((nx-1)/2, (ny-1)/2, 0, temp[0]);
	scale( -nx/(xmax-xmin), -ny/(ymax-ymin), 1, temp[1]);
	translate((-xmax-xmin)/2, (-ymax-ymin)/2, 0, temp[2]);
	mult(temp[1], temp[2], temp[3]);
	mult(temp[0], temp[3], m);
    }
    
public static void persp(double n, double f, prematrix m)
    {
	m.M[0][0] = 1;
	m.M[0][1] = 0;
	m.M[0][2] = 0;
	m.M[0][3] = 0;
	m.M[1][0] = 0;
	m.M[1][1] = 1;
	m.M[1][2] = 0;
	m.M[1][3] = 0;
	m.M[2][0] = 0;
	m.M[2][1] = 0;
	m.M[2][2] = (n+f)/n;
	m.M[2][3] = -f;
	m.M[3][0] = 0;
	m.M[3][1] = 0;
	m.M[3][2] = 1/n;	
	m.M[3][3] = 0;	
    }

public static void view(vector lf, vector la, vector Vup, prematrix m)
    {
	vector.subu(la, lf, tempv[0]);             // tempv[0] = w
	vector.crossu(Vup, tempv[0], tempv[1]);    // tempv[1] = u
	vector.cross(tempv[0], tempv[1], tempv[2]);// tempv[2] = v

	temp[4].M[0][0] = tempv[1].xyz[0];
	temp[4].M[0][1] = tempv[1].xyz[1];
	temp[4].M[0][2] = tempv[1].xyz[2];
	temp[4].M[0][3] = 0;

	temp[4].M[1][0] = tempv[2].xyz[0];
	temp[4].M[1][1] = tempv[2].xyz[1];
	temp[4].M[1][2] = tempv[2].xyz[2];
	temp[4].M[1][3] = 0;

	temp[4].M[2][0] = tempv[0].xyz[0];
	temp[4].M[2][1] = tempv[0].xyz[1];
	temp[4].M[2][2] = tempv[0].xyz[2];
	temp[4].M[2][3] = 0;

	temp[4].M[3][0] = 0;
	temp[4].M[3][1] = 0;
	temp[4].M[3][2] = 0;
	temp[4].M[3][3] = 1;

	translate(-lf.xyz[0], -lf.xyz[1], -lf.xyz[2], temp[5]);
	mult(temp[4], temp[5], m);
    }
/*
public static void make_mom(vector lf,  vector la,   vector Vup, double n,    double f,
			    double xmax, double  xmin, double ymax, double ymin, double nx,
			    double ny,   double  sx,   double sy,   double sz,   prematrix m)
    {
	scale(sx, sy, sz, temp[6]);                     // mscal
	view(lf, la, Vup, temp[7]);                     // mview
	persp(n, f, temp[8]);                           // mpers
	screen(nx, ny, xmin, ymin, xmax, ymax, temp[9]);// mscrn
	mult(temp[7], temp[6], temp[0]);                // temp[0] = mview * mscale
	mult(temp[8], temp[0], temp[1]);                // temp[1] = mpers * mview * mscale
	mult(temp[9], temp[1], m);                      // m       = mscrn * mpers * mview * mscale
    }

public static void make_mom(prematrix mscal, prematrix mview, prematrix mpers, prematrix mscrn, prematrix m)
{
	mult(mview, mscal, temp[6]);   // temp[6] = mview * mscale
	mult(mpers, temp[6], temp[7]); // temp[7] = mpers * mview * mscale
	mult(mscrn, temp[7], m);       // m       = mscrn * mpers * mview * mscale
    }
    */
public static void make_mom(prematrix mview, prematrix mpers, prematrix mscrn, prematrix m)
    {
	mult(mpers, mview, temp[6]); // temp[6] = mpers * mview
	mult(mscrn, temp[6], m);     // m       = mscrn * mpers * mview
    }
    
// this function auto homogenizes the answer so it fits into 3x1 vector
public static void mult(prematrix a, vector b, vector r)
    {
	double x = a.M[0][0]*b.xyz[0] + a.M[0][1]*b.xyz[1] + a.M[0][2]*b.xyz[2] +  a.M[0][3];
	double y = a.M[1][0]*b.xyz[0] + a.M[1][1]*b.xyz[1] + a.M[1][2]*b.xyz[2] +  a.M[1][3];
	double z = a.M[2][0]*b.xyz[0] + a.M[2][1]*b.xyz[1] + a.M[2][2]*b.xyz[2] +  a.M[2][3];
	double h = a.M[3][0]*b.xyz[0] + a.M[3][1]*b.xyz[1] + a.M[3][2]*b.xyz[2] +  a.M[3][3];
	r.xyz[0] = x/h;
	r.xyz[1] = y/h;
	r.xyz[2] = z/h;
    }

public static void assign(prematrix a, prematrix b) // a <- b
    {
	a.M[0][0] = b.M[0][0];
	a.M[0][1] = b.M[0][1];
	a.M[0][2] = b.M[0][2];
	a.M[0][3] = b.M[0][3];
	a.M[1][0] = b.M[1][0];
	a.M[1][1] = b.M[1][1];
	a.M[1][2] = b.M[1][2];
	a.M[1][3] = b.M[1][3];
	a.M[2][0] = b.M[2][0];
	a.M[2][1] = b.M[2][1];
	a.M[2][2] = b.M[2][2];
	a.M[2][3] = b.M[2][3];
	a.M[3][0] = b.M[3][0];
	a.M[3][1] = b.M[3][1];
	a.M[3][2] = b.M[3][2];	
	a.M[3][3] = b.M[3][3];
    }
    
public void dump(String mesg)
    {
	System.out.println(mesg);
	System.out.println("|"+M[0][0]+"\t"+M[0][1]+"\t"+M[0][2]+"\t"+M[0][3]+"\t|");
	System.out.println("|"+M[1][0]+"\t"+M[1][1]+"\t"+M[1][2]+"\t"+M[1][3]+"\t|");
	System.out.println("|"+M[2][0]+"\t"+M[2][1]+"\t"+M[2][2]+"\t"+M[2][3]+"\t|");
	System.out.println("|"+M[3][0]+"\t"+M[3][1]+"\t"+M[3][2]+"\t"+M[3][3]+"\t|");
	System.out.println();
    }

public abstract void recalc(double t);
}