#region Written by Tamas Budavari (budavari@jhu.edu)
/**
 * SphericalLib.Cartesian
 * 
 * SphericalLib is a .NET library to deal with geometry on the unit sphere.
 * It is intented to be used for complex survey footprint descriptions in
 * standalone applications as well as in web services and perhaps Yukon.
 * 
 * Written by Tamas Budavari (budavari@jhu.edu)
 * See bottom of file for revision history
 * 
 * Current revision:
 *   ID:          $Id: Cartesian.cs,v 1.1.1.1 2006/11/14 02:56:39 gupchup Exp $
 *   Revision:    $Revision: 1.1.1.1 $
 *   Date:        $Date: 2006/11/14 02:56:39 $
 */
using System;
using System.Xml.Serialization;
using System.Collections.Generic;
#endregion

namespace Spherical
{
	/// <summary>
	/// 3D cartesian vectors with members x, y, z.
    /// <remarks>It is a mutable value type.</remarks>
	/// </summary>
    [Serializable, XmlRoot(Namespace = Const.NameSpace, IsNullable=false)]
    [CLSCompliant(true)]
	public struct Cartesian
    {
        #region Private members x, y, z
        private double x,y,z;
        #endregion
        #region Public properties
        
        [XmlAttribute]
        public double X
        {
            get { return x; }
            set { x = value; }
        }

        [XmlAttribute]
        public double Y
        {
            get { return y; }
            set { y = value; }
        }

        [XmlAttribute]
        public double Z
        {
            get { return z; }
            set { z = value; }
        }

        /// <summary>
        /// Compute length of x,y,z vector
        /// </summary>
        public double Norm
        {
            get { return Math.Sqrt(x * x + y * y + z * z); }
        }
        
		/// <summary>
		/// Get property of angular coordinate RA, 
		/// which is computed on the fly
		/// </summary>
		public double RA
		{
			get 
			{ 
				double ra;
				if (Math.Abs(z) < 1) 
				{
					ra = Math.Atan2(y,x) * Const.Radian2Degree;
					if (ra<0) ra += 360;
				}
				else 
				{
					ra = 0;
				}
				return ra;
			}
		}

		/// <summary>
		/// Get property of angular coordinate Dec, 
		/// which is computed on the fly
		/// </summary>
		public double Dec
		{
			get 
			{	
				if (Math.Abs(z) < 1) 
				{
					return Math.Asin(z) * Const.Radian2Degree;
				}
				else
				{
					return 90.0 * Math.Sign(z);
				}
			}
		}

        #endregion
		#region Constructors of various kinds
		public Cartesian(Cartesian other, bool normalize) {
			this.x = other.x;
			this.y = other.y;
			this.z = other.z;
			if (normalize) {
				this.Normalize();
			}
		}

		/// <summary>
		/// Create a point using a vector (x,y,z), which should get 
		/// normalized unless you really know what you are doing
		/// </summary>
		/// <param name="x"></param>
		/// <param name="y"></param>
		/// <param name="z"></param>
		/// <param name="normalize"></param>
		public Cartesian (double x, double y, double z, bool normalize)
		{
			this.x = x;
			this.y = y;
			this.z = z;
			if (normalize)
				this.Normalize();
		}

		/// <summary>
		/// Create a point from the angular coordinates (RA,Dec)
		/// x, y, z = cos cos, cos sin, sin
		/// </summary>
		/// <param name="ra"></param>
		/// <param name="dec"></param>
		public Cartesian (double ra, double dec)
		{
			Cartesian.Radec2Xyz (ra, dec, out x, out y, out z);
		}

		#endregion
		#region Set methods
		public void Set(Cartesian v, bool normalize) {
			x = v.x;
			y = v.y;
			z = v.z;
			if (normalize) {
				this.Normalize();
			}
		}

		/// <summary>
		/// Set point by a direction vector (x,y,z),
		/// which should get normalized unless you 
		/// realy know what you are doing
		/// </summary>
		/// <param name="x"></param>
		/// <param name="y"></param>
		/// <param name="z"></param>
		/// <param name="normalize"></param>
		public void Set (double x, double y, double z, bool normalize)
		{
			this.x = x;
			this.y = y;
			this.z = z;
			if (normalize)
			{
				this.Normalize();
			}
		}

		/// <summary>
		/// Set point from angular position (RA,Dec)
		/// </summary>
		/// <param name="ra"></param>
		/// <param name="dec"></param>
		public void Set(double ra, double dec)
		{
			Cartesian.Radec2Xyz (ra, dec, out x, out y, out z);
		}

        /// <summary>
        /// Set point to be the middle point of two given points
        /// and normalize if needed
        /// </summary>
        /// <param name="p"></param>
        /// <param name="q"></param>
        /// <param name="normalize"></param>
        public void SetMiddlePoint(Cartesian p, Cartesian q, bool normalize)
        {
            this.x = p.x + q.x;
            this.y = p.y + q.y;
            this.z = p.z + q.z;
            if (normalize)
            {
                this.Normalize();
            }
        }

		#endregion
		#region Operations on self

		/// <summary>
		/// Normalize the point
		/// </summary>
		/// <returns></returns>
		public double Normalize()
		{
            double L = this.Norm;
			this.x /= L;
			this.y /= L;
			this.z /= L;
			return L;
		}

        /// <summary>
		/// Mirror the point: v *= -1
		/// </summary>
		public void Mirror()
		{
			this.x *= -1;
			this.y *= -1;
			this.z *= -1;
		}

		/// <summary>
		/// Return mirrored point: -v
		/// </summary>
		/// <returns></returns>
		public Cartesian Mirrored()
		{
			return new Cartesian(-x,-y,-z,false);
		}

        /// <summary>
        /// Multiply components with given number
        /// </summary>
        /// <param name="s"></param>
        public void Scale(double s)
        {
            this.x *= s;
            this.y *= s;
            this.z *= s;
        }

        /// <summary>
        /// Return a vector that is scaled by the given number
        /// </summary>
        /// <param name="s"></param>
        /// <returns></returns>
        public Cartesian Scaled(double s)
        {
            Cartesian r = this;
            r.Scale(s);
            return r;
        }

        #endregion
		#region Geometrical relations

		/// <summary>
		/// Return middle point: m = (p+q) / |p+q|
		/// </summary>
		/// <param name="p"></param>
		/// <returns></returns>
		public Cartesian MiddlePoint (Cartesian p, bool normalize)
		{
			return new Cartesian (x+p.x, y+p.y, z+p.z, normalize);
		}

		/// <summary>
		/// Cosine of angle to an another point.
		/// It is the 3D dot product of the 2 vectors
		/// </summary>
		/// <param name="p"></param>
		/// <returns></returns>
		public double Dot (Cartesian p)
		{
			return x*p.x + y*p.y + z*p.z;
		}

		/// <summary>
		/// Add to another vector, result is a new vector
		/// </summary>
		/// <param name="p"></param>
		/// <returns></returns>
		public Cartesian Add(Cartesian p) 
        {
            return new Cartesian(x + p.x, y + p.y, z + p.z, false);
		}

		/// <summary>
		/// Subtract another vector, result is a new vector
		/// </summary>
		/// <param name="p"></param>
		/// <returns></returns>
		public Cartesian Sub(Cartesian p) 
        {
			return new Cartesian(x - p.x, y - p.y, z - p.z, false);
		}

		/// <summary>
		/// Cross product with another direction/point, result is again on the unit sphere
		/// </summary>
		/// <param name="p"></param>
		/// <returns></returns>
		public Cartesian Cross (Cartesian p, bool normalize)
		{
			return new Cartesian
				(y*p.z-z*p.y, z*p.x-x*p.z, x*p.y-y*p.x, normalize);
		}

		/// <summary>
		/// Angle to an another point
		/// </summary>
		/// <param name="p"></param>
		/// <returns></returns>
		public double AngleInRadian(Cartesian p)
		{
			double a, d, dx,dy,dz;
			dx = x - p.x;
			dy = y - p.y;
			dz = z - p.z;
			d = 0.5 * Math.Sqrt(dx*dx + dy*dy + dz*dz);
			if (d<1) a = 2*Math.Asin(d);
			else     a = Math.PI;
			return a;
		}

		/// <summary>
		/// Angle to another point
		/// </summary>
		/// <param name="p"></param>
		/// <returns></returns>
		public double AngleInDegree(Cartesian p)
		{
			return AngleInRadian(p) * Const.Radian2Degree;
		}

		/// <summary>
		/// Angle to another point
		/// </summary>
		/// <param name="p"></param>
		/// <returns></returns>
		public double AngleInArcmin(Cartesian p)
		{
			return AngleInRadian(p) * Const.Radian2Arcmin;
		}

		/// <summary>
		/// Area of spherical triangle defined this and other two point.
		/// </summary>
		/// <param name="p"></param>
		/// <param name="q"></param>
		/// <returns></returns>
		public double GirardArea(Cartesian p, Cartesian q)
		{
			return Cartesian.GirardArea(this,p,q);
		}

        /// <summary>
        /// Check if Point is identical to another within a given tolerance
        /// </summary>
        /// <param name="p"></param>
        /// <param name="epsilon"></param>
        /// <returns></returns>
        public bool Same(Cartesian p, double epsilon)
        {
            bool eq = false;
            if (Math.Abs(x - p.x) < epsilon
             && Math.Abs(y - p.y) < epsilon
             && Math.Abs(z - p.z) < epsilon)
                eq = true;
            return eq;
        }

        /// <summary>
        /// Check if Point is identical to another within Const.Epsilon
        /// </summary>
        /// <param name="p"></param>
        /// <returns></returns>
        public bool Same(Cartesian p)
        {
            return Same(p, Const.Epsilon);
        }

        /// <summary>
        /// Test whether or not this vector is the opposite
        /// of another vector.
        /// Opposite means same magnitude, opposite direction.
        /// </summary>
        /// <param name="p">The other vector</param>
        /// <returns>true if other vector is close enough to 
        /// being opposite, false otherwise</returns>
        public bool Opposite(Cartesian p, double epsilon)
        {
            if (Math.Abs(p.x + x) > epsilon 
             || Math.Abs(p.y + y) > epsilon
             || Math.Abs(p.z + z) > epsilon)
                return false;
            return true;
        }

        /// <summary>
        /// Test whether or not this vector is the opposite
        /// of another vector within Const.Epsilon
        /// Opposite means same magnitude, opposite direction.
        /// </summary>
        /// <param name="p">The other vector</param>
        /// <returns>true if other vector is close enough to 
        /// being opposite, false otherwise</returns>
        public bool Opposite(Cartesian p)
        {
            return this.Opposite(p, Const.Epsilon);
        }

#if false
        /// <summary>
        /// Check if given halfspace includes this point - uses given tolerance
        /// </summary>
        /// <param name="h"></param>
        /// <returns></returns>
        public bool CoveredByHalfspace(Halfspace h, double tolerance)
        {
            return x * h.x + y * h.y + z * h.z > h.cos0 - tolerance;
        }

        /// <summary>
        /// Check if given halfspace includes this point - uses default tolerance
        /// </summary>
        /// <param name="h"></param>
        /// <returns></returns>
        public bool CoveredByHalfspace(Halfspace h)
        {
            return x * h.x + y * h.y + z * h.z > h.cos0 - Const.Tolerance;
        }

        /// <summary>
        /// Check if given convex includes this point - uses default tolerance
        /// </summary>
        /// <param name="c"></param>
        /// <returns></returns>
        public bool CoveredByConvex(Convex c)
        {
            return c.HalfspaceList.TrueForAll(this.CoveredByHalfspace);
        }

        /// <summary>
        /// Check if given region includes this point - uses default tolerance
        /// </summary>
        /// <param name="r"></param>
        /// <returns></returns>
        public bool CoveredByRegion(Region r)
        {
            return r.ConvexList.Exists(this.CoveredByConvex);
        }
#endif

		#endregion
		#region Static methods

        /// <summary>
        /// Triple product of three vectors is the volume of the parallelepipedon.
        /// The sign shows the winding of the points.
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <param name="p3"></param>
        /// <returns></returns>
        public static double TripleProduct(Cartesian p1, Cartesian p2, Cartesian p3)
        {
            double x, y, z; // cross product of 1 and 2
            x = p1.y * p2.z - p1.z * p2.y;
            y = p1.z * p2.x - p1.x * p2.z;
            z = p1.x * p2.y - p1.y * p2.x;
            return x * p3.x + y * p3.y + z * p3.z;
        }

        /// <summary>
        /// Convert (x,y,z) to RA,DEC.
        /// (x,y,z) must be normalized.
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="z"></param>
        /// <param name="ra">out in degrees</param>
        /// <param name="dec">out in degrees</param>
        public static void Xyz2Radec (double x, double y, double z, 
            out double ra, out double dec)
        {
            dec = Math.Asin(z) * Const.Radian2Degree;
            double cd = Math.Cos(dec * Const.Degree2Radian);
            if (cd > Const.Epsilon || cd < -Const.Epsilon)
                if (y > Const.Epsilon || y < -Const.Epsilon)
                    if (y < 0.0)
                        ra = 360 - Math.Acos(x / cd) * Const.Radian2Degree;
                    else
                        ra = Math.Acos(x / cd) * Const.Radian2Degree;
                else
                    ra = (x < 0.0 ? 180.0 : 0.0);
            else
                ra = 0.0;
            return;
        }

        /// <summary>
        /// Convert (RA, DEC) to (x, y, z).
        /// RA, DEC must be well defined.
        /// </summary>
        /// <param name="ra">in degrees</param>
        /// <param name="dec">in degrees</param>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="z"></param>
        static public void Radec2Xyz (double ra, double dec, 
            out double x, out double y, out double z)
        {
            double diff;
            double cd = Math.Cos(dec * Const.Degree2Radian);
            diff = 90.0 - dec;
            /// First, compute Z, consider cases, where declination is almost
            /// +/- 90 degrees
            if (diff < Const.Epsilon && diff > -Const.Epsilon)
            {
                x = 0.0;
                y = 0.0;
                z = 1.0;
                return;
            }
            diff = -90.0 - dec;
            if (diff < Const.Epsilon && diff > -Const.Epsilon)
            {
                x = 0.0;
                y = 0.0;
                z = -1.0;
                return;
            }
            z = Math.Sin(dec * Const.Degree2Radian);
            ///
            /// If we get here, then 
            /// at least z is not singular
            ///
            double quadrant;
            double qint;
            int iint;
            quadrant = ra / 90.0; // how close is it to an integer?
            // if quadrant is (almost) an integer, force x, y to particular
            // values of quad:
            // quad,   (x,y)
            // 0       (1,0)
            // 1,      (0,1)
            // 2,      (-1,0)
            // 3,      (0,-1)
            // q>3, make q = q mod 4, and reduce to above
            // q mod 4 should be 0.
            qint = (double)((int)quadrant);
            if (Math.Abs(qint - quadrant) < Const.Epsilon)
            {
                iint = (int)qint;
                iint %= 4;
                if (iint < 0) iint += 4;
                switch (iint)
                {
                    case 0:
                        x = cd;
                        y = 0.0;
                        break;
                    case 1:
                        x = 0.0;
                        y = cd;
                        break;
                    case 2:
                        x = -cd;
                        y = 0.0;
                        break;
                    case 3:
                    default:
                        x = 0.0;
                        y = -cd;
                        break;
                }
            }
            else
            {
                x = Math.Cos(ra * Const.Degree2Radian) * cd;
                y = Math.Sin(ra * Const.Degree2Radian) * cd;
            }
            return;
        }

        /// <summary>
        /// Static method to convert point to (RA,Dec)
        /// </summary>
        /// <param name="p"></param>
        /// <param name="ra"></param>
        /// <param name="dec"></param>
        public static void Cartesian2Radec(Cartesian p, out double ra, out double dec)
        {
            Xyz2Radec(p.x, p.y, p.z, out ra, out dec);
        }

        /// <summary>
        /// Select and return a subset of points. Copies references of points.
        /// </summary>
        /// <param name="points"></param>
        /// <param name="idx"></param>
        /// <returns></returns>
        public static Cartesian[] SelectPoints(List<Cartesian> points, List<int> idx)
        {
            Cartesian[] sub = new Cartesian[idx.Count];
            for (int i = 0; i < sub.Length; i++)
                sub[i] = points[idx[i]];
            return sub;
        }

        /// <summary>
        /// Select and return a subset of root points. Copies references of points.
        /// </summary>
        /// <param name="points"></param>
        /// <param name="idx"></param>
        /// <returns></returns>
        internal static Cartesian[] SelectPoints(List<Root> roots, List<int> idx)
        {
            Cartesian[] sub = new Cartesian[idx.Count];
            for (int i = 0; i < sub.Length; i++)
                sub[i] = roots[idx[i]].Point;
            return sub;
        }

        /// <summary>
        /// Compute the center of mass of the points in the list
        /// </summary>
        /// <returns></returns>
        public static Cartesian CenterOfMass(List<Cartesian> plist, bool normalize)
        {
            Cartesian m = new Cartesian(0, 0, 0, false);
            foreach (Cartesian p in plist)
            {
                m.x += p.x;
                m.y += p.y;
                m.z += p.z;
            }
            if (normalize)
                m.Normalize();
            return m;
        }

        /// <summary>
        /// Brute-force implementation of the minimal enclosing circle
        /// takes three points.
        /// </summary>
        /// <returns></returns>
        public static Halfspace MinimalEnclosingCircle(Cartesian p0, Cartesian p1, Cartesian p2)
        {
            Halfspace mec = new Halfspace(new Cartesian(1, 0, 0, false), -1);
            Halfspace c;
            
            // check diagonals
            // 0
            c = new Halfspace(p0, p1);  
            if (c.Cos0 > mec.Cos0 && c.Covers(p2,Const.Tolerance)) 
                mec = c;
            // 1
            c = new Halfspace(p0, p2);
            if (c.Cos0 > mec.Cos0 && c.Covers(p1, Const.Tolerance))
                mec = c;
            // 2
            c = new Halfspace(p1, p2);
            if (c.Cos0 > mec.Cos0 && c.Covers(p0, Const.Tolerance))
                mec = c;
            
            // check triangle
            c = new Halfspace(p0, p1, p2); 
            if (c.Cos0 > mec.Cos0) mec = c;

            return mec;
        }

        /// <summary>
        /// Brute-force implementation of the minimal enclosing circle, 
        /// which scales with N^4 (means slow ;-)
        /// </summary>
        /// <returns></returns>
        public static Halfspace MinimalEnclosingCircle(List<Cartesian> plist)
        {
            Halfspace mec = new Halfspace(new Cartesian(1, 0, 0, false), -1);
            double angle, minangle;
            bool found = false;
            minangle = Math.PI;

            // if a single point -> single point
            if (plist.Count == 1)
                return new Halfspace(plist[0], 1);

            // check diagonals
            Cartesian Pi, Pj, c;
            for (int i = 0; i < plist.Count - 1; i++)
            {
                Pi = plist[i];
                for (int j = i + 1; j < plist.Count; j++)
                {
                    Pj = plist[j];
                    // middle point (normalized, of course)
                    if (Pi.Same(Pj.Mirrored(), Const.Tolerance))
                        c = Cartesian.CenterOfMass(plist, true);
                    else
                        c = plist[i].MiddlePoint(plist[j], true); 
                    Halfspace h = new Halfspace (c, c.Dot(Pj));
                    angle = h.RadiusInRadian;
                    if (angle < minangle && h.Covers(plist,Const.Tolerance))
                    {
                        found = true;
                        mec = h;
                        minangle = angle;
                    }
                }
            }
            // check triangles
            for (int i = 0; i < plist.Count - 1; i++)
            {
                for (int j = i + 1; j < plist.Count; j++)
                {
                    for (int k = 0; k < plist.Count; k++)
                    {
                        if (k == i || k == j) continue;
                        Halfspace h = new Halfspace(plist[i], plist[j], plist[k]);
                        angle = h.RadiusInRadian;
                        if (angle < minangle && h.Covers(plist,Const.Tolerance))
                        {
                            found = true;
                            mec = h;
                            minangle = angle;
                        }
                        else
                        {
                            h.Invert();
                            angle = h.RadiusInRadian;
                            if (angle < minangle && h.Covers(plist,Const.Tolerance))
                            {
                                found = true;
                                mec = h;
                                minangle = angle;
                            }
                        }
                    }
                }
            }
            if (!found)
                throw new ArgumentException("MinimalEnclosingCircle(): invalid input?");

            return mec;
        }

        #endregion
        #region Girard area computation

        /// <summary>
        /// Compute area of triangle defined by 3 points connected with great circles.
        /// </summary>
        /// <param name="p"></param>
        /// <param name="q"></param>
        /// <param name="r"></param>
        /// <returns></returns>
        public static double GirardArea(Cartesian p1, Cartesian p2, Cartesian p3)
        {
            //
            // Girard's theorem
            //
            double area = 0, det, a1 = 0, a2 = 0, a3 = 0, s = 0, dx, dy, dz;
            Cartesian v1 = new Cartesian(p1.x - p3.x, p1.y - p3.y, p1.z - p3.z, false);
            Cartesian v2 = new Cartesian(p2.x - p3.x, p2.y - p3.y, p2.z - p3.z, false);
            Cartesian v3 = p3;
            int branch = 0;
            // determinant is the volume
            det = Cartesian.TripleProduct(v1,v2,v3);

            // if zero then area = 0, else...
            if (Math.Abs(det) > 1e-15)
            {
                branch |= 1;
                a1 = 2 * Math.Asin(0.5 * Math.Sqrt(v2.x * v2.x + v2.y * v2.y + v2.z * v2.z));
                a2 = 2 * Math.Asin(0.5 * Math.Sqrt(v1.x * v1.x + v1.y * v1.y + v1.z * v1.z));
                dx = v1.x - v2.x;
                dy = v1.y - v2.y;
                dz = v1.z - v2.z;
                a3 = 2 * Math.Asin(0.5 * Math.Sqrt(dx * dx + dy * dy + dz * dz));

                s = 0.5 * (a1 + a2 + a3);

#if true
                area = -Math.PI
                    + 2 * Math.Asin(Math.Sqrt(Math.Sin(s - a2) * Math.Sin(s - a3) / (Math.Sin(a2) * Math.Sin(a3))))
                    + 2 * Math.Asin(Math.Sqrt(Math.Sin(s - a1) * Math.Sin(s - a3) / (Math.Sin(a1) * Math.Sin(a3))))
                    + 2 * Math.Asin(Math.Sqrt(Math.Sin(s - a2) * Math.Sin(s - a1) / (Math.Sin(a2) * Math.Sin(a1))))
                    ;
                branch |= 2;
#else
                if ( (s - a1) > 1e-8 && (s - a2) > 1e-8 && (s - a3) > 1e-8 )
                {
                    area = -Math.PI
                        + 2 * Math.Asin(Math.Sqrt(Math.Sin(s - a2) * Math.Sin(s - a3) / (Math.Sin(a2) * Math.Sin(a3))))
                        + 2 * Math.Asin(Math.Sqrt(Math.Sin(s - a1) * Math.Sin(s - a3) / (Math.Sin(a1) * Math.Sin(a3))))
                        + 2 * Math.Asin(Math.Sqrt(Math.Sin(s - a2) * Math.Sin(s - a1) / (Math.Sin(a2) * Math.Sin(a1))))
                        ;
                    branch |= 2;
                }
                else
                {
                    area =
                        + Math.Asin( det / (Math.Sin(a2) * Math.Sin(a3)))
                        + Math.Asin( det / (Math.Sin(a1) * Math.Sin(a3)))
                        + Math.Asin( det / (Math.Sin(a1) * Math.Sin(a2)));
                    branch |= 4;
                }
#endif
                area *= Const.SquareRadian2SquareDegree;

                // it should not be negative (maybe rounding error)
                if (area < 0)
                {
#if false
                    throw new ArgumentException
                        (String.Format("Point.GirardArea(): negative area: {0}", area));
#else
                    area = 0;
#endif
                }
                if (det < 0) area *= -1;
            }
#if false
            Console.Out.WriteLine("G -> det: {0}", det);
            Console.Out.WriteLine("G -> branch: {0}", branch);
            Console.Out.WriteLine("G -> a1 a2 a3 s: {0} {1} {2} {3}",
                a1 * Const.Radian2Degree, a2 * Const.Radian2Degree, a3 * Const.Radian2Degree, s * Const.Radian2Degree);
            Console.Out.WriteLine("G -> s1 s2 s3: {0} {1} {2}",
                +2 * Math.Asin(Math.Sqrt(Math.Sin(s - a2) * Math.Sin(s - a3) / (Math.Sin(a2) * Math.Sin(a3))))
                   * Const.Radian2Degree,
                +2 * Math.Asin(Math.Sqrt(Math.Sin(s - a1) * Math.Sin(s - a3) / (Math.Sin(a1) * Math.Sin(a3))))
                   * Const.Radian2Degree,
                +2 * Math.Asin(Math.Sqrt(Math.Sin(s - a2) * Math.Sin(s - a1) / (Math.Sin(a2) * Math.Sin(a1))))
                   * Const.Radian2Degree
               );
            Console.Out.WriteLine("G -> n1 n2 n3: {0} {1} {2}",
                 +Math.Asin(det / (Math.Sin(a2) * Math.Sin(a3))) * Const.Radian2Degree,
                 +Math.Asin(det / (Math.Sin(a1) * Math.Sin(a3))) * Const.Radian2Degree,
                 +Math.Asin(det / (Math.Sin(a1) * Math.Sin(a2))) * Const.Radian2Degree
               );
            Console.Out.WriteLine("G -> area: {0}", area);
#endif
            return area;
        }

        #endregion
        #region Predefined readonly points

        public static readonly Cartesian NaN = new Cartesian(double.NaN,double.NaN,double.NaN,false);
        public static readonly Cartesian Xaxis = new Cartesian(1, 0, 0, false);
        public static readonly Cartesian Yaxis = new Cartesian(0, 1, 0, false);
        public static readonly Cartesian Zaxis = new Cartesian(0, 0, 1, false); 

		#endregion
        #region Verification

        /// <summary>
        /// Check if point is really on the unit sphere:
        /// Length is close to one within the tolerance limit
        /// </summary>
        /// <param name="tolerance"></param>
        /// <returns></returns>
        public bool VerifyOkay(double tolerance)
        {
            return Math.Abs(1 - x * x - y * y - z * z) < tolerance;
        }

        #endregion
        #region Equality related methods...

        public override bool Equals(object right)
        {
            if (right == null)
                return false;

            if (object.ReferenceEquals(this, right))
                return true;

            if (this.GetType() != right.GetType())
                return false;

            Cartesian p = (Cartesian)right;
            return (x == p.x) && (y == p.y) && (z == p.z);
        }

        public override int GetHashCode()
        {
            return base.GetHashCode();
        }

        public static bool operator ==(Cartesian left, Cartesian right)
        {
            return (left.x == right.x) && (left.y == right.y) && (left.z == right.z);
        }

        public static bool operator !=(Cartesian left, Cartesian right)
        {
            return (left.x != right.x) || (left.y != right.y) || (left.z != right.z);
        }

        #endregion
        #region Parsing, ToString and revision

        /// <summary>
        /// Parse an input string of "x y z" into a Point. 
        /// </summary>
        /// <param name="s"></param>
        /// <returns></returns>
        public static Cartesian Parse(string repr, bool normalize)
        {
            char[] sep = new char[] { ' ' };

            // only spaces after this (hopefully)
            string proc = repr.Trim().Replace(',', ' ').Replace("\r\n", " ").Replace('\n', ' ').Replace('\t', ' ');

            // replace duplicate spaces with a single one
            string trim = "";
            while (trim != proc)
            {
                trim = proc;
                proc = trim.Replace("  ", " ");
            }

            string[] t = trim.Split(sep, 3);
            double x = double.Parse(t[0]);
            double y = double.Parse(t[1]);
            double z = double.Parse(t[2]);
            return new Cartesian(x, y, z, normalize);
        }
        /// <summary>
		/// String representation is just the 3 coordinates: x y z
		/// </summary>
		/// <returns></returns>
		public override string ToString()
		{
			return String.Format("{0:R} {1:R} {2:R}", x, y, z);
		}

		/// <summary>
		/// Revision from CVS
		/// </summary>
		public static string Revision
		{
			get	{ return "$Revision: 1.1.1.1 $"; }
        }
        #endregion
    }
}

#region Revision History
/* Revision History

        $Log: Cartesian.cs,v $
        Revision 1.1.1.1  2006/11/14 02:56:39  gupchup
        no message

        Revision 1.12  2006/05/22 18:55:54  budavari
        BinaryFormatter now works with Region and co

        Revision 1.11  2006/04/19 16:51:01  budavari
        EVerything is now CLS compliant

        Revision 1.10  2006/04/10 19:58:26  budavari
        Explicit not nullable attrib added

        Revision 1.9  2006/03/31 15:27:00  budavari
        small additions and typo fix

        Revision 1.8  2006/03/09 20:30:07  budavari
        Equality checks for value types

        Revision 1.7  2006/03/07 14:54:37  budavari
        Added set to public properties...

        Revision 1.6  2006/03/03 21:21:15  gyuri
        Added a Set and a copying constructor

        Revision 1.5  2006/02/28 16:25:28  budavari
        Major revision: struct for speed

        Revision 1.4  2006/01/31 15:47:28  budavari
        Added missing point checks in MEC(1,2,3)

        Revision 1.3  2006/01/30 19:09:50  budavari
        MinimalEnclosingCircle w/ 3 points and Halfspace.ctor w/ 2 points

        Revision 1.2  2006/01/27 21:49:59  budavari
        Need to add fast MEC for 3 points

        Revision 1.1  2006/01/24 15:56:05  budavari
        Re-structured directory layout to accommodate more projects

        Revision 1.12  2006/01/20 17:39:03  budavari
        MiddlePoint() has optional normalize

        Revision 1.11  2006/01/20 16:43:22  budavari
        Added Set and SetMiddlePoint

        Revision 1.10  2005/12/29 14:26:17  budavari
        Optional normalize in Parse()

        Revision 1.9  2005/12/14 22:38:16  budavari
        Epsilon changes in Same()

        Revision 1.8  2005/12/08 14:05:35  gyuri
        Added Add and Sub methods which create new instances of Cartesians similarly to Cross

        Revision 1.7  2005/12/05 21:43:55  budavari
        Working Outline implementation :-) Hurray! Except for the fly in the soup:
        I did find an example, Samples.Saw3(), which makes the code throw up as
        it turns the outline into a bowtie. Coming up next...

        Revision 1.6  2005/12/04 16:50:55  budavari
        Small features and fixes. Also added Outline class - not function, yet

        Revision 1.5  2005/12/02 18:03:43  budavari
        Code changes to use the new TripleProduct()

        Revision 1.4  2005/12/02 15:37:49  budavari
        Static Cartesian.TripleProduct() added

        Revision 1.3  2005/12/02 14:49:53  budavari
        New radec <-> xyz conversion a'la HTM code

        Revision 1.2  2005/12/02 14:40:13  budavari
        Okay... struct Area is gone - not needed really (I think ;-)

        Revision 1.1  2005/12/01 22:34:57  budavari
        Area struct added and fully functional

        Revision 1.16  2005/12/01 19:02:20  budavari
        Middle of adding an Area class

        Revision 1.15  2005/11/30 20:50:28  budavari
        Cosmetics... #region specifications etc added

        Revision 1.14  2005/11/30 19:50:37  budavari
        Dilate() methods added. Also explicit tolerance required everywhere for clarity

        Revision 1.13  2005/11/30 15:32:37  budavari
        Header changes

        Revision 1.12  2005/11/30 15:20:10  budavari
        Exact translation of Alex's Girard formula

        Revision 1.11  2005/11/24 16:27:16  budavari
        Region areas compared to SQL code and releated fixes of degenerate arcs.
        Looking good except maybe for one patch (convex 0) of the 'footprint',
        which has 2 degenerate triangles and patch area changes when shifting
        pivot point.

        Revision 1.10  2005/11/22 14:45:03  budavari
        Small fixes and more samples/tests

        Revision 1.9  2005/11/21 14:26:28  budavari
        Fixed smaller bugs after adding NUnit tests

        Revision 1.8  2005/11/18 05:58:09  budavari
        Working on boolean operations on regions...

        Revision 1.7  2005/11/17 01:58:54  budavari
        Done with the convex stuff... Hurray!!! It can even deal with the "orange
        slices" to some extent: the problem is the area - coming up next...

        Revision 1.6  2005/11/15 22:20:20  budavari
        Using Generics - almost there

        Revision 1.5  2005/11/15 15:25:36  budavari
        Everything is a class, things are coming together for area computation, etc. Will moved to C# 2.0

        Revision 1.4  2005/11/05 17:47:44  budavari
        Fixed a bug in Convex caused by renaming things and did some tidying up...

        Revision 1.3  2005/11/04 22:49:36  budavari
        Re-write to have Point and Halfspace as struct again and other tinkerings. XML is still fine...

        Revision 1.2  2005/11/03 15:22:01  budavari
        Fairly major revision to make sure the objects serialize to XML nicely:
          - Added separate containers of halfspaces and convexes
          - Customized XML output for compactness and friendly proxy interface
          - Implemented simple .Parse() methods taking cartesion representation

        Revision 1.1.1.1  2005/10/28 14:13:43  budavari
        Initial import of the SphericalLib library. It is pretty complete in terms
        of the classes but lacking the simplication and other advanced methods.
		
*/
#endregion