using System;using ESRI.ArcGIS.Client.Geometry;namespace GISProject.Extensions{    /// 
    /// Extension methods for geodesic calculations.    /// 
    public static class Geodesic    {        private const double EarthRadius = 6378.137; //kilometers. Change to miles to return all values in miles instead        /// 
        /// Gets the distance between two points in Kilometers.        /// 
        /// The start point.        /// The end point.        /// 
             public static double GetSphericalDistance(this MapPoint start, MapPoint end)        {            double lon1 = start.X / 180 * Math.PI;            double lon2 = end.X / 180 * Math.PI;            double lat1 = start.Y / 180 * Math.PI;            double lat2 = end.Y / 180 * Math.PI;            return 2 * Math.Asin(Math.Sqrt(Math.Pow((Math.Sin((lat1 - lat2) / 2)), 2) +             Math.Cos(lat1) * Math.Cos(lat2) * Math.Pow(Math.Sin((lon1 - lon2) / 2), 2))) * EarthRadius;        }        /// 
        /// Returns a polygon with a constant distance from the center point measured on the sphere.        /// 
        /// The center.        /// Radius in kilometers.        /// 
             public static Polygon GetRadiusAsPolygon(this MapPoint center, double distKM)        {            Polyline line = GetRadius(center, distKM);            Polygon poly = new Polygon();            if (line.Paths.Count > 1)            {                PointCollection ring = line.Paths[0];                MapPoint last = ring[ring.Count - 1];                for (int i = 1; i < line.Paths.Count; i++)                {                    PointCollection pnts = line.Paths[i];                    ring.Add(new MapPoint(180 * Math.Sign(last.X), 90 * Math.Sign(center.Y)));                    last = pnts[0];                    ring.Add(new MapPoint(180 * Math.Sign(last.X), 90 * Math.Sign(center.Y)));                    foreach (MapPoint p in pnts)                        ring.Add(p);                    last = pnts[pnts.Count - 1];                }                poly.Rings.Add(ring);                //pnts.Add(first);            }            else            {                poly.Rings.Add(line.Paths[0]);            }            if (distKM > EarthRadius * Math.PI / 2 && line.Paths.Count != 2)            {                PointCollection pnts = new PointCollection();                pnts.Add(new MapPoint(-180, -90));                pnts.Add(new MapPoint(180, -90));                pnts.Add(new MapPoint(180, 90));                pnts.Add(new MapPoint(-180, 90));                pnts.Add(new MapPoint(-180, -90));                poly.Rings.Add(pnts); //Exterior            }            return poly;        }        /// 
        /// Returns a polyline with a constant distance from the center point measured on the sphere.        /// 
        /// The center.        /// Radius in kilometers.        // 
             public static Polyline GetRadius(this MapPoint center, double distKM)        {            Polyline line = new Polyline();            PointCollection pnts = new PointCollection();            line.Paths.Add(pnts);            for (int i = 0; i < 360; i++)            {                //double angle = i / 180.0 * Math.PI;                MapPoint p = GetPointFromHeading(center, distKM, i);                if (pnts.Count > 0)                {                    MapPoint lastPoint = pnts[pnts.Count - 1];                    int sign = Math.Sign(p.X);                    if (Math.Abs(p.X - lastPoint.X) > 180)                    {   //We crossed the date line                        double lat = LatitudeAtLongitude(lastPoint, p, sign * -180);                        pnts.Add(new MapPoint(sign * -180, lat));                        pnts = new PointCollection();                        line.Paths.Add(pnts);                        pnts.Add(new MapPoint(sign * 180, lat));                    }                }                pnts.Add(p);            }            pnts.Add(line.Paths[0][0]);            return line;        }        /// 
        /// Gets the shortest path line between two points. THe line will be following the great        /// circle described by the two points.        /// 
        /// The start point.        /// The end point.        /// 
             public static Polyline GetGeodesicLine(this MapPoint start, MapPoint end)        {            Polyline line = new Polyline();            if (Math.Abs(end.X - start.X) <= 180) // Doesn't cross dateline            {                PointCollection pnts = GetGeodesicPoints(start, end);                line.Paths.Add(pnts);            }            else            {                double lon1 = start.X / 180 * Math.PI;                double lon2 = end.X / 180 * Math.PI;                double lat1 = start.Y / 180 * Math.PI;                double lat2 = end.Y / 180 * Math.PI;                double latA = LatitudeAtLongitude(lat1, lon1, lat2, lon2, Math.PI) / Math.PI * 180;                //double latB = LatitudeAtLongitude(lat1, lon1, lat2, lon2, -180) / Math.PI * 180;                line.Paths.Add(GetGeodesicPoints(start, new MapPoint(start.X < 0 ? -180 : 180, latA)));                line.Paths.Add(GetGeodesicPoints(new MapPoint(start.X < 0 ? 180 : -180, latA), end));            }            return line;        }        private static PointCollection GetGeodesicPoints(MapPoint start, MapPoint end)        {            double lon1 = start.X / 180 * Math.PI;            double lon2 = end.X / 180 * Math.PI;            double lat1 = start.Y / 180 * Math.PI;            double lat2 = end.Y / 180 * Math.PI;            double dX = end.X - start.X;            int points = (int)Math.Floor(Math.Abs(dX));            dX = lon2 - lon1;            PointCollection pnts = new PointCollection();            pnts.Add(start);            for (int i = 1; i < points; i++)            {                double lon = lon1 + dX / points * i;                double lat = LatitudeAtLongitude(lat1, lon1, lat2, lon2, lon);                pnts.Add(new MapPoint(lon / Math.PI * 180, lat / Math.PI * 180));            }            pnts.Add(end);            return pnts;        }        /// 
        /// Gets the latitude at a specific longitude for a great circle defined by p1 and p2.        /// 
        /// The p1.        /// The p2.        /// The longitude in degrees.        /// 
             private static double LatitudeAtLongitude(MapPoint p1, MapPoint p2, double lon)        {            double lon1 = p1.X / 180 * Math.PI;            double lon2 = p2.X / 180 * Math.PI;            double lat1 = p1.Y / 180 * Math.PI;            double lat2 = p2.Y / 180 * Math.PI;            lon = lon / 180 * Math.PI;            return LatitudeAtLongitude(lat1, lon1, lat2, lon2, lon) / Math.PI * 180;        }        /// 
        /// Gets the latitude at a specific longitude for a great circle defined by lat1,lon1 and lat2,lon2.        /// 
        /// The start latitude in radians.        /// The start longitude in radians.        /// The end latitude in radians.        /// The end longitude in radians.        /// The longitude in radians for where the latitude is.        /// 
             private static double LatitudeAtLongitude(double lat1, double lon1, double lat2, double lon2, double lon)        {            return Math.Atan((Math.Sin(lat1) * Math.Cos(lat2) * Math.Sin(lon - lon2)     - Math.Sin(lat2) * Math.Cos(lat1) * Math.Sin(lon - lon1)) / (Math.Cos(lat1) * Math.Cos(lat2) * Math.Sin(lon1 - lon2)));        }        /// 
        /// Gets the true bearing at a distance from the start point towards the new point.        /// 
        /// The start point.        /// The point to get the bearing towards.        /// The distance in kilometers travelled between start and end.        /// 
             public static double GetTrueBearing(MapPoint start, MapPoint end, double distanceKM)        {            double d = distanceKM / EarthRadius; //Angular distance in radians            double lon1 = start.X / 180 * Math.PI;            double lat1 = start.Y / 180 * Math.PI;            double lon2 = end.X / 180 * Math.PI;            double lat2 = end.Y / 180 * Math.PI;            double tc1;            if (Math.Sin(lon2 - lon1) < 0)                tc1 = Math.Acos((Math.Sin(lat2) - Math.Sin(lat1) * Math.Cos(d)) / (Math.Sin(d) * Math.Cos(lat1)));            else                tc1 = 2 * Math.PI - Math.Acos((Math.Sin(lat2) - Math.Sin(lat1) * Math.Cos(d)) / (Math.Sin(d) * Math.Cos(lat1)));            return tc1 / Math.PI * 180;        }        /// 
        /// Gets the point based on a start point, a heading and a distance.        /// 
        /// The start.        /// The distance KM.        /// The heading.        /// 
             public static MapPoint GetPointFromHeading(MapPoint start, double distanceKM, double heading)        {            double brng = heading / 180 * Math.PI;            double lon1 = start.X / 180 * Math.PI;            double lat1 = start.Y / 180 * Math.PI;            double dR = distanceKM / 6378.137; //Angular distance in radians            double lat2 = Math.Asin(Math.Sin(lat1) * Math.Cos(dR) + Math.Cos(lat1) * Math.Sin(dR) * Math.Cos(brng));            double lon2 = lon1 + Math.Atan2(Math.Sin(brng) * Math.Sin(dR) * Math.Cos(lat1), Math.Cos(dR) - Math.Sin(lat1) * Math.Sin(lat2));            double lon = lon2 / Math.PI * 180;            double lat = lat2 / Math.PI * 180;            while (lon < -180) lon += 360;            while (lat < -90) lat += 180;            while (lon > 180) lon -= 360;            while (lat > 90) lat -= 180;            return new MapPoint(lon, lat);        }    }}