/*
 * @(#)Matrix4D.java
 *
 * Name: Matrix4D - A 4-dimensional Matrix Object
 * Version: Java Applet V1.2, last revised: 2007-02-01
 * Author: Harry J. Smith, 19628 Via Monte Dr., Saratoga CA 95070
 * Copyright (c) 1999-2007 by author, All rights reserved
 *
 * Based on Matrix3D.java by Sun Microsystems, Inc.
 */

/* A fairly conventional 4D matrix object that can transform sets of
   4D points and perform a variety of manipulations on the transform */
class Matrix4D {
	float xx, xy, xz, xh, xo;
	float yx, yy, yz, yh, yo;
	float zx, zy, zz, zh, zo;
	float hx, hy, hz, hh, ho;
	static final double pi = 4 * Math.atan(1); // 3.1415926535898

	// Create a new unit matrix
	Matrix4D () {
		xx = 1.0f;  yy = 1.0f;  zz = 1.0f;  hh = 1.0f;
	}

	// Scale by f in all dimensions
	void scale(float f) {
		xx *= f;  xy *= f;  xz *= f;  xh *= f;  xo *= f;
		yx *= f;  yy *= f;  yz *= f;  yh *= f;  yo *= f;
		zx *= f;  zy *= f;  zz *= f;  zh *= f;  zo *= f;
		hx *= f;  hy *= f;  hz *= f;  hh *= f;  ho *= f;
	}

	// Scale along each axis independently
	void scale(float xf, float yf, float zf, float hf) {
		xx *= xf;  xy *= xf;  xz *= xf;  xh *= xf;  xo *= xf;
		yx *= yf;  yy *= yf;  yz *= yf;  yh *= yf;  yo *= yf;
		zx *= zf;  zy *= zf;  zz *= zf;  zh *= zf;  zo *= zf;
		hx *= hf;  hy *= hf;  hz *= hf;  hh *= hf;  ho *= hf;
	}

	// Translate the origin
	void translate(float x, float y, float z, float h) {
		xo += x;  yo += y;  zo += z;  ho += h;
	}

	// rotate theta degrees in the xz plane
	void xzrot(double theta) {
		theta *= (pi / 180);
		double ct = Math.cos(theta);
		double st = Math.sin(theta);

		float Nxx = (float) (xx * ct + zx * st);
		float Nxy = (float) (xy * ct + zy * st);
		float Nxz = (float) (xz * ct + zz * st);
		float Nxh = (float) (xh * ct + zh * st);
		float Nxo = (float) (xo * ct + zo * st);

		float Nzx = (float) (zx * ct - xx * st);
		float Nzy = (float) (zy * ct - xy * st);
		float Nzz = (float) (zz * ct - xz * st);
		float Nzh = (float) (zh * ct - xh * st);
		float Nzo = (float) (zo * ct - xo * st);

		xx = Nxx;  xy = Nxy;  xz = Nxz;  xh = Nxh;  xo = Nxo;
		zx = Nzx;  zy = Nzy;  zz = Nzz;  zh = Nzh;  zo = Nzo;
	}

	// rotate theta degrees in the yz plane
	void yzrot(double theta) {
		theta *= (pi / 180);
		double ct = Math.cos(theta);
		double st = Math.sin(theta);

		float Nyx = (float) (yx * ct + zx * st);
		float Nyy = (float) (yy * ct + zy * st);
		float Nyz = (float) (yz * ct + zz * st);
		float Nyh = (float) (yh * ct + zh * st);
		float Nyo = (float) (yo * ct + zo * st);

		float Nzx = (float) (zx * ct - yx * st);
		float Nzy = (float) (zy * ct - yy * st);
		float Nzz = (float) (zz * ct - yz * st);
		float Nzh = (float) (zh * ct - yh * st);
		float Nzo = (float) (zo * ct - yo * st);

		yx = Nyx;  yy = Nyy;  yz = Nyz;  yh = Nyh;  yo = Nyo;
		zx = Nzx;  zy = Nzy;  zz = Nzz;  zh = Nzh;  zo = Nzo;
	}

