using System;
using System.Collections.Generic;
using System.Text;
using System.Drawing;
using Spherical.HTM;
using Spherical;
namespace hcalc
{
	/// <summary>
	/// 
	/// </summary>
	class Camera3D
	{

		private int[] viewport;
		private double[] datawindow;
		private double _scalex, _scaley;
		double[,] MRTH;

		private void init_components() {
			MRTH = new double[3, 3];
			M_setrot(new Cartesian(0.0, 0.0, 1.0, false), 0.0);
			viewport = new int[4];
			viewport[0] = 0;	// initial size; allows clipping later
			viewport[1] = 0;
			viewport[2] = 320;
			viewport[3] = 320;
			datawindow = new double[4];
			// make it so that origin is centered 1 is actually
			// and +/-1 touches the top/bottom of viewport
			// 
			datawindow[0] = -1.0;
			datawindow[1] = -1.0;
			datawindow[2] = 1.0;
			datawindow[3] = 1.0;
			_scalex = (viewport[2] - viewport[0]) / (datawindow[2] - datawindow[0]);
			_scaley = (viewport[3] - viewport[1]) / (datawindow[3] - datawindow[1]);
		}
		public void setZoom(double zoom) {
			datawindow[0] = -zoom;
			datawindow[1] = -zoom;
			datawindow[2] = zoom;
			datawindow[3] = zoom;
			_scalex = (viewport[2] - viewport[0]) / (datawindow[2] - datawindow[0]);
			_scaley = (viewport[3] - viewport[1]) / (datawindow[3] - datawindow[1]);
		}
		void XXXparallelProject(double x, double y, double z, out double u, out double v) {
			u = (x - datawindow[0]) * _scalex + viewport[0];
			v = (y - datawindow[1]) * _scaley + viewport[1];
		}
		public Camera3D() {
			init_components();
		}
		public void reset() {
			M_setrot(new Cartesian(0.0, 0.0, 1.0, false), 0.0);
		}
		public void set_size(int w, int h) {

			viewport[2] = viewport[0] + w;
			viewport[3] = viewport[1] + h;
			return;
		}
		public void set_camera(double ra, double dec) {
			double x, y, z;
			Cartesian.Radec2Xyz(ra, dec, out x, out y, out z);
			// make a matrix that rotates ra, dec into (0.0, 90) or (0, 0, 1)
			double[] P = { x, y, z };
			double[] Q = {0.0, 0.0, 1.0};
			M_set_rotinto(P, Q);
			return;
		}
		public void vp2dw(int u, int v, out double x, out double y) {
			v = viewport[3] - v;
			x = (u - viewport[0]) / _scalex + datawindow[0];
			y = (v - viewport[1]) / _scaley + datawindow[1];
		}
		public void dw2vp(Cartesian v, ref Point p) {
			p.X = (int)((v.X - datawindow[0]) * _scalex + viewport[0]);
			p.Y = viewport[3] - ((int)((v.Y - datawindow[1]) * _scaley + viewport[1]));
		}
		public void dw2vp(double x, double y, out int ox, out int oy) {
			ox = (int)((x - datawindow[0]) * _scalex + viewport[0]);
			oy = viewport[3] - ((int)((y - datawindow[1]) * _scaley + viewport[1]));
		}

		public static void rotate(Cartesian pin, double[,] mat, ref Cartesian pout) {
			double x = mat[0, 0] * pin.X + mat[0, 1] * pin.Y + mat[0, 2] * pin.Z;
			double y = mat[1, 0] * pin.X + mat[1, 1] * pin.Y + mat[1, 2] * pin.Z;
			double z = mat[2, 0] * pin.X + mat[2, 1] * pin.Y + mat[2, 2] * pin.Z;
			pout.Set(x, y, z, false);
			return;
		}
		/// <summary>
		/// Rotate input point into output point
		/// 
		/// </summary>
		/// <param name="pin"></param>
		/// <param name="pout"></param>
		public void rotate(double[] pin, double[] pout) {
			pout[0] = MRTH[0, 0] * pin[0] + MRTH[0, 1] * pin[1] + MRTH[0, 2] * pin[2];
			pout[1] = MRTH[1, 0] * pin[0] + MRTH[1, 1] * pin[1] + MRTH[1, 2] * pin[2];
			pout[2] = MRTH[2, 0] * pin[0] + MRTH[2, 1] * pin[1] + MRTH[2, 2] * pin[2];
			return;
		}
		public void rotate(double x, double y, double z, ref Cartesian pout) {
			double ox, oy, oz;
			ox = MRTH[0, 0] * x + MRTH[0, 1] * y + MRTH[0, 2] * z;
			oy = MRTH[1, 0] * x + MRTH[1, 1] * y + MRTH[1, 2] * z;
			oz = MRTH[2, 0] * x + MRTH[2, 1] * y + MRTH[2, 2] * z;
			pout.Set(ox, oy, oz, false);
			return;
		}
		public void rotate(Cartesian vin, ref Cartesian vout) {
			double[] P = { vin.X, vin.Y, vin.Z };
			double[] Q = new double[3];
			rotate(P, Q);
			vout.Set(Q[0], Q[1], Q[2], false);
			return;
		}

		public void unrotate(double[] pin) {
			double x, y, z;
			unrotate(out x, out y, out z, pin);
			pin[0] = x;
			pin[1] = y;
			pin[2] = z;
		}
		public void unrotate(out double x, out double y, out double z, double[] pin) {
			x = MRTH[0, 0] * pin[0] + MRTH[1, 0] * pin[1] + MRTH[2, 0] * pin[2];
			y = MRTH[0, 1] * pin[0] + MRTH[1, 1] * pin[1] + MRTH[2, 1] * pin[2];
			z = MRTH[0, 2] * pin[0] + MRTH[1, 2] * pin[1] + MRTH[2, 2] * pin[2];
		}

		public void make_deccircle(double dec, double[] xarr, double[] yarr, double[] zarr, int n) {
			double ra, dra;
			double x, y, z;
			Cartesian pout = new Cartesian();
			ra = 0.0;
			dra = 360.0 / (n - 1);

			for (int i = 0; i < n; i++) {
				Cartesian.Radec2Xyz(ra, dec, out x, out y, out z);
				rotate(x, y, z, ref pout);
				xarr[i] = pout.X;
				yarr[i] = pout.Y;
				zarr[i] = pout.Z;
				ra += dra;
			}
		}
		public void make_racircle(double ra, int n, double[] xarr, double[] yarr, double[] zarr) {
			double dec, ddec;
			double x, y, z;
			Cartesian pout = new Cartesian();
			dec = -90.0;
			ddec = 180.0 / (n - 1);

			for (int i = 0; i < n; i++) {
				Cartesian.Radec2Xyz(ra, dec, out x, out y, out z);
				rotate(x, y, z, ref pout);
				xarr[i] = pout.X;
				yarr[i] = pout.Y;
				zarr[i] = pout.Z;
				dec += ddec;
			}
		}
		public static void getRotationMatrix(double[,] mat, Cartesian R, double theta) {
			double co, si;
			co = Math.Cos(theta);
			si = Math.Sin(theta);
			mat[0, 0] = co + (1.0 - co) * R.X * R.X;
			mat[0, 1] = (1.0 - co) * R.X * R.Y - si * R.Z;
			mat[0, 2] = (1.0 - co) * R.X * R.Z + si * R.Y;

			mat[1, 0] = (1.0 - co) * R.X * R.Y + si * R.Z;
			mat[1, 1] = co + (1.0 - co) * R.Y * R.Y;
			mat[1, 2] = (1.0 - co) * R.Y * R.Z - si * R.X;

			mat[2, 0] = (1.0 - co) * R.Z * R.X - si * R.Y;
			mat[2, 1] = (1.0 - co) * R.Z * R.Y + si * R.X;
			mat[2, 2] = co + (1.0 - co) * R.Z * R.Z;

		}

		//
		// axis and angle is an axis (x,y,z) followed by an angle 
		public static void getAxisandAngle(Cartesian P, Cartesian Q, ref Cartesian axis, out double angle) {
			Cartesian Pc = new Cartesian(P, false);
			Cartesian Qc = new Cartesian(Q, false);

			double n1, n2;
			n1 = Pc.Norm;
			n2 = Qc.Norm;

			angle = Math.Acos(Pc.Dot(Qc));

			Cartesian Tc;
			Tc = Pc.Cross(Qc, true);
			axis.Set(Tc.X, Tc.Y, Tc.Z, false);
		}
		//
		// M_stuff changes global rotation matrix
		//
		public void rotinto(double[,] mat, double[] P, double[] Q) {
			Cartesian Pc = new Cartesian(P[0], P[1], P[2], false);
			Cartesian Qc = new Cartesian(Q[0], Q[1], Q[2], false);
			double n1, n2;
			n1 = Pc.Norm;
			n2 = Qc.Norm;

			double theta = Math.Acos(Pc.Dot(Qc));
			Cartesian Tc = Pc.Cross(Qc, true);
			//Pc.crossMe(Qc);
			//Pc.Normalize();
			setrot(mat, Tc, theta);
		}
		public void M_set_rotinto(double[] P, double[] Q) {
			// make a matrixk that rotates P into Q
			Cartesian Pc = new Cartesian(P[0], P[1], P[2], false);
			Cartesian Qc = new Cartesian(Q[0], Q[1], Q[2], false);
			double n1, n2;
			n1 = Pc.Norm;
			n2 = Qc.Norm;

			double theta = Math.Acos(Pc.Dot(Qc));

			Cartesian Tc = Pc.Cross(Qc, true);
			//Pc.crossMe(Qc);
			//Pc.Normalize();
			M_setrot(Tc, theta); 
		}
		public void M_setrot(Cartesian R, double theta) {
			double co, si;
			co = Math.Cos(theta);
			si = Math.Sin(theta);
			MRTH[0, 0] = co + (1.0 - co) * R.X * R.X;
			MRTH[0, 1] = (1.0 - co) * R.X * R.Y - si * R.Z;
			MRTH[0, 2] = (1.0 - co) * R.X * R.Z + si * R.Y;

			MRTH[1, 0] = (1.0 - co) * R.X* R.Y + si * R.Z;
			MRTH[1, 1] = co + (1.0 - co) * R.Y * R.Y;
			MRTH[1, 2] = (1.0 - co) * R.Y * R.Z - si * R.X;

			MRTH[2, 0] = (1.0 - co) * R.Z * R.X - si * R.Y;
			MRTH[2, 1] = (1.0 - co) * R.Z * R.Y + si * R.X;
			MRTH[2, 2] = co + (1.0 - co) * R.Z * R.Z;


		}
		public void M_pre_multiply(double[,] mat) {
			double[,] tm = new double[3, 3];
			double s = 0.0;
			for (int i = 0; i < 3; i++) {
				for (int j = 0; j < 3; j++) {
					s = 0.0;
					for (int k = 0; k < 3; k++) {
						s += mat[i, k] * MRTH[k, j];
					}
					tm[i, j] = s;
				}
			}
			for (int i = 0; i < 3; i++) {
				for (int j = 0; j < 3; j++) {
					MRTH[i, j] = tm[i, j];
				}
			}

		}
		public void M_post_multiply(double[,] mat) {
			double[,] tm = new double[3, 3];
			double s = 0.0;
			for (int i = 0; i < 3; i++) {
				for (int j = 0; j < 3; j++) {
					s = 0.0;
					for (int k = 0; k < 3; k++) {
						s += MRTH[i, k] * mat[k, j];
					}
					tm[i, j] = s;
				}
			}
			for (int i = 0; i < 3; i++) {
				for (int j = 0; j < 3; j++) {
					MRTH[i, j] = tm[i, j];
				}
			}

		}
		public void setrot(double[,] mat, Cartesian R, double theta) {
			double co, si;
			co = Math.Cos(theta);
			si = Math.Sin(theta);
			mat[0, 0] = co + (1.0 - co) * R.X * R.X;
			mat[0, 1] = (1.0 - co) * R.X * R.Y - si * R.Z;
			mat[0, 2] = (1.0 - co) * R.X * R.Z + si * R.Y;

			mat[1, 0] = (1.0 - co) * R.X * R.Y + si * R.Z;
			mat[1, 1] = co + (1.0 - co) * R.Y * R.Y;
			mat[1, 2] = (1.0 - co) * R.Y * R.Z - si * R.X;

			mat[2, 0] = (1.0 - co) * R.Z * R.X - si * R.Y;
			mat[2, 1] = (1.0 - co) * R.Z * R.Y + si * R.X;
			mat[2, 2] = co + (1.0 - co) * R.Z * R.Z;


		}
	}
}
