using System;
using System.IO;
using System.Configuration;
using System.Collections.Generic;
using System.Text;
using System.Runtime.Serialization;
using System.Runtime.Serialization.Formatters.Binary;

using NUnit.Framework;
using Spherical;
using Spherical.IO;
using Spherical.Visualizer;

namespace Spherical.Test
{
    [TestFixture]
    public class HtmTests
    {
        [Test]
        public void Name()
        {
            RandomPoint rnd = new RandomPoint();
            Cartesian p = rnd.Next();
            Console.Out.WriteLine("P: {0}", p);
            int level = 24;

            string s = HTM.Trixel.CartesianToString(p.X, p.Y, p.Z, level);
            Console.Out.WriteLine("Level {0}. Name: {1}", level, s);

            long id = HTM.Trixel.CartesianToHid (p.X, p.Y, p.Z, level);
            Console.Out.WriteLine("Level {0}.  HID: {1}", level, id);

            Cartesian a,b,c;
            HTM.Trixel.ToTriangle(id, out a, out b, out c);
            Console.Out.WriteLine("Triangle:");
            Console.Out.WriteLine(a);
            Console.Out.WriteLine(b);
            Console.Out.WriteLine(c);
            Console.Out.WriteLine("Area: {0:R}", Cartesian.GirardArea(a,b,c));

        }
    }

    [TestFixture]
    public class IO
    {
        string conn = "Data Source=sdssdb001; Initial Catalog=bestdr4; User=test; Password=riuzg54; Timeout=60;";

        [Test]
        public void USA()
        {
            string filename = @"C:\Tamas\Visual Studio 2005\Projects\Spherical\Geometry\USA\st24_d00.dat";
            DatReader dat = new DatReader(filename);
            Cartesian mid = new Cartesian();
            Region res = new Region();

            List<Cartesian> pts = new List<Cartesian>();
            //dat.NextPointList(pts, 10, out mid);
            //pts.Clear();

            for (int i = 0; i < 3; i++)
            {
                Console.Out.Write("Starting {0}...", i);
                int id = dat.NextPointList(pts, 100, out mid);
                Outline o = new Outline(pts);
                Region r = o.Region(Const.Tolerance);
                Console.Out.WriteLine(" done: {0}", id);
                res.Union(r);
                pts.Clear();
            }
            RegionImage.SaveGif(@"C:\Temp\rState.gif", res, mid.RA, mid.Dec, 300);
        }

        [Test]
        public void BinIO()
        {
            Region w = Samples.SdssDr4();
            FileInfo dr4 = new FileInfo(@"C:\Temp\dr4.bin");

            // write
            using (FileStream dr4stream = dr4.Create())
            {
                w.Simplify(); // for sorting 
                foreach (Convex c in w.ConvexList)
                    c.Simplified = false;
                Binary.Writer(dr4stream, w);
                w.Simplify();
                Binary.Writer(dr4stream, w);
            }

            // read simple region nd test
            using (FileStream dr4stream = dr4.OpenRead())
            {
                Region r = Binary.Reader(dr4stream);
                Assert.AreEqual(w.ConvexList.Count, r.ConvexList.Count, "# of convex (plain)");
                for (int i = 0; i < w.ConvexList.Count; i++)
                {
                    Assert.AreEqual
                        (w.ConvexList[i].HalfspaceList.Count, r.ConvexList[i].HalfspaceList.Count,
                         String.Format("# of halfspaces in convex {0}", i));
                    for (int j = 0; j < w.ConvexList[i].HalfspaceList.Count; j++)
                    {
                        Assert.IsTrue(w.ConvexList[i].HalfspaceList[j].Same
                            (r.ConvexList[i].HalfspaceList[j], double.Epsilon),
                            String.Format("convex '{0}' halfspace '{1}' different", i, j));
                    }
                }

                // read full region and test
                r = Binary.Reader(dr4stream);
                Assert.AreEqual(w.ConvexList.Count, r.ConvexList.Count, "# of convex (full)");
                Assert.AreEqual(w.Area.Value, r.Area.Value, "Area");
            }
        }

        [Test]
        public void SqlDr4Region6599()
        {
            SqlIO sql = new SqlIO(conn);
            Region r = sql.ReadRegion(6599);
            sql.Close();
            sql.Dispose();
            r.SimpleSimplify(Const.Tolerance);

            //Console.Out.WriteLine("{0}:\n{1}", r.name, r);
            Assert.AreEqual(12, r.ConvexList.Count, "# of convexes in region");
            for (int i = 0; i < 12; i++)
                Assert.AreEqual(4, r.ConvexList[i].HalfspaceList.Count,
                    String.Format("# of halfspace in convex '{0}'", i));

            r.Simplify(false, Const.Tolerance, Const.MaximumIteration, false);
            Console.Out.WriteLine("Area: {0}", r.Area);
        }

        [Test]
        public void SqlDr4RegionString6599()
        {
            SqlIO sql = new SqlIO(conn);
            Region r = sql.ReadRegionString(6599);
            sql.Close();
            sql.Dispose();
            r.SimpleSimplify(Const.Tolerance);

            //Console.Out.WriteLine("{0}:\n{1}", r.name, r);
            Assert.AreEqual(12, r.ConvexList.Count, "# of convexes in region");
            for (int i = 0; i < 12; i++)
                Assert.AreEqual(4, r.ConvexList[i].HalfspaceList.Count,
                    String.Format("# of halfspace in convex '{0}'", i));

            r.Simplify(false, Const.Tolerance, Const.MaximumIteration, false);
            Console.Out.WriteLine("Area: {0}", r.Area);
        }

        [Test]
        public void SqlDr4RegionString()
        {
            SqlIO sql = new SqlIO(conn);
            List<long> id = sql.QueryRegionId("select regionid from region");
            double area;
            string type;
            byte ismask;
            int counter = 0, problem = 0;
            foreach (long i in id)
            {
                Region r = sql.ReadRegionString(i, out area, out type, out ismask);
                //Console.Out.WriteLine("{0} \t{1} \t{2} \t{3} \t{4}", i, area, r.name, type, ismask);
                counter++;
                r.Simplify(false, Const.Tolerance, Const.MaximumIteration, false);
                double delta = area - r.Area.Value;
                if (Math.Abs(delta) / area > 1e-4)
                {
                    Console.Error.WriteLine("{0} \t{1} \t{2} \t{3} \t{4} \t{5} ",
                        i, area, r.Area.Value, delta, delta / area, r.Name);
                    problem++;
                }
            }
            Console.Error.WriteLine("# of regions with different area: {0} / {1}", problem, counter);
        }

        [Test]
        public void SqlDr4Region()
        {
            SqlIO sql = new SqlIO(conn);
            List<long> id = sql.QueryRegionId("select regionid from region");
            double area;
            string type;
            byte ismask;
            int counter = 0, problem = 0;
            foreach (long i in id)
            {
                sql.ReadRegionString(i, out area, out type, out ismask);
                Region r = sql.ReadRegion(i);
                //Console.Out.WriteLine("{0}\n{1}", r.name, r);
                //Console.Out.WriteLine("{0} \t{1} \t{2} \t{3} \t{4}", i, area, r.name, type, ismask);
                counter++;
                r.Simplify(false, Const.Tolerance, Const.MaximumIteration, false);
                double delta = area - r.Area.Value;
                Console.Out.WriteLine("{0} \t{1} \t{2} \t{3} \t{4} \t{5} ",
                        i, area, r.Area, delta, delta / area, r.Name);
                if (Math.Abs(delta) / area > 1e-4)
                {
                    Console.Error.WriteLine("{0} \t{1} \t{2} \t{3} \t{4} \t{5} ",
                        i, area, r.Area, delta, delta / area, r.Name);
                    problem++;
                }
                //Assert.AreEqual(area, s.area, 1e-5, String.Format("area of {0}", r.name));
            }
            Console.Error.WriteLine("# of regions with different area: {0} / {1}", problem, counter);
        }

        [Test]
        public void SqlDr4Region41458problem()
        {
            SqlIO sql = new SqlIO(conn);
            Region r = sql.ReadRegion(41458);
            sql.Close();
            sql.Dispose();
            //r.SimpleSimplify(Const.Tolerance);

            Console.Out.WriteLine("{0}:\n{1}", r.Name, r);
            r.Simplify(false, Const.Tolerance, Const.MaximumIteration, false);
            Console.Out.WriteLine("Area: {0}", r.Area);
        }

        [Test]
        public void SqlRegionWrite()
        {
            using(SqlIO sql = new SqlIO("Data Source=sdssdb001; Initial Catalog=GalexCollab; Integrated Security=SSPI;"))
            {
                Region r = new Region(Samples.Clover(false), false);
                r.Simplify(false, Const.Tolerance, Const.MaximumIteration, false);
                long id = sql.Write(r);
                Console.Out.WriteLine("Region ID: {0}", id);
                sql.Close();
            }            
        }

        [Test]
        public void BinaryFormatter()
        {
            FileInfo tempfile = new FileInfo(@"C:\Temp\regionBinaryFormatter.bin");
            double area;

            using (FileStream sw = tempfile.Create())
            {
                Region r = new Region(Samples.Clover(false), false);
                r.Simplify(false, Const.Tolerance, Const.MaximumIteration, false);
                area = r.Area.Value;
                Console.Out.WriteLine("Area before: {0}", area);
                Binary.Writer(sw, r);
            }

            using (FileStream rd = tempfile.OpenRead())
            {
                Region r = Binary.Reader(rd);
                Console.Out.WriteLine("Area after: {0}", area);
                Assert.AreEqual(area, r.Area.Value, double.Epsilon, "area after roundtrip");
            }

            tempfile.Delete();
        }

        [Test]
        public void POV()
        {
            double area;
            Region r = Samples.Box(out area);
#if false
            FileInfo f = new FileInfo(@"c:\temp\spherical.pov");
            PovWriter pov = new PovWriter(f.Create());
#else
            PovWriter pov = new PovWriter(Console.Out);
#endif
            pov.Write(r, new Cartesian(0, 0).Scaled(1.05));
            pov.Close();
#if false
            Console.Out.WriteLine(f.OpenText().ReadToEnd());
#endif
        }
    }

    [TestFixture]
    public class Outlines
    {
        [Test]
        public void Concave()
        {
            double e = 0.01;
            Cartesian a = new Cartesian(e, 0);
            Cartesian b = new Cartesian(e, e);
            Cartesian c = new Cartesian(0, e);
            Cartesian d = new Cartesian(0, 2 * e);
            bool bcon = Outline.IsConcave(a, b, c, Const.Tolerance);
            bool ccon = Outline.IsConcave(b, c, d, Const.Tolerance);
            bool icon = Outline.IsConcave(c, b, a, Const.Tolerance);
            Console.Out.WriteLine("B is concave? {0}", bcon);
            Console.Out.WriteLine("C is concave? {0}", ccon);
            Console.Out.WriteLine("Brev is concave? {0}", icon);
            Assert.IsFalse(bcon, "B");
            Assert.IsTrue(ccon, "C");
            Assert.IsTrue(icon, "Brev");
        }

        [Test]
        public void Greek()
        {
            Outline o = Samples.Greek(false);
            Region r = o.Region();
            Console.Out.WriteLine(r);
            Console.Out.WriteLine("Area: {0}", r.Area);
            //r.Simplify(false);
            Assert.AreEqual(5, r.ConvexList.Count, "convex");
            Assert.IsTrue(r.Area < Const.WholeSphereInSquareDegree, "less area than whole sky");
            RegionImage.SaveGif(@"C:\Temp\rGreek.gif", r, 0, 0, 200);
        }

        [Test]
        public void GreekInverse()
        {
            Outline o = Samples.Greek(true);
            Region r = o.Region();
            Console.Out.WriteLine(r);
            Console.Out.WriteLine("Area: {0}", r.Area);
            //Assert.AreEqual(5, r.ConvexList.Count, "convex");
            Assert.IsTrue(r.Area < Const.WholeSphereInSquareDegree, "less area than whole sky");
            RegionImage.SaveGif(@"C:\Temp\rGreekInv.gif", r, 0, 0, 200);
        }

        [Test]
        public void Star()
        {
            Outline o = Samples.Star(false);
            Region r = o.Region();
            Console.Out.WriteLine(r);
            Console.Out.WriteLine("Area: {0}", r.Area);
            Assert.AreEqual(4, r.ConvexList.Count, "convex");
            Assert.IsTrue(r.Area < Const.WholeSphereInSquareDegree, "less area than whole sky");
            RegionImage.SaveGif(@"C:\Temp\rStar.gif", r, 0, 0, 300);
        }
        
        [Test]
        public void StarInverse()
        {
            Outline o = Samples.Star(true);
            Region r = o.Region();
            Console.Out.WriteLine(r);
            Console.Out.WriteLine("Area: {0}", r.Area);
            Assert.AreEqual(8, r.ConvexList.Count, "convex");
            Assert.IsTrue(r.Area < Const.WholeSphereInSquareDegree, "less area than whole sky");
            RegionImage.SaveGif(@"C:\Temp\rStarInv.gif", r, 0, 0, 300);
        }

        [Test]
        public void Saw()
        {
            Outline o = Samples.Saw(false);
            Region r = o.Region();
            Console.Out.WriteLine(r);
            Console.Out.WriteLine("Area: {0}", r.Area);
            Assert.AreEqual(3, r.ConvexList.Count, "convex");
            RegionImage.SaveGif(@"C:\Temp\rSaw.gif", r, 0, 0, 1500);
        }

        [Test]
        public void Arrow()
        {
            Outline o = Samples.Arrow(false);
            Region r = o.Region();
            Console.Out.WriteLine("Region:\n{0}", r);
            Console.Out.WriteLine("Area: {0}", r.Area);
            Assert.AreEqual(2, r.ConvexList.Count, "convex");
            Assert.AreEqual(3, r.ConvexList[0].HalfspaceList.Count, "1st convex");
            Assert.AreEqual(3, r.ConvexList[1].HalfspaceList.Count, "2nd convex");
            Console.Out.WriteLine("Area: {0}", r.Area);
            RegionImage.SaveGif(@"C:\Temp\rArrow.gif", r, 0, 0, 1000);
        }

        [Test]
        public void Pacman()
        {
            Outline o = Samples.Pacman(false);
            Region r = o.Region();
            Console.Out.WriteLine("Region:\n{0}", r);
            Console.Out.WriteLine("Area: {0}", r.Area);
            Assert.AreEqual(2, r.ConvexList.Count, "2 convex");
            Assert.IsTrue(r.Area < Const.WholeSphereInSquareDegree, "less area than whole sky");
            RegionImage.SaveGif(@"C:\Temp\rPacman.gif", r, 0, 0, 1000);
        }
        
        [Test]
        public void PacmanInverse()
        {
            Outline o = Samples.Pacman(true);
            //Console.Out.WriteLine("Outline:\n{0}", o);
            Region r = o.Region();
            Console.Out.WriteLine("Region:\n{0}", r);
            foreach (Convex c in r.ConvexList) Console.Out.WriteLine("Convex area: {0}", c.Area);
            Console.Out.WriteLine("Area: {0}", r.Area);
            Assert.IsTrue(r.Area < Const.WholeSphereInSquareDegree, "less area than whole sky");
            RegionImage.SaveGif(@"C:\Temp\rPacmanInv.gif", r, 0, 0, 1000);
        }

        [Test]
        public void Square()
        {
            Outline o = Samples.Square(false);
            Region r = o.Region(Const.Tolerance);
            Console.Out.WriteLine("Area: {0}", r.Area);
            Assert.AreEqual(1, r.ConvexList.Count, "one convex");
            Assert.AreEqual(4, r.ConvexList[0].HalfspaceList.Count, "four halfspaces");
            Assert.AreEqual(0.1 * 0.1, r.Area.Value, 1e-5, "area");
            RegionImage.SaveGif(@"C:\Temp\rSquare.gif", r, 0, 0, 500);
        }

        [Test]
        public void SquareInverse()
        {
            Outline o = Samples.Square(true);
            Region r = o.Region();
            Console.Out.WriteLine("# of Convex: {0}", r.ConvexList.Count);
            Console.Out.WriteLine("Area: {0}", r.Area);
            foreach (Convex c in r.ConvexList) Console.Out.WriteLine("Convex area: {0}", c.Area);
            Assert.IsTrue(r.Area < Const.WholeSphereInSquareDegree, "less area than whole sky");
            RegionImage.SaveGif(@"C:\Temp\rSquareInv.gif", r, 0, 0, 500);
        }
        
    }

    [TestFixture]
    public class Areas
    {
        [Test]
        public void RegionTest()
        {
            double ref_area;
            Region r = Samples.Box(out ref_area);
            r.Simplify();
            Console.Out.WriteLine("Super Box area: {0}", r.Area);
            double area = r.CalcArea(new Cartesian(120, 0));
            Console.Out.WriteLine(" shd be same as: {0}", area);

            Assert.AreEqual(4, r.ConvexList.Count, "four convexes");
            Assert.AreEqual(ref_area, r.Area.Value, 1e-5, "area");
        }

        [Test]
        public void ConvexTest()
        {
            double ref_area;
            Convex c = Samples.OrangeSlices(out ref_area);
            c.Simplify();
            Console.Out.WriteLine("Super OraneSlices area: {0}", c.Area);
            double area = c.CalcArea(new Cartesian(120, 0));
            Console.Out.WriteLine("  which shd be same as: {0}", area);

            Assert.AreEqual(2, c.HalfspaceList.Count, "two halfspaces");
            Assert.AreEqual(ref_area, c.Area.Value, 1e-5, "area");
        }

        [Test]
        public void ConvexIntersect()
        {
            Convex c1 = new Convex(new Halfspace(180, 0, 60));
            Convex c2 = new Convex(new Halfspace(181, 0, 60));
            Convex cc1 = new Convex(c1);

            c1.Intersect(c2);
            c1.Simplify();

            c2.Simplify();
            cc1.Simplify();
            
            cc1.SmartIntersect(c2,Const.Tolerance,Const.MaximumIteration);
            //cc1.Simplify();

            Console.Out.WriteLine("Areas: {0} ? {1}", c1.Area, cc1.Area);
            Assert.AreEqual(c1.Area.Value, cc1.Area.Value, 1e-7, "areas");
        }

        [Test]
        public void ConvexIntersectEmpty()
        {
            Convex c1 = new Convex(new Halfspace(180, 0, 60));
            Convex c2 = new Convex(new Halfspace(183, 0, 60));
            Convex cc1 = new Convex(c1);

            c1.Intersect(c2);
            c1.Simplify();
            c2.Simplify();
            cc1.Simplify();
            cc1.SmartIntersect(c2, Const.Tolerance, Const.MaximumIteration);
            //cc1.Simplify();

            Console.Out.WriteLine("Areas: {0} ? {1}", c1.Area, cc1.Area);
            Assert.AreEqual(c1.Area.Value, cc1.Area.Value, 1e-7, "areas");
        }

        [Test]
        public void RegionIntersect()
        {
            Region r1 = new Region(180, 0, 60);
            Region r2 = new Region(181, 0, 60);

            Region ri = r1.Intersect(r2);
            ri.Simplify();
            r1.Simplify();
            r2.Simplify();
            Region si = r1.SmartIntersect(r2, Const.Tolerance, Const.MaximumIteration);

            Console.Out.WriteLine(ri);

            Console.Out.WriteLine("Areas: {0} ? {1}", ri.Area, si.Area);
            Assert.AreEqual(ri.Area.Value, si.Area.Value, 1e-7, "areas");
        }

        [Test]
        public void RegionIntersectEmpty()
        {
            Region r1 = new Region(180, 0, 60);
            Region r2 = new Region(183, 0, 60);
            Region ri = r1.Intersect(r2);
            ri.Simplify();
            r1.Simplify();
            r2.Simplify();
            Region si = r1.SmartIntersect(r2, Const.Tolerance, Const.MaximumIteration);

            Console.Out.WriteLine("Areas: {0} ? {1}", ri.Area, si.Area);
            Assert.AreEqual(ri.Area.Value, si.Area.Value, 1e-7, "areas");
        }

        [Test]
        public void DisjointConvexes()
        {
            Region r = new Region(2);
            r.Add(new Convex(new Halfspace(180, 0, 60)));
            r.Add(new Convex(new Halfspace(181, 0, 60)));
            Console.Out.WriteLine("Before: {0}", r);
            r.Simplify();
            Console.Out.WriteLine("After: {0}", r);
            Assert.AreEqual(2, r.ConvexList.Count, "2 disjoint convexes");
            Assert.IsTrue(r.Area.Value < 6, "area less than 2*Pi");
        }

        [Test]
        public void AreaSortingConvex()
        {
            Region r = new Region(4);
            r.Add(new Convex(new Halfspace(180, 0, 10)));
            r.Add(new Convex(new Halfspace(180, 0, 20)));
            r.Add(new Convex(new Halfspace(180, 0, 30)));
            r.Add(new Convex(new Halfspace(180, 0, 10)));
            foreach (Convex c in r.ConvexList) c.Simplify();

            r.Sort(Convex.ComparisonAreaDescending);
            Console.Out.Write("Areas in desc order: ");
            foreach (Convex c in r.ConvexList)
                Console.Out.Write(" {0}", c.Area);
            Console.Out.WriteLine();
            for (int i = 1; i < r.ConvexList.Count; i++)
                Assert.IsTrue(r.ConvexList[i - 1].Area.Value >= r.ConvexList[i].Area.Value,
                    String.Format("area order not desc: {0} ? {1}", i - 1, i));
        }

    }

    [TestFixture]
    public class Basics
    {
        [Test]
        public void EquatorRing()
        {
            Convex c = new Convex();
            c.Add(new Halfspace(1, 0, 0, false, -0.01));
            c.Add(new Halfspace(-1, 0, 0, false, -0.01));
            c.Simplify();
            Assert.AreEqual(2, c.PatchList.Count, "2 patches");
            Console.Out.WriteLine("Area: {0}", c.Area);
            foreach (Patch p in c.PatchList)
                Console.Out.WriteLine(" - patch area: {0}", p.Area);
        }

        [Test]
        public void Moon()
        {
            Convex c = Samples.Moon();
            c.Simplify(Const.Tolerance);
            Console.Out.WriteLine("Area: {0}", c.Area);
        }

        [Test]
        public void SimplifyDisjoint()
        {
            Region r1 = Region.Parse("REGION CONVEX 0.9 0.2 0 0.999",true);
            Region r2 = Region.Parse("REGION CONVEX 1 0 0 0.998",true);
            Region i = r1.Intersect(r2);
            Console.Out.WriteLine("Before: {0}", i);
            i.Simplify();
            Console.Out.WriteLine("After: {0}", i);
            Assert.AreEqual(0, i.ConvexList.Count, "no convex");
        }

        [Test]
        public void Circle()
        {
            Region r = new Region(0, 0, 60);
            r.Simplify();
            Console.Out.WriteLine("Region area: {0}\n", r.Area);
            Assert.AreEqual(1, r.ConvexList.Count, "1 convex");
            Assert.AreEqual(1, r.ConvexList[0].PatchList.Count, "1 patch");
            Assert.AreEqual(1, r.ConvexList[0].PatchList[0].ArcList.Count, "1 arc");
            Console.Out.WriteLine("Arc: {0}", r.ConvexList[0].PatchList[0].ArcList[0]);
            Console.Out.WriteLine("Patch length: {0}", r.ConvexList[0].PatchList[0].Length);
        }

        [Test]
        public void CircleUnion()
        {
            Region r = new Region(0, 0, 60);
            r.Union(new Convex(new Halfspace(0, 1, 60)));
            Console.Out.WriteLine("Before: {0}", r);
            r.Simplify(false);
            Console.Out.WriteLine("Area of two: {0}", r.Area);
            r.Simplify();
            Console.Out.WriteLine("Disjoint area: {0}", r.Area);
        }

        [Test]
        public void MasterCard()
        {
            Region u2circ = new Region(2);
            u2circ.Add(new Convex(new Halfspace(0, 0, 60)));
            u2circ.ConvexList[0].Add(new Halfspace(1, 0, 60).Inverse());
            u2circ.Add(new Convex(new Halfspace(1, 0, 60)));

            Convex i2circ = new Convex(2);
            i2circ.Add(new Halfspace(0, 0, 60));
            i2circ.Add(new Halfspace(1, 0, 60));
            Region r = u2circ.Difference(i2circ);

            Console.Out.WriteLine("Union: {0}\n", u2circ);
            Console.Out.WriteLine("Intersection: {0}\n", i2circ);
            Console.Out.WriteLine("Raw Diff: {0}\n", r);

            //foreach (Convex c in r.convexList) c.SimpleSimplify(Const.Tolerance);
            //Console.Out.WriteLine("Convex simple: {0}\n", r);
            r.SimpleSimplify(Const.Tolerance);
            Console.Out.WriteLine("Region simple: {0}\n", r);
        }

        [Test]
        public void ArcLength()
        {
            // Great circle tests
            Arc agc = new Arc(new Cartesian(0, 0), new Cartesian(90, 0));
            Console.Out.WriteLine("Arc GC90 length: {0}", agc.Length);
            Assert.AreEqual(0, agc.Area, "semilune area of great circle arc shd be 0");
            Assert.AreEqual(90, agc.Length, 1e-7, "length of great circle arc");

            // Small circle tests
            Halfspace h = new Halfspace(new Cartesian(0, 0, 1, false), Math.Sqrt(3)/2);
            Halfspace g = new Halfspace(new Cartesian(1, 0, 0, false), 0.0);
            Cartesian a, b;
            h.Roots(g, out a, out b);
            Arc ha = new Arc(h, a, b);
            Arc hi = new Arc(h, b, a);
            Arc ga = new Arc(g, a, b);
            Arc gi = new Arc(g, b, a);
            Console.Out.WriteLine("Arc Ha length: {0}", ha.Length);
            Console.Out.WriteLine("Arc Hi length: {0}", hi.Length);
            Console.Out.WriteLine("Arc Ga length: {0}", ga.Length);
            Console.Out.WriteLine("Arc Gi length: {0}", gi.Length);
            Assert.AreEqual(90, ha.Length, 1e-7, "small circle arc");
            Assert.AreEqual(90, hi.Length, 1e-7, "small circle arc complementer");
            Assert.AreEqual(60, ga.Length, 1e-7, "great circle arc");
            Assert.AreEqual(300, gi.Length, 1e-7, "great circle arc complementer");
        }

        [Test]
        public void RegionSimpleSimplify()
        {
            Convex c = new Convex();
            c.Add(new Halfspace(0, 0, 60));
            c.Add(new Halfspace(0, 1, 60));
            c.Add(new Halfspace(0.1, 0, 60));
            Convex k = new Convex();
            k.Add(new Halfspace(0.1, 0, 60));
            k.Add(new Halfspace(0, 1, 60).Inverse());
            k.Add(new Halfspace(0, 0, 60));
            Region r = new Region();
            r.Add(c);
            r.Add(k);
            Console.Out.WriteLine("Before: {0}", r);
            r.SimpleSimplify(Const.Tolerance);
            Console.Out.WriteLine("After: {0}", r);
            Assert.AreEqual(1, r.ConvexList.Count, "1 convex");
            Assert.AreEqual(2, r.ConvexList[0].HalfspaceList.Count, "2 halfspace");
        }

        [Test]
        public void RegionSimpleSimplify2()
        {
            Convex c = new Convex();
            c.Add(new Halfspace(0, 0, 60));
            c.Add(new Halfspace(0, 1, 60));
            c.Add(new Halfspace(0.1, 0, 60));
            Convex k = new Convex();
            k.Add(new Halfspace(0.1, 0, 60));
            k.Add(new Halfspace(0, 1, 60));
            k.Add(new Halfspace(0, 0, 60));
            Region r = new Region();
            r.Add(c);
            r.Add(k);
            Console.Out.WriteLine("Before: {0}", r);
            r.SimpleSimplify(Const.Tolerance);
            Console.Out.WriteLine("After: {0}", r);
            Assert.AreEqual(1, r.ConvexList.Count, "1 convex");
            Assert.AreEqual(3, r.ConvexList[0].HalfspaceList.Count, "3 halfspace");
        }

        [Test]
        public void RegionSimpleSimplify3()
        {
            double a;
            Region r = Samples.Box(out a);
            Region rr = new Region(r);
            r.Union(r.ConvexList[0]);
            Console.Out.WriteLine("Before: {0}", r);
            r.SimpleSimplify(Const.Tolerance);
            Console.Out.WriteLine("After: {0}", r);
            Assert.AreEqual(rr.ConvexList.Count, r.ConvexList.Count, "# of convexes");
        }

        [Test]
        public void Dr4Simplify()
        {
            Region r = Region.Parse(Samples.dr4repr, true);
            r.SimpleSimplify(Const.Tolerance);
            r.Simplify(true, Const.Tolerance, Const.MaximumIteration, true);
        }

        [Test]
        public void HalfspaceSort()
        {
            Convex c = new Convex(3);
            c.Add(new Halfspace(0, 0, 100));
            c.Add(new Halfspace(0, 0, 10));
            c.Add(new Halfspace(0, 0, 50));
            Console.Out.WriteLine(c);
            c.Sort(Halfspace.ComparisonRadiusAscending);
            Console.Out.WriteLine(c);
        }

        [Test]
        public void SqlDr4Region44675problem()
        {
#if true
            string repr = @"CONVEX
  -0.858720331930258 0.511384230539128 0.0329478418943097 0.999661878489243
  -0.0564770543906784 -0.00494110200852136 -0.998391670557349 0
  0.848093384554534 -0.529075305017794 -0.028582034482598 -0.999661878489243
  0.857721947847771 -0.514111977414764 -0.00139099205869182 -0.999661878489243
  0.862921419021446 -0.500686996330638 -0.0684043587749869 -0.999661878489243
  0.873697811880468 -0.485337704293679 -0.0331578996051319 -0.999661878489243
";
            Convex c = Convex.Parse(repr, true);
#else
            SqlIO sql = new SqlIO(conn);
            Region r = sql.ReadRegion(44675);
            sql.Close(); sql.Dispose();
            Convex c = r.convexList[0];
#endif
            Console.Out.WriteLine("Before:\n{0}", c);
            c.Simplify();
            Console.Out.WriteLine("After 1st simplify:\n{0}", c);
            int np1 = c.PatchList.Count;

            c.Simplify();
            Console.Out.WriteLine("After 2nd simplify:\n{0}", c);
            int np2 = c.PatchList.Count;

            Console.Out.WriteLine("PatchList counts: {0} ? {1}", np1, np2);
            Assert.AreEqual(np1, np2, "# of patches after 2nd simplification");
        }
    
        [Test]
        public void Simple()
        {
            Convex c = new Convex(2);
            c.Add(new Halfspace(1, 0, 0, false, 0));
            c.Add(new Halfspace(0, 1, 0, false, 0.99));
            c.Simplify(Const.Tolerance);
            Console.Out.WriteLine("Area: {0}", c.Area);
            Assert.AreEqual(1, c.PatchList.Count, "one patch");
            Assert.IsTrue(c.Area > 0, "area is not > 0");
        }

        [Test]
        public void ArcConvex()
        {
            Arc a = new Arc(new Cartesian(180, 0), new Cartesian(181, 0));
            Convex c = a.WideConvex(6);
            Console.Out.WriteLine("Arc: {0}", a);
            Console.Out.WriteLine("Convex:\n{0}", c);
            c.Simplify();
            Assert.AreEqual(0.1, c.Area.Value, 1e-7, "area");
        }

        [Test]
        public void Empty2()
        {
            Convex c = new Convex("empty", 2);
            Halfspace h = new Halfspace(1, 0, 0, false, 0.2);
            c.Add(h);
            c.Add(h.Inverse());
            bool exist = c.SimpleSimplify(Const.Tolerance);
            Assert.AreEqual(false, exist, "nothing left");
            Assert.AreEqual(0, c.HalfspaceList.Count, "no halfspace");
        }

        [Test]
        public void Empty()
        {
            Convex c = new Convex("empty", 2);
            c.Add(new Halfspace(1, 0, 0, false, 0.2));
            c.Add(new Halfspace(-1, 0, 0, false, 0.2));
            bool exist = c.SimpleSimplify(Const.Tolerance);
            Assert.AreEqual(false, exist, "nothing left");
            Assert.AreEqual(0, c.HalfspaceList.Count, "no halfspace");
        }

        [Test]
        public void WholeSky()
        {
            Convex c = new Convex(new Halfspace(1, 0, 0, false, -1));
            c.SimpleSimplify(Const.Tolerance);
            Assert.AreEqual(1, c.HalfspaceList.Count, "one halfspace");
            c.Simplify(Const.Tolerance);
            Assert.AreEqual(1, c.PatchList.Count, "one patch");
            Assert.AreEqual(1, c.PatchList[0].ArcList.Count, "one arc");
            Console.Out.WriteLine("All sky area: {0}", c.Area);
            Assert.AreEqual(Const.WholeSphereInSquareDegree, c.Area.Value, 1e-5, "area");

            double area = c.CalcArea(new Cartesian(120, 0));
            Console.Out.WriteLine(" shd be same as: {0}", area);
            Assert.AreEqual(Const.WholeSphereInSquareDegree, area, 1e-5, "area");
        }

        [Test]
        public void Degenerate()
        {
            double e = 1e-4 ;
            Console.Out.WriteLine("Epsilon: {0}", e);

            List<Cartesian> pts = new List<Cartesian>(3);
            pts.Add(new Cartesian(0,-40));
            pts.Add(new Cartesian(e, 0));
            pts.Add(new Cartesian(0, 40));
            Convex c = new Convex(pts, PointOrder.Safe);
            Assert.AreEqual(3, c.HalfspaceList.Count, "three halfspaces");

            c.Simplify();
            Assert.AreEqual(1, c.PatchList.Count, "one patch");
            Assert.AreEqual(3, c.PatchList[0].ArcList.Count, "three arcs");
            Assert.AreEqual(3, c.HalfspaceList.Count, "still 3 halfspaces");

            double area, est_area, girard_area;
            girard_area = Cartesian.GirardArea(pts[0], pts[1], pts[2]);
            est_area = 2* e * Const.Degree2Radian * Const.SquareRadian2SquareDegree * 80 / 180 * 0.8;

            area = c.PatchList[0].Area;
            Console.Out.WriteLine("Area (def): {0}  ??  Girard: {1} ?? Est: {2}", area, girard_area, est_area);
            Assert.AreEqual(girard_area, area, 1e-5, "area w/ default center vs girard");

            RandomPoint rnd = new RandomPoint();
            Cartesian center = rnd.Next();
            area = c.PatchList[0].CalcArea(center);
            Console.Out.WriteLine("Area (rnd): {0}  ??  Girard: {1} ?? Est: {2}", area, girard_area, est_area);
            Assert.AreEqual(girard_area, area, 1e-5, "area w/ random center vs girard");
        }

        [Test]
        public void Semilune()
        {
            double ref_area = Math.PI / 3 - Math.Sqrt(3) / 4;
            Halfspace h = new Halfspace(180, 0, 60);
            Halfspace g = new Halfspace(181, 0, 60);
            Cartesian m = new Cartesian(181, 0);
            
            Cartesian p1, p2;
            h.Roots(g, out p1, out p2);
            Arc a = new Arc(h, p2, p1);
            Assert.AreEqual(ref_area, a.Area, 1e-6, "area");
            Assert.IsTrue(a.Middle.Same(m, Const.Tolerance), "middle Cartesian");

            p2 = p1;
            a = new Arc(h, p1, p2);
            Assert.AreEqual(0, a.Area , 1e-6, "zero area");
            Assert.IsTrue(a.Middle.Same(p1, Const.Tolerance), "zero middle Cartesian");
        }

        [Test]
        public void Cap()
        {
            Halfspace h = new Halfspace(180, 0, 6);
            Convex c = new Convex(h);

            double ref_area = 0.01 * Math.PI, area = 0;
            Console.Out.WriteLine(c);
            c.Simplify(Const.Tolerance);
            Assert.AreEqual(1, c.PatchList.Count, "one patch");
            Assert.AreEqual(1, c.PatchList[0].ArcList.Count, "one arc");
            foreach (Patch e in c.PatchList)
            {
                area += e.Area;
                Console.Out.WriteLine
                    ("Patch and total area {0} -> {1}", e.Area, area);
                Console.Out.WriteLine("Bounding Circle: {0}", 
                    e.Mec);
            }
            Assert.AreEqual(ref_area , area , 1e-7, "area"); 
        }

        [Test]
        public void Hole()
        {
            Halfspace h = new Halfspace(180, 0, 6);
            Convex c = new Convex(h.Inverse());

            double ref_area, area = 0;
            ref_area = Const.WholeSphereInSquareDegree - 0.01 * Math.PI;
            Console.Out.WriteLine(c);
            c.Simplify(Const.Tolerance);
            Assert.AreEqual(1, c.PatchList.Count, "one patch");
            Assert.AreEqual(1, c.PatchList[0].ArcList.Count, "one arc");
            foreach (Patch e in c.PatchList)
            {
                area += e.Area;
                Console.Out.WriteLine
                    ("Patch and total area {0} -> {1}", e.Area, area);
                Console.Out.WriteLine("Bounding Circle: {0}",
                    e.Mec);
            }
            Assert.AreEqual(ref_area, Const.WholeSphereInSquareDegree+area, 1e-7, "area");
        }

        [Test]
        public void Holes()
        {
            Convex c = new Convex(2);
            c.Add(new Halfspace(180, 0, 6).Inverse());
            c.Add(new Halfspace(180.19999, 0, 6).Inverse());

            double ref_area, area = 0;
            ref_area = Const.WholeSphereInSquareDegree - 2 * 0.01 * Math.PI;
            Console.Out.WriteLine(c);
            c.Simplify(Const.Tolerance);
            Assert.AreEqual(1, c.PatchList.Count, "one patch");
            Assert.AreEqual(2, c.PatchList[0].ArcList.Count, "two arcs");
            foreach (Patch e in c.PatchList)
            {
                area += e.Area;
                Console.Out.WriteLine
                    ("Patch and total area {0} -> {1}", e.Area, area);
                Console.Out.WriteLine("Bounding Circle: {0}",
                    e.Mec);
            }
            Assert.AreEqual(ref_area, Const.WholeSphereInSquareDegree+area, 1e-7, "area");
        }

        [Test]
        public void Ring()
        {
            double area = 0;
            double ref_area;
            Convex c = Samples.Ring(out ref_area);
            Console.Out.WriteLine(c);
            c.Simplify(Const.Tolerance);
            Assert.AreEqual(2, c.PatchList.Count, "two patch");
            Assert.AreEqual(1, c.PatchList[0].ArcList.Count, "one arc");
            Assert.AreEqual(1, c.PatchList[1].ArcList.Count, "one arc");
            foreach (Patch e in c.PatchList)
            {
                area += e.Area;
                Console.Out.WriteLine
                    ("Patch and total area {0} -> {1}", e.Area, area);
                Console.Out.WriteLine("Bounding Circle: {0}", e.Mec);
            }
            Assert.AreEqual(ref_area, area , 1e-7, "area"); 
        }

        [Test]
        public void OrangeSlice()
        {
            double ref_area; // = 4*Math.PI/8*Const.SquareRadian2SquareDegree;
            Convex c = Samples.OrangeSlice(out ref_area); // 45 deg angle
            double area = 0;
            c.Simplify(Const.Tolerance);
            Assert.AreEqual(2, c.HalfspaceList.Count, "2 halfspace");
            Assert.AreEqual(1, c.PatchList.Count, "only one patch");
            Assert.AreEqual(2, c.PatchList[0].ArcList.Count, "two arcs");
            foreach (Patch e in c.PatchList)
            {
                area += e.Area;
                Console.Out.WriteLine
                    ("Patch and total area {0} -> {1}", e.Area, area);
                Console.Out.WriteLine("Bounding Circle: {0}", e.Mec);
            }
            Assert.AreEqual(ref_area, area , 1e-5, "area");
        }

        [Test]
        public void OrangeSlices()
        {
            double ref_area; // = 4*Math.PI/8*Const.SquareRadian2SquareDegree;
            Convex c = Samples.OrangeSlices(out ref_area); // 45 deg angle
            double area = 0;
            c.Simplify(Const.Tolerance);
            Assert.AreEqual(1, c.PatchList.Count, "only one patch");
            Assert.AreEqual(2, c.PatchList[0].ArcList.Count, "two arcs");
            Assert.AreEqual(2, c.HalfspaceList.Count, "2 halfspace");
            foreach (Patch e in c.PatchList)
            {
                area += e.Area;
                Console.Out.WriteLine
                    ("Patch and total area {0} -> {1}", e.Area, area);
                 Console.Out.WriteLine("Bounding Circle: {0}", e.Mec);
            }
            Assert.AreEqual(ref_area, area , 1e-5, "area");
        }

        [Test]
        public void OrangeSlices2()
        {
            double ref_area; // = 4 * Math.PI / 4 * Const.SquareRadian2SquareDegree;
            Convex c = Samples.OrangeSlices2(out ref_area); // 90 deg angle
            double area = 0;
            c.Simplify(Const.Tolerance);
            Assert.AreEqual(1, c.PatchList.Count, "only one patch");
            Assert.AreEqual(2, c.PatchList[0].ArcList.Count, "two arcs");
            Assert.AreEqual(2, c.HalfspaceList.Count, "2 halfspace");
            foreach (Patch e in c.PatchList)
            {
                area += e.Area;
                Console.Out.WriteLine
                    ("Patch and total area {0} -> {1}", e.Area, area);
            }
            Assert.AreEqual(ref_area, area , 1e-5, "area");
        }

        [Test]
        public void OrangeSlicesEps()
        {
            double ref_area;
            Convex c = Samples.OrangeSlicesEps(out ref_area);
            c.SimpleSimplify(Const.Tolerance);
            Assert.AreEqual(3, c.HalfspaceList.Count, "orig 3 halfspace");
            double area = 0;
            c.Simplify(Const.Tolerance);
            Assert.AreEqual(2, c.HalfspaceList.Count, "2 halfspace");
            Assert.AreEqual(1, c.PatchList.Count, "only one patch");
            Assert.AreEqual(2, c.PatchList[0].ArcList.Count, "two arcs");
            foreach (Patch e in c.PatchList)
            {
                area += e.Area;
                Console.Out.WriteLine
                    ("Patch and total area {0} -> {1}", e.Area, area);
            }
            Assert.AreEqual(ref_area, area , 1e-5, "area");
        }

        [Test]
        public void Box()
        {
            double ref_area;
            Region r = Samples.Box(out ref_area);
            Assert.AreEqual(4, r.ConvexList.Count, "four convex");
            double area = 0;
            foreach (Convex c in r.ConvexList)
            {
                c.Simplify(Const.Tolerance);
                foreach (Patch e in c.PatchList)
                {
                    area += e.Area;
                    Console.Out.WriteLine
                        ("Patch and total area {0} -> {1}", e.Area, area);
                }
            }
            Assert.AreEqual(ref_area, area , 1e-5, "area");
        }

        [Test]
        public void Clover()
        {
            Convex c = Samples.Clover(true);
            Console.Out.WriteLine("Before:\n{0}", c);
            int Ns = c.HalfspaceList.Count;

            c.SimpleSimplify(Const.Tolerance);
            c.Simplify(Const.Tolerance);
            Console.Out.WriteLine("After:\n{0}", c);
            int Ne = c.HalfspaceList.Count;
            
            Assert.AreEqual(Ns-1, Ne, "one redudant halfspace");
            Assert.AreEqual(2, c.PatchList.Count, "two patches");
            Assert.IsTrue(c.PatchList[0].Area < 0 || c.PatchList[1].Area < 0, "one of them is a hole");

            double area = 0;
            foreach (Patch e in c.PatchList)
            {
                for (int i = 0; i < e.ArcList.Count; i++)
                {
                    // Console.Out.WriteLine("Arc: {0}", e[i]);
                }
                area += e.Area;
                Console.Out.WriteLine
                    ("Patch and total area {0} -> {1}", e.Area, area);
            }
        }

        [Test]
        public void Poly()
        {
            List<Cartesian> pts = new List<Cartesian>(3);
            pts.Add(new Cartesian(0, 0));
            pts.Add(new Cartesian(0.1, 0));
            pts.Add(new Cartesian(0.1, 0.1));
            pts.Add(new Cartesian(0, 0.1));
            Convex c = new Convex(pts, PointOrder.Safe);
            Assert.AreEqual(4, c.HalfspaceList.Count, "four halfspaces");

            Region r = new Region(c, true);
            c.Clear();
            double area = 0;
            foreach (Convex k in r.ConvexList)
            {
                k.Simplify(Const.Tolerance);
                foreach (Patch e in k.PatchList)
                {
                    area += e.Area;
                    Console.Out.WriteLine
                        ("Patch and total area {0} -> {1}", e.Area, area);
                }
            }
            Assert.AreEqual(0.01, area , 1e-5, "area of square");
        }

        [Test]
        public void ConvexCopy()
        {
            Convex c = new Convex(new Halfspace(180, 0, 60));
            Convex k = new Convex(c);
            c.HalfspaceList[0].Vector.Set(1, 0, 0, false);
            Assert.AreEqual(-1,k.HalfspaceList[0].Vector.X,Const.Tolerance);
        }

        [Test]
        public void GirardArea()
        {
            // big triangle
            Cartesian p = new Cartesian(0,0);
            Cartesian q = new Cartesian(90,0);
            Cartesian r = new Cartesian(0,90);

            double area = Cartesian.GirardArea(p,q,r);
            Console.Out.WriteLine("Area of triangle [{0}][{1}][{2}] is {3}", p, q, r, area);
            Assert.AreEqual(Math.PI/2*Const.SquareRadian2SquareDegree,area ,1e-6,"big");
            // small triange
            q = new Cartesian(0.1, 0);
            r = new Cartesian(0, 0.1);
            area = Cartesian.GirardArea(p, q, r);
            Console.Out.WriteLine("Area of triangle [{0}][{1}][{2}] is {3}", p, q, r, area);
            Assert.AreEqual(0.005, area , 1e-6,"small");
            // smaller triangle
            q = new Cartesian(0.00001, 0);
            r = new Cartesian(0, 0.00001);
            area = Cartesian.GirardArea(p, q, r);
            Console.Out.WriteLine("Area of triangle [{0}][{1}][{2}] is {3}", p, q, r, area);
            Assert.AreEqual(5e-11, area , 1e-11,"smaller");
            // empty triangle
            r = q;
            area = Cartesian.GirardArea(p, q, r);
            Console.Out.WriteLine("Area of triangle [{0}][{1}][{2}] is {3}", p, q, r, area);
            Assert.AreEqual(0, area , "empty");
            
        }

        [Test]
        public void GirardOrder()
        {
            RandomPoint r = new RandomPoint();
            Cartesian[] p = r.NextArray(3);
            double a1, a2, a3, b1, b2, b3;
            a1 = Cartesian.GirardArea(p[0], p[1], p[2]);
            a2 = Cartesian.GirardArea(p[1], p[2], p[0]);
            a3 = Cartesian.GirardArea(p[2], p[0], p[1]);
            b1 = Cartesian.GirardArea(p[0], p[2], p[1]);
            b2 = Cartesian.GirardArea(p[1], p[0], p[2]);
            b3 = Cartesian.GirardArea(p[2], p[1], p[0]);
            Assert.AreEqual(a1 , a2 , 1e-8);
            Assert.AreEqual(a1 , a3 , 1e-8);
            Assert.AreEqual(a1 , -b1 , 1e-8);
            Assert.AreEqual(a1 , -b2 , 1e-8);
            Assert.AreEqual(a1 , -b2 , 1e-8);
        }

        [Test]
        public void GirardCenter()
        {
            List<Cartesian> pts = new List<Cartesian>(3);
            pts.Add(new Cartesian(0, 0));
            pts.Add(new Cartesian(0.1, 0));
            pts.Add(new Cartesian(0.1, 0.1));
            pts.Add(new Cartesian(0, 0.1));
            Convex c = new Convex(pts, PointOrder.Safe);
            c.Simplify(Const.Tolerance);
            Assert.AreEqual(1,c.PatchList.Count,"1 patch");
            
            RandomPoint r = new RandomPoint();
            Cartesian m = r.Next();
            Console.Out.WriteLine("Random center Cartesian: {0}", m);

            double area = 0;
            foreach (Arc a in c.PatchList[0].ArcList)
                {
                    area += a.CalcArea(m,Const.Tolerance);
                    Console.Out.WriteLine
                        ("Triangle and total patch area {0} -> {1}", a.CalcArea(m,Const.Tolerance), area);
            }      
            Assert.AreEqual(0.01, area , 1e-5, "area of square");
        }

        [Test]
        public void RegionCopy()
        {
            Region r = new Region(180, 0, 60);
            Region q = new Region(r); // test if really clones it
            r.ConvexList[0].Clear();
            r.Clear();
            Assert.AreEqual(-1, q.ConvexList[0].HalfspaceList[0].Vector.X, 
                Const.Tolerance, "same if not removed");
        }

        [Test]
        public void CartesianCopy()
        {
            Cartesian p = new Cartesian(1, 0, 0, false);
            Cartesian q = p;
            p.Set(0, 0, 0, false);
            Assert.AreEqual(1, q.X);
        }

        [Test]
        public void CartesianRadec()
        {
            Cartesian p = new Cartesian(180, 0);
            Assert.AreEqual(-1, p.X, Const.Tolerance, "x coordinate");
            Assert.AreEqual(0, p.Y, Const.Tolerance, "y coordinate");
            Assert.AreEqual(0, p.Z, Const.Tolerance, "z coordinate");
        }

        [Test]
        public void CubePatches()
        {
            Convex c = Samples.Cube();
            c.Simplify(Const.Tolerance);
            Assert.AreEqual(8, c.PatchList.Count, "eight patches");
        }

        [Test]
        public void Growing()
        {
            Region r = new Region(0, 0, 60);
            r.Grow(60);
            r.Simplify();
            Console.Out.WriteLine("Area: {0}", r.Area);
            Assert.IsTrue(r.Area.Value > 12, "area ~ 4Pi");
        }

        [Test]
        public void ContainsPoint()
        {
            Region r = new Region(0, 0, 60);
            Cartesian c = new Cartesian(0, 0);
            bool isin = r.Covers(c, 0);
            Console.Out.WriteLine("Contains? {0}", isin);
            Assert.IsTrue(isin, "contains");
        }
    }

    [TestFixture]
    public class Problems
    {
        [Test]
        public void Intersect()
        {
            Convex c = Convex.Parse(@"CONVEX
  -0.76130275695572069 0.59849748704895045 -0.24943710679790604 0
  -0.35894347350399164 -0.78306424243025574 -0.50790744733269944 0
  0.76105698321837123 -0.59863910672181841 0.24984692953457566 0
  0.35844840531296612 0.78277010082156373 0.50870984852702694 0
  -0.76130586562456548 0.59849187858396891 -0.24944107567652724 0
  -0.35896747157119147 -0.7830476106002332 -0.50791612879199988 0
  0.76105387060565166 -0.59864472701805582 0.24984294436632318 0
  0.35842439981317614 0.78278674799164438 0.50870114683105461 0", false);
            
            bool b = false;
            b = c.SimpleSimplify(Const.Tolerance);
            Console.Out.WriteLine("Empty? {0}", !b);
            Console.Out.WriteLine("# of Halfspaces: {0}", c.HalfspaceList.Count);
            Console.Out.WriteLine("C:\n {0}", c);
            b = c.Simplify(false, Const.Tolerance, Const.MaximumIteration);
            Console.Out.WriteLine("Empty? {0}", !b);
            Console.Out.WriteLine("# of Halfspaces: {0}", c.HalfspaceList.Count);
            Console.Out.WriteLine("C:\n {0}", c);
            Console.Out.WriteLine("Area: {0}", c.Area);
        }
    }

    [TestFixture]
    public class RegionCSV
    {
        static string filename = @"C:\Tamas\Visual Studio 2005\Projects\Spherical\SphericalLibTest\regionarea.csv";

        [Test]
        public void AreaConvexes()
        {
            Console.Out.WriteLine("Filename: {0}", filename);
            StreamReader rdr = new StreamReader(filename);

            string line;
            int counter = 0, problem = 0;
            while ((line = rdr.ReadLine()) != null)
            {
                counter++;
                if (line.StartsWith("regionid")) continue;
                string[] tok = line.Split(new char[] { ',' }, 3);
                //Console.Error.WriteLine("{0} :: {1} :: {2}", tok[0], tok[1], tok[2]);

                Region r = Region.Parse(tok[2], true);
                r.Name = tok[0].Trim();
                double ref_area = double.Parse(tok[1]);

                // Console.Out.WriteLine("\t{0}   {1}", counter, r.Name);

                double area = 0;
                foreach (Convex c in r.ConvexList)
                {
                    bool cont = c.Simplify();
                    if (cont)
                    {
                        foreach (Patch e in c.PatchList)
                        {
                            area += e.Area;
                        }
                    }
                }
                //Assert.IsTrue(area >= 0, "not hole");

                double delta = Math.Abs(ref_area - area);
                double rel = delta / ref_area;
#if true
                if (rel > 1e-4)
                {
                    Console.Error.WriteLine("{0}  {1} \t {4:0.0000000}  {5:0.00000} \t {2:e}   {3:e}",
                        counter, r.Name, ref_area, area, delta, rel);
                    problem++;
                }
#endif
            }
            rdr.Close();
            Console.Error.WriteLine("# of regions w/ different area measurements: {0}",
                problem);
        }

        [Test]
        public void AreaRegions()
        {
            Console.Out.WriteLine("Filename: {0}", filename);
            StreamReader rdr = new StreamReader(filename);

            string line;
            int counter = 0, problem = 0;
            while ((line = rdr.ReadLine()) != null)
            {
                counter++;
                if (line.StartsWith("regionid")) continue;
                string[] tok = line.Split(new char[] { ',' }, 3);
                //Console.Error.WriteLine("{0} :: {1} :: {2}", tok[0], tok[1], tok[2]);

                Region r = Region.Parse(tok[2], true);
                r.Name = tok[0].Trim();
                double ref_area = double.Parse(tok[1]);

                //if (r.Name == "178685")
                //{
                //    Console.Out.WriteLine(r.Name);
                //}
                // Console.Out.WriteLine("\t{0}   {1}", counter, r.Name);

                r.Simplify(true, Const.Tolerance, Const.MaximumIteration, true);
                double area = r.Area.Value;

                double delta = Math.Abs(ref_area - area);
                double rel = delta / ref_area;
#if true
                if (rel > 1e-4)
                {
                    Console.Error.WriteLine("{0}  {1} \t {4:0.0000000}  {5:0.00000} \t {2:e}   {3:e}",
                        counter, r.Name, ref_area, area, delta, rel);
                    problem++;
                }
#endif
            }
            rdr.Close();
            Console.Error.WriteLine("# of regions w/ different area measurements: {0}",
                problem);
        }
    }

    [TestFixture]
    public class Dr4Footprint
    {
        [Test]
        public void Footprint()
        {
            Region r = Samples.SdssDr4();
            r.SimpleSimplify(Const.Tolerance);

            RandomPoint rnd = new RandomPoint();
            double parea, garea, tarea;
            foreach (Convex c in r.ConvexList)
            {
                c.Simplify(Const.Tolerance);
                Assert.AreEqual(1, c.PatchList.Count, "1 patch");
                Patch p = c.PatchList[0];
                Cartesian center = new Cartesian(183.70622459853212, 14.996375195638675);
                center = rnd.Next();
                Console.Error.WriteLine("Center/pivot Cartesian: {0}", center);
                parea = garea = 0;
                foreach (Arc a in p.ArcList)
                {
                    tarea = a.CalcArea(center,Const.Tolerance);
                    //Console.Out.WriteLine("Arc: \n{0}Center: {1} {2}\n\t{3}\nTarea: {4}\n",
                    //    a, center.RA, center.Dec, center, tarea);
                    parea += tarea;
                    garea += Cartesian.GirardArea(center, a.Point1, a.Point2);
                }
                Console.Error.WriteLine("Total patch area: {0}", parea);
                Console.Error.WriteLine("Total girard area: {0}\n", garea);
            }
        }

        [Test]
        public void Patch0()
        {
            RandomPoint rnd = new RandomPoint();
            Region r = Samples.SdssDr4();
            double area = 0, garea = 0, tarea;
            r.SimpleSimplify(Const.Tolerance);
            foreach (Convex c in r.ConvexList)
            {
                c.Simplify(Const.Tolerance);
                Patch p = c.PatchList[0];
                Cartesian center = rnd.Next();
                //center = new Cartesian(183.70622459853212, 14.996375195638675);
                Console.Error.WriteLine("Random center/pivot Cartesian: {0}", center);
                foreach (Arc a in p.ArcList)
                {
                    tarea = a.CalcArea(center,Const.Tolerance);
                    Console.Out.WriteLine("Arc: \n{0}Center: {1} {2}\n\t{3}\nTarea: {4}\n", 
                        a, center.RA,center.Dec, center, tarea); 
                    area += tarea;
                    garea += Cartesian.GirardArea(center, a.Point1, a.Point2);
                }
                break;
            }
            Console.Error.WriteLine("Total patch area: {0}", area);
            Console.Error.WriteLine("Total girard area: {0}", garea);
        }

        [Test]
        public void Patch0_Reduced()
        {
            RandomPoint rnd = new RandomPoint();
            Region r = Samples.SdssDr4();
            double area = 0, garea = 0, tarea;
            r.SimpleSimplify(Const.Tolerance);
            foreach (Convex c in r.ConvexList)
            {
                c.Simplify(Const.Tolerance);
                Patch p = c.PatchList[0];

                Arc[] arcs = new Arc[4];
                arcs[0] = p.ArcList[0];
                arcs[1] = p.ArcList[1];
                arcs[2] = p.ArcList[2];
                arcs[3] = new Arc(p.ArcList[2].Point2,p.ArcList[0].Point1);
                Patch n = new Patch(arcs,Const.Tolerance);
                Cartesian center = rnd.Next();
                //center = new Cartesian(183.70622459853212, 14.996375195638675);
                Console.Error.WriteLine("Random center/pivot Cartesian: {0}", center);
                foreach (Arc a in n.ArcList)
                {
                    tarea = a.CalcArea(center,Const.Tolerance);
                    Console.Out.WriteLine("Arc: \n{0}Center: {1} {2}\n\t{3}\nTarea: {4}\n",
                        a, center.RA, center.Dec, center, tarea);
                    area += tarea;
                    garea += Cartesian.GirardArea(center, a.Point1, a.Point2);
                }
                break;
            }
            Console.Error.WriteLine("Total patch area: {0}", area);
            Console.Error.WriteLine("Total girard area: {0}", garea);
        }

        [Test]
        public void Patch0_Reduced2()
        {
            RandomPoint rnd = new RandomPoint();
            Region r = Samples.SdssDr4();
            double area = 0, garea = 0, tarea;
            r.SimpleSimplify(Const.Tolerance);
            foreach (Convex c in r.ConvexList)
            {
                c.Simplify(Const.Tolerance);
                Patch p = c.PatchList[0];

                Arc[] arcs = new Arc[3];
                arcs[0] = p.ArcList[1];
                arcs[1] = p.ArcList[2];
                arcs[2] = new Arc(p.ArcList[2].Point2, p.ArcList[1].Point1);
                Patch n = new Patch(arcs,Const.Tolerance);
                Cartesian center = rnd.Next();
                //center = new Cartesian(183.70622459853212, 14.996375195638675);
                Console.Error.WriteLine("Random center/pivot Cartesian: {0}", center);
                foreach (Arc a in n.ArcList)
                {
                    tarea = a.CalcArea(center,Const.Tolerance);
                    Console.Out.WriteLine("Arc: \n{0}Center: {1} {2}\n\t{3}\nTarea: {4}\n",
                        a, center.RA, center.Dec, center, tarea);
                    area += tarea;
                    garea += Cartesian.GirardArea(center, a.Point1, a.Point2);
                }
                break;
            }
            Console.Error.WriteLine("Total patch area: {0}", area);
            Console.Error.WriteLine("Total girard area: {0}", garea);
        }

        [Test]
        public void Patch0_Rest()
        {
            RandomPoint rnd = new RandomPoint();
            Region r = Samples.SdssDr4();
            double area = 0, garea = 0, tarea;
            r.SimpleSimplify(Const.Tolerance);
            foreach (Convex c in r.ConvexList)
            {
                c.Simplify(Const.Tolerance);
                Patch p = c.PatchList[0];

                Arc[] arcs = new Arc[p.ArcList.Count-1];
                arcs[0] = p.ArcList[0];
                arcs[1] = new Arc(p.ArcList[0].Point2, p.ArcList[3].Point1);
                for (int i = 3; i < p.ArcList.Count; i++)
                    arcs[i - 1] = p.ArcList[i];
                Patch n = new Patch(arcs,Const.Tolerance);
                Cartesian center = rnd.Next();
                //center = new Cartesian(183.70622459853212, 14.996375195638675);
                Console.Error.WriteLine("Random center/pivot Cartesian: {0}", center);
                foreach (Arc a in n.ArcList)
                {
                    tarea = a.CalcArea(center,Const.Tolerance);
                    Console.Out.WriteLine("Arc: \n{0}Center: {1} {2}\n\t{3}\nTarea: {4}\n",
                        a, center.RA, center.Dec, center, tarea);
                    area += tarea;
                    garea += Cartesian.GirardArea(center, a.Point1, a.Point2);
                }
                break;
            }
            Console.Error.WriteLine("Total patch area: {0}", area);
            Console.Error.WriteLine("Total girard area: {0}", garea);
        }

        [Test]
        public void Patch0_Triangle1_Girard()
        {
            Cartesian p, q, r;
            double area;
            // center
            r = new Cartesian(183.706224598532, 14.9963751956387);

            // arc '1' (2nd)
            p = new Cartesian(246.58837, 7.89564000000001);
            q = new Cartesian(133.418029967136, 10.2678933742546);
            Console.Out.WriteLine("Cartesians:\n\t{0}\n\t{1}\n\t{2}", p, q, r);
            Console.Out.WriteLine("\t1: {0} {1}", p.RA, p.Dec);
            Console.Out.WriteLine("\t2: {0} {1}", q.RA, r.Dec);
            Console.Out.WriteLine("\t3: {0} {1}", r.RA, q.Dec);
            area = Cartesian.GirardArea(p, q, r);
            Console.Out.WriteLine("T1: {0}", area);
            Assert.AreEqual(75.76979765, area , 1e-5);

            // arc '1' w/ xyz
            //p = Cartesian.Parse("-0.393567409974194 -0.908974369568987 0.137369171479144");
            q = Cartesian.Parse("-0.676308803094126 0.71472584741758 0.178250851023529", true);
            Console.Out.WriteLine("Cartesians:\n\t{0}\n\t{1}\n\t{2}", p, q, r);
            Console.Out.WriteLine("\t1: {0} {1}", p.RA, p.Dec);
            Console.Out.WriteLine("\t2: {0} {1}", q.RA, r.Dec);
            Console.Out.WriteLine("\t3: {0} {1}", r.RA, q.Dec);
            area = Cartesian.GirardArea(p, q, r);
            Console.Out.WriteLine("T1: {0}", area);
            Assert.AreEqual(75.76979765, area , 1e-5);
            
            // other
            //q = Cartesian.Parse("-0.676428168546963 0.714851993305448 0.17825085102353", true);
            q = new Cartesian(q.RA, q.Dec);
            Console.Out.WriteLine("Cartesians:\n\t{0}\n\t{1}\n\t{2}", p, q, r);
            Console.Out.WriteLine("\t1: {0} {1}", p.RA, p.Dec);
            Console.Out.WriteLine("\t2: {0} {1}", q.RA, r.Dec);
            Console.Out.WriteLine("\t3: {0} {1}", r.RA, q.Dec);
            area = Cartesian.GirardArea(p, q, r);
            Console.Out.WriteLine("T1: {0}", area);
            Assert.AreEqual(75.76979765, area , 1e-5);

        }

        [Test]
        public void Region_179922()
        {
            string repr = "REGION CONVEX 0.087155742747658 -0.996194698091746 0.0 0.429175990306221 -0.871116348225734 0.48573730257896 0.072218977680232 0.0 -0.09980681579004 -0.008731964924249 -0.994968518250912 0.0 0.065154250577678 0.005700258306074 0.997858923238603 0.0 -0.087155742747658 0.996194698091746 0.0 -0.828619840776703";
            Region r = Region.Parse(repr,true);
            Console.Out.WriteLine(r);

            double ref_area = 45.771099356918953, area = 0;
            foreach (Convex c in r.ConvexList)
            {
                c.Simplify(Const.Tolerance);
#if true
                foreach (Patch e in c.PatchList)
                {
                    area += e.Area;
                    Console.Out.WriteLine
                        ("Patch and total area {0} -> {1}", e.Area, area);
                }
#else
                f.Test();
#endif
            }
            Assert.AreEqual(ref_area , area , 1e-5, "area");
        }

        [Test]
        public void Region_179920()
        {
            string repr = "REGION CONVEX -0.108452603525345 -0.009488373338339 -0.99405633822232 0.0 0.065154250577678 0.005700258306074 0.997858923238603 0.0 0.871116348225734 -0.48573730257896 -0.072218977680232 0.0 -0.087155742747658 0.996194698091746 0.0 -0.828619840776703";
            Region r = Region.Parse(repr,true);
            Console.Out.WriteLine(r);

            double ref_area = 2.9510634609093336E-3, area = 0;
            foreach (Convex c in r.ConvexList)
            {
                c.Simplify(Const.Tolerance);
                foreach (Patch e in c.PatchList)
                {
                    area += e.Area;
                    Console.Out.WriteLine
                        ("Patch and total area {0} -> {1}", e.Area, area);
                }
            }
            Assert.AreEqual(ref_area , area , 1e-5, "area");
        }
    }
}