	// rotate theta degrees in the xy plane
	void xyrot(double theta) {
		theta *= (pi / 180);
		double ct = Math.cos(theta);
		double st = Math.sin(theta);

		float Nyx = (float) (yx * ct + xx * st);
		float Nyy = (float) (yy * ct + xy * st);
		float Nyz = (float) (yz * ct + xz * st);
		float Nyh = (float) (yh * ct + xh * st);
		float Nyo = (float) (yo * ct + xo * st);

		float Nxx = (float) (xx * ct - yx * st);
		float Nxy = (float) (xy * ct - yy * st);
		float Nxz = (float) (xz * ct - yz * st);
		float Nxh = (float) (xh * ct - yh * st);
		float Nxo = (float) (xo * ct - yo * st);

		yx = Nyx;  yy = Nyy;  yz = Nyz;  yh = Nyh;  yo = Nyo;
		xx = Nxx;  xy = Nxy;  xz = Nxz;  xh = Nxh;  xo = Nxo;
	}

	// rotate theta degrees in the xh plane
	void xhrot(double theta) {
		theta *= (pi / 180);
		double ct = Math.cos(theta);
		double st = Math.sin(theta);

		float Nhx = (float) (hx * ct + xx * st);
		float Nhy = (float) (hy * ct + xy * st);
		float Nhz = (float) (hz * ct + xz * st);
		float Nhh = (float) (hh * ct + xh * st);
		float Nho = (float) (ho * ct + xo * st);

		float Nxx = (float) (xx * ct - hx * st);
		float Nxy = (float) (xy * ct - hy * st);
		float Nxz = (float) (xz * ct - hz * st);
		float Nxh = (float) (xh * ct - hh * st);
		float Nxo = (float) (xo * ct - ho * st);

		hx = Nhx;  hy = Nhy;  hz = Nhz;  hh = Nhh;  ho = Nho;
		xx = Nxx;  xy = Nxy;  xz = Nxz;  xh = Nxh;  xo = Nxo;
	}

	// rotate theta degrees in the yh plane
	void yhrot(double theta) {
		theta *= (pi / 180);
		double ct = Math.cos(theta);
		double st = Math.sin(theta);

		float Nhx = (float) (hx * ct + yx * st);
		float Nhy = (float) (hy * ct + yy * st);
		float Nhz = (float) (hz * ct + yz * st);
		float Nhh = (float) (hh * ct + yh * st);
		float Nho = (float) (ho * ct + yo * st);

		float Nyx = (float) (yx * ct - hx * st);
		float Nyy = (float) (yy * ct - hy * st);
		float Nyz = (float) (yz * ct - hz * st);
		float Nyh = (float) (yh * ct - hh * st);
		float Nyo = (float) (yo * ct - ho * st);

		hx = Nhx;  hy = Nhy;  hz = Nhz;  hh = Nhh;  ho = Nho;
		yx = Nyx;  yy = Nyy;  yz = Nyz;  yh = Nyh;  yo = Nyo;
	}

