#region Written by Tamas Budavari (budavari@jhu.edu)
/**
 * SphericalLib.Arc
 * 
 * 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: Arc.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.Collections.Generic;
#endregion

namespace Spherical
{
	/// <summary>
	/// An arclet is a segment of a circle on the unit sphere.
	/// </summary>
    [CLSCompliant(true)]
    [Serializable]
	public class Arc
    {
        #region Private members

        private Halfspace circle;
        private Cartesian point1;
        private Cartesian point2;
        private Cartesian? middle;
        private double? area;
        private double? length;

        #endregion
        #region Public properties

        public Halfspace Circle
        {
            get { return circle; }
        }

        public Cartesian Point1
        {
            get { return point1; }
        }

        public Cartesian Point2
        {
            get { return point2; }
        }

        public Cartesian Middle
        {
            get
            {
                if (middle == null)
                {
                    if (point1 == point2)
                    {
                        this.middle = point1;
                        this.area = 0;
                        this.length = 0;
                    }
                    else
                    {
                        this.middle = Arc.MiddlePoint(circle, point1, point2);
                    }
                }
                return middle.Value;
            }
        }

        public double Area
        {
            get
            {
                if (area == null)
                {
                    if (point1 == point2)
                    {
                        this.middle = point1;
                        this.area = 0;
                        this.length = 0;
                    }
                    else
                    {
                        this.area = Arc.AreaSemilune(circle, point1, point2, 
                            out this.length);
                    }
                }
                return area.Value;
            }
        }

        public double Length
        {
            get
            {
                if (length == null)
                {
                    if (point1 == point2)
                    {
                        this.middle = point1;
                        this.area = 0;
                        this.length = 0;
                    }
                    else
                    {
                        this.area = Arc.AreaSemilune(circle, point1, point2, 
                            out this.length);
                    }
                }
                return length.Value;
            }
        }

        #endregion
        #region Constructors

        /// <summary>
		/// Create an arc
		/// </summary>
		/// <param name="circle"></param>
		/// <param name="p1"></param>
		/// <param name="p2"></param>
		public Arc (Halfspace circle, Cartesian p1, Cartesian p2)
		{
			this.circle = circle;
			this.point1 = p1;
			this.point2 = p2;
        }

        /// <summary>
        /// Great circle arc
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        public Arc(Cartesian p1, Cartesian p2)
        {
            Cartesian x = p1.Cross(p2, true);
            this.circle = new Halfspace(x, 0);

            this.point1 = p1;
            this.point2 = p2;
            
            this.middle = p1.MiddlePoint(p2, true);
            this.area = 0;
            this.length = p1.AngleInDegree(p2);
        }
        
        /// <summary>
        /// Copy constructor
        /// </summary>
        /// <param name="a"></param>
        public Arc(Arc a)
        {
            this.point1 = a.point1;
            this.point2 = a.point2;
            this.circle = a.circle;
            this.middle = a.middle;
            this.length = a.length;
            this.area = a.area;
        }

        #endregion
        #region Area computation and verification

        /// <summary>
        /// Compute the area of the triangle defined by the arc and a given 
        /// point which is connected to the endpoints with great circles.
        /// </summary>
        /// <param name="p"></param>
        /// <returns></returns>
        public double CalcArea(Cartesian p, double tolerance)
        {
            if (p.Equals(point1) || p.Equals(point2) || p.Same(point1, tolerance) || p.Same(point2, tolerance)) 
                return this.area.Value;

            // 
            // 3 contributions: girard area, semilune and orance slice
            //
            //
            // Girard area: 0 if 2 points are the same
            double area = Cartesian.GirardArea(p, this.point1, this.point2);

            //
            // slice area is 2*angle
            if (point1.Same(point2.Mirrored(), tolerance))
            {
                Cartesian n = p.Cross(point1,true); // normal of plane with p
                area += 2 * (Math.PI-n.AngleInRadian(circle.Vector)) * Const.SquareRadian2SquareDegree;
                // set sign
                double det = Cartesian.TripleProduct(n,circle.Vector,point1);
                if (det < 0)
                    area *= -1;
            }
            //
            // semilune area
            area += this.Area;
            return area;
        }

        /// <summary>
        /// Verify if point are on circle within the given tolerance
        /// </summary>
        /// <param name="tolerance"></param>
        /// <returns></returns>
        public bool VerifyOkay(double tolerance)
        {
            if (Math.Abs(circle.Vector.Dot(point1) - circle.Cos0) < tolerance &&
                Math.Abs(circle.Vector.Dot(point2) - circle.Cos0) < tolerance)
                return true;
            else
                return false;
        }
       
        #endregion
        #region Static methods for area and middle point

        /// <summary>
		/// Compute area of the semilune
		/// </summary>
		/// <returns></returns>
        public static double AreaSemilune(Halfspace circle, Cartesian p1, Cartesian p2,
            out double? length)
		{
            double a, b, dx, dy, dz, dist, theta, beta, det, area, ratio, sinTheta;
            Cartesian center = circle.Vector;
            det = Cartesian.TripleProduct(center, p1, p2);

            // if great circle -> zero area
            if (Math.Abs(circle.Cos0) < double.Epsilon)
            {
                length = p1.AngleInDegree(p2);
                if (det < 0) length = 360 - length;
                return 0;
            }
            
            // otherwise
            if (circle.Cos0 < 0) center.Mirror();
            // angles
            theta = center.AngleInRadian(p1);               // great circle dist of p1 and center 
            sinTheta = Math.Sin(theta);
            dx = p1.X - p2.X;
            dy = p1.Y - p2.Y;
            dz = p1.Z - p2.Z;
            dist = 0.5 * Math.Sqrt(dx * dx + dy * dy + dz * dz);
            beta = 2*Math.Asin(dist);                       // great circ dist of p1 and p2
            ratio = Math.Tan(beta / 2) / Math.Tan(theta);
            if (ratio > 1) ratio = 1;
            a = Math.Asin(ratio);
            ratio = dist / Math.Sin(theta);
            if (ratio > 1) ratio = 1;
            b = 2*Math.Asin(ratio);
            length = b*Const.Radian2Degree*sinTheta;                 // output arc length
            // precise area
            area = 2 * a - b * Math.Cos(theta);
            if (det < 0)
            {
                area = 2 * Math.PI * (1 - Math.Cos(theta)) - area;
                length = 360 * sinTheta - length;
            }
            area *= Const.SquareRadian2SquareDegree;
            if (circle.Cos0 < 0) area *= -1;
            return area;
		}

        /// <summary>
        /// Compute middle point of arclet
        /// </summary>
        /// <returns></returns>
        public static Cartesian MiddlePoint(Halfspace circle, Cartesian p1, Cartesian p2)
        {
            Halfspace g = new Halfspace(new Cartesian(p2.X - p1.X, p2.Y - p1.Y, p2.Z - p1.Z, true), 0);
            return g.Root(circle);
        }

        /// <summary>
        /// Derive the minimal enclosing circle for a list of arcs
        /// </summary>
        /// <param name="arclist"></param>
        /// <returns></returns>
        public static Halfspace MinimalEnclosingCircle(Arc[] arclist)
        {
            Halfspace mec;
            Arc a0 = arclist[0];
            if (arclist.Length == 1)
            {
                mec = a0.circle;
            }
            else
            {
                List<Cartesian> pts = new List<Cartesian>(arclist.Length);
                foreach (Arc a in arclist)
                {
                    // end points may already be in list
                    if (!pts.Contains(a.point1)) pts.Add(a.point1);
                    if (!pts.Contains(a.point2)) pts.Add(a.point2);
                    // middle point cannot
                    if (a.circle.Cos0 > -double.Epsilon)
                        pts.Add(a.Middle);
                }
                mec = Cartesian.MinimalEnclosingCircle(pts);
            }
            return mec;
        }

        #endregion
        #region Check endpoints

        public bool EndpointsCoveredByHalfspace(Halfspace h)
        {
            return h.Covers(point1, Const.Tolerance) && h.Covers(point2, Const.Tolerance);
        }

        public bool EndpointsCoveredByConvex(Convex c)
        {
            return c.TrueForAll(this.EndpointsCoveredByHalfspace);
        }

        public bool EndpointsCoveredByRegion(Region r)
        {
            return r.Exists(this.EndpointsCoveredByConvex);
        }

        #endregion
        #region Morphological operations

        /// <summary>
        /// Create a convex with given width along the arc
        /// </summary>
        /// <param name="arcmin"></param>
        /// <returns></returns>
        public Convex WideConvex(double arcmin)
        {
            Convex c = new Convex(4);
            arcmin /= 2;

            // one side
            Halfspace one = this.circle; // copy it
            one.Grow(arcmin);            // grow it
            c.Add(one);
            
            // other side
            Halfspace two = this.circle; // copy it
            two.Invert();                // invert it
            two.Grow(arcmin);            // grow
            c.Add(two);
            
            // ends
            c.Add(new Halfspace(circle.Vector.Cross(point1, true), 0));
            c.Add(new Halfspace(point2.Cross(circle.Vector, true), 0));

            return c;
        }

        #endregion
        #region ToString and revision

        /// <summary>
		/// String representation
		/// </summary>
		/// <returns></returns>
		public override string ToString()
		{
#if false
            return String.Format("  {0}\n\t start:  {1}\n\t end:    {2}", 
                this.circle, this.p1, this.p2);
#else
            string r = "";
            r += String.Format("   Area: {0}\n", Area);
            r += String.Format("   Length: {0}\n", Length);
            r += String.Format("   Circle: {0}\n", Circle);
            r += String.Format("   Start:  {0}\n", Point1);
            r += String.Format("   End:    {0}\n", Point2);
            r += String.Format("   Middle: {0}\n", Middle);
            r += String.Format("   -- RADEC --\n");
            r += String.Format("   Circle: {0} {1} / {2}\n", Circle.Vector.RA, Circle.Vector.Dec, Circle.RadiusInRadian);
            r += String.Format("   Start  point: {0} {1}\n", Point1.RA, Point1.Dec);
            r += String.Format("   End    point: {0} {1}\n", Point2.RA, Point2.Dec);
            r += String.Format("   Middle point: {0} {1}\n", Middle.RA, Middle.Dec);
            return r;
#endif
		}

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

#region Revision History
/* Revision History

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

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

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

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

        Revision 1.3  2006/02/28 16:26:59  budavari
        Major revision: public properties

        Revision 1.2  2006/01/24 20:53:40  budavari
        Arc has length

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

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

        Revision 1.20  2005/12/29 14:22:53  budavari
        Check argument of asin

        Revision 1.19  2005/12/20 02:14:40  budavari
        Build convex from an arc with given width

        Revision 1.18  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.17  2005/12/04 16:50:55  budavari
        Small features and fixes. Also added Outline class - not function, yet

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

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

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

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

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

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

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

        Revision 1.9  2005/11/29 15:18:31  budavari
        Separate MinimalEnclosingCircle() method for Arc[]

        Revision 1.8  2005/11/29 13:57:24  budavari
        Added constructor for great circle arc (takes two points)

        Revision 1.7  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.6  2005/11/22 14:45:03  budavari
        Small fixes and more samples/tests

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

        Revision 1.4  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.3  2005/11/15 23:23:39  budavari
        Format string changed

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

        Revision 1.1  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

*/
#endregion