	// Multiply this matrix by a second: M = M*R
	void mult(Matrix4D rhs) {
		float lxx = xx * rhs.xx + yx * rhs.xy + zx * rhs.xz + hx * rhs.xh;
		float lxy = xy * rhs.xx + yy * rhs.xy + zy * rhs.xz + hy * rhs.xh;
		float lxz = xz * rhs.xx + yz * rhs.xy + zz * rhs.xz + hz * rhs.xh;
		float lxh = xh * rhs.xx + yh * rhs.xy + zh * rhs.xz + hh * rhs.xh;
		float lxo = xo * rhs.xx + yo * rhs.xy + zo * rhs.xz + ho * rhs.xh + rhs.xo;

		float lyx = xx * rhs.yx + yx * rhs.yy + zx * rhs.yz + hx * rhs.yh;
		float lyy = xy * rhs.yx + yy * rhs.yy + zy * rhs.yz + hy * rhs.yh;
		float lyz = xz * rhs.yx + yz * rhs.yy + zz * rhs.yz + hz * rhs.yh;
		float lyh = xh * rhs.yx + yh * rhs.yy + zh * rhs.yz + hh * rhs.yh;
		float lyo = xo * rhs.yx + yo * rhs.yy + zo * rhs.yz + ho * rhs.yh + rhs.yo;

		float lzx = xx * rhs.zx + yx * rhs.zy + zx * rhs.zz + hx * rhs.zh;
		float lzy = xy * rhs.zx + yy * rhs.zy + zy * rhs.zz + hy * rhs.zh;
		float lzz = xz * rhs.zx + yz * rhs.zy + zz * rhs.zz + hz * rhs.zh;
		float lzh = xh * rhs.zx + yh * rhs.zy + zh * rhs.zz + hh * rhs.zh;
		float lzo = xo * rhs.zx + yo * rhs.zy + zo * rhs.zz + ho * rhs.zh + rhs.zo;

		float lhx = xx * rhs.hx + yx * rhs.hy + zx * rhs.hz + hx * rhs.hh;
		float lhy = xy * rhs.hx + yy * rhs.hy + zy * rhs.hz + hy * rhs.hh;
		float lhz = xz * rhs.hx + yz * rhs.hy + zz * rhs.hz + hz * rhs.hh;
		float lhh = xh * rhs.hx + yh * rhs.hy + zh * rhs.hz + hh * rhs.hh;
		float lho = xo * rhs.hx + yo * rhs.hy + zo * rhs.hz + ho * rhs.hh + rhs.ho;

		xx = lxx;  xy = lxy;  xz = lxz;  xh = lxh;  xo = lxo;
		yx = lyx;  yy = lyy;  yz = lyz;  yh = lyh;  yo = lyo;
		zx = lzx;  zy = lzy;  zz = lzz;  zh = lzh;  zo = lzo;
		hx = lhx;  hy = lhy;  hz = lhz;  hh = lhh;  ho = lho;
	}

	// Reinitialize to the unit matrix
	void unit() {
		xx = 1;  xy = 0;  xz = 0;  xh = 0;  xo = 0;
		yx = 0;  yy = 1;  yz = 0;  yh = 0;  yo = 0;
		zx = 0;  zy = 0;  zz = 1;  zh = 0;  zo = 0;
		hx = 0;  hy = 0;  hz = 0;  hh = 1;  ho = 0;
	}

	/* Transform nvert points from v into tv. v contains the input
	   coordinates in floating point. Three successive entries in
	   the array constitute a point. tv ends up holding the transformed
	   points as integers; three successive entries per point */
	void transform(float v[], int tv[], int nvert) {
		float lxx = xx, lxy = xy, lxz = xz, lxh = xh, lxo = xo;
		float lyx = yx, lyy = yy, lyz = yz, lyh = yh, lyo = yo;
		float lzx = zx, lzy = zy, lzz = zz, lzh = zh, lzo = zo;
		float lhx = hx, lhy = hy, lhz = hz, lhh = hh, lho = ho;
		for (int i = nvert * 4; (i -= 4) >= 0;) {
			float x = v[i];
			float y = v[i + 1];
			float z = v[i + 2];
			float h = v[i + 3];
			tv[i	] = (int) (x * lxx + y * lxy + z * lxz + h * lxh + lxo);
			tv[i + 1] = (int) (x * lyx + y * lyy + z * lyz + h * lyh + lyo);
			tv[i + 2] = (int) (x * lzx + y * lzy + z * lzz + h * lzh + lzo);
			tv[i + 3] = (int) (x * lhx + y * lhy + z * lhz + h * lhh + lho);
		}
	}

	public String toString() {
		return ("[" + xx + "," + xy + "," + xz + "," + xh + "," + xo + ";"
					+ yx + "," + yy + "," + yz + "," + yh + "," + yo + ";"
					+ zx + "," + zy + "," + zz + "," + zh + "," + zo + ";"
					+ hx + "," + hy + "," + hz + "," + hh + "," + ho + "]");
	}
} // end of class Matrix4D