search

airbreather / Gavaghan.Geodesy

A C# implementation of Vincenty's formulae
Other
9 stars 8 forks source link

Remove GlobalCoordinates and Generalize #6

Open juliusfriedman opened 3 years ago

juliusfriedman commented 3 years ago

Define derived types of Angle and a Location class which contains them.

See also: https://github.com/airbreather/Gavaghan.Geodesy/issues/2 https://github.com/airbreather/Gavaghan.Geodesy/issues/3

/// <summary>Represents a latitude ("y" axis) co-ordinate.</summary>
    public sealed class Latitude : Angle
    {
        const double PIOver2 = Math.PI / 2;

        /// <summary>Initializes a new instance of the Latitude class.</summary>
        /// <param name="angle">The angle of the latitude.</param>
        /// <exception cref="ArgumentNullException">angle is null.</exception>
        /// <exception cref="ArgumentOutOfRangeException">
        /// angle is greater than 90 degrees or less than -90 degrees.
        /// </exception>
        public Latitude(Angle angle)
            : base((angle ?? new Angle(0)).Radians) // Prevent null reference access
        {
            if (angle == null)
            {
                throw new ArgumentNullException("angle");
            }
            ValidateRange("angle", angle.Radians, -PIOver2, PIOver2);
        }

        private Latitude(double radians)
            : base(radians)
        {
        }

        /// <summary>
        /// Gets a value indicating whether this instance represents a north or
        /// south latitude.
        /// </summary>
        public CompassDirection Direction => Radians < 0 ? CompassDirection.S : CompassDirection.N;

        /// <summary>Creates a new Latitude from an angle in degrees.</summary>
        /// <param name="degrees">The angle of the latitude in degrees.</param>
        /// <returns>A new Latitude representing the specified value.</returns>
        /// <exception cref="ArgumentOutOfRangeException">
        /// degrees is greater than 90 or less than -90.
        /// </exception>
        public static new Latitude FromDegrees(double degrees)
        {
            ValidateRange("degrees", degrees, -90, 90);
            return new Latitude(Angle.FromDegrees(degrees).Radians);
        }

        /// <summary>
        /// Creates a new Latitude from an angle in degrees and minutes.
        /// </summary>
        /// <param name="degrees">The amount of degrees.</param>
        /// <param name="minutes">The amount of minutes.</param>
        /// <returns>A new Latitude representing the specified value.</returns>
        /// <exception cref="ArgumentOutOfRangeException">
        /// The specified angle (degrees + minutes) is greater than 90 or less
        /// than -90.
        /// </exception>
        public static new Latitude FromDegrees(double degrees, double minutes)
        {
            var angle = Angle.FromDegrees(degrees, minutes);
            ValidateRange("angle", angle.TotalDegrees, -90, 90);

            return new Latitude(angle.Radians);
        }

        /// <summary>
        /// Creates a new Latitude from an angle in degrees, minutes and seconds.
        /// </summary>
        /// <param name="degrees">The amount of degrees.</param>
        /// <param name="minutes">The amount of minutes.</param>
        /// <param name="seconds">The amount of seconds.</param>
        /// <returns>A new Latitude representing the specified value.</returns>
        /// <exception cref="ArgumentOutOfRangeException">
        /// The specified angle (degrees + minutes + seconds) is greater than
        /// 90 or less than -90.
        /// </exception>
        public static new Latitude FromDegrees(double degrees, double minutes, double seconds)
        {
            var angle = Angle.FromDegrees(degrees, minutes, seconds);
            ValidateRange("angle", angle.TotalDegrees, -90, 90);

            return new Latitude(angle.Radians);
        }

        /// <summary>Creates a new Latitude from an amount in radians.</summary>
        /// <param name="radians">The angle of the latitude in radians.</param>
        /// <returns>A new Latitude representing the specified value.</returns>
        /// <exception cref="ArgumentOutOfRangeException">
        /// radians is greater than PI/2 or less than -PI/2.
        /// </exception>
        public static new Latitude FromRadians(double radians)
        {
            ValidateRange("radians", radians, -PIOver2, PIOver2);
            return new Latitude(radians);
        }

        /// <summary>
        /// Formats the value of the current instance using the specified format.
        /// </summary>
        /// <param name="format">
        /// The format to use (see remarks) or null to use the default format.
        /// </param>
        /// <param name="formatProvider">
        /// The provider to use to format the value or null to use the format
        /// information from the current locale setting of the operating system.
        /// </param>
        /// <returns>
        /// The value of the current instance in the specified format.
        /// </returns>
        /// <exception cref="ArgumentException">format is unknown.</exception>
        /// <remarks>
        /// Valid format strings are those for
        /// <see cref="Angle.ToString(string, IFormatProvider)"/> plus "ISO"
        /// (without any precision specifier), which returns the angle in
        /// ISO 6709 compatible format.
        /// </remarks>
        public override string ToString(string format, IFormatProvider formatProvider)
        {
            if (format == "ISO")
            {
                char sign = this.Radians < 0 ? '-' : '+';
                return string.Format(
                    CultureInfo.InvariantCulture, // ISO defines the punctuation
                    "{0}{1:00}{2:00}{3:00.####}",
                    sign,
                    Math.Abs(this.Degrees),
                    Math.Abs(this.Minutes),
                    Math.Abs(this.Seconds));
            }

            string formatted = base.ToString(format, formatProvider);

            // We're going to remove the negative sign, but find out what a
            // negative sign is in the current format provider
            var numberFormat = NumberFormatInfo.GetInstance(formatProvider);
            string negativeSign = numberFormat.NegativeSign;
            if (formatted.StartsWith(negativeSign, StringComparison.Ordinal))
            {
                formatted = formatted.Substring(negativeSign.Length);
            }

            return formatted + " " + this.Direction;
        }
    }

 /// <summary>Represents a longitude ("x" axis) co-ordinate.</summary>
    public sealed class Longitude : Angle
    {
        /// <summary>Initializes a new instance of the Longitude class.</summary>
        /// <param name="angle">The angle of the longitude.</param>
        /// <exception cref="ArgumentNullException">angle is null.</exception>
        /// <exception cref="ArgumentOutOfRangeException">
        /// angle is greater than 180 degrees or less than -180 degrees.
        /// </exception>
        public Longitude(Angle angle)
            : this((angle ?? new Angle(0)).Radians) // Prevent null reference access, we'll validate later
        {
            if (angle == null)
            {
                throw new ArgumentNullException("angle");
            }
            ValidateRange("angle", angle.Radians, -Math.PI, Math.PI);
        }

        private Longitude(double radians)
            : base(radians == Math.PI ? -Math.PI : radians)
        {
            // The above test is a special case. According to the ISO 6709, the
            // 180th meridian (180 degrees == Pi radians) is always -180 degrees.
            // Instead of throwing an exception we'll just change the value.
        }

        /// <summary>
        /// Gets a value indicating whether this instance represents am east or
        /// west longitude.
        /// </summary>
        /// <summary>
        public CompassDirection Direction => Radians < 0 ? CompassDirection.W : CompassDirection.E;

        /// <summary>Creates a new Longitude from an angle in degrees.</summary>
        /// <param name="degrees">The angle of the longitude in degrees.</param>
        /// <returns>A new Longitude representing the specified value.</returns>
        /// <exception cref="ArgumentOutOfRangeException">
        /// degrees is greater than 180 or less than -180.
        /// </exception>
        public static new Longitude FromDegrees(double degrees)
        {
            ValidateRange("degrees", degrees, -180, 180);
            return new Longitude(Angle.FromDegrees(degrees).Radians);
        }

        /// <summary>
        /// Creates a new Longitude from an angle in degrees and minutes.
        /// </summary>
        /// <param name="degrees">The amount of degrees.</param>
        /// <param name="minutes">The amount of minutes.</param>
        /// <returns>A new Longitude representing the specified value.</returns>
        /// <exception cref="ArgumentOutOfRangeException">
        /// The specified angle (degrees + minutes) is greater than 180 or less
        /// than -180.
        /// </exception>
        public static new Longitude FromDegrees(double degrees, double minutes)
        {
            var angle = Angle.FromDegrees(degrees, minutes);
            ValidateRange("angle", angle.TotalDegrees, -180, 180);

            return new Longitude(angle.Radians);
        }

        /// <summary>
        /// Creates a new Longitude from an angle in degrees, minutes and seconds.
        /// </summary>
        /// <param name="degrees">The amount of degrees.</param>
        /// <param name="minutes">The amount of minutes.</param>
        /// <param name="seconds">The amount of seconds.</param>
        /// <returns>A new Longitude representing the specified value.</returns>
        /// <exception cref="ArgumentOutOfRangeException">
        /// The specified angle (degrees + minutes + seconds) is greater than
        /// 180 or less than -180.
        /// </exception>
        public static new Longitude FromDegrees(double degrees, double minutes, double seconds)
        {
            var angle = Angle.FromDegrees(degrees, minutes, seconds);
            ValidateRange("angle", angle.TotalDegrees, -180, 180);

            return new Longitude(angle.Radians);
        }

        /// <summary>Creates a new Longitude from an amount in radians.</summary>
        /// <param name="radians">The angle of the longitude in radians.</param>
        /// <returns>A new Longitude representing the specified value.</returns>
        /// <exception cref="ArgumentOutOfRangeException">
        /// radians is greater than PI or less than -PI.
        /// </exception>
        public static new Longitude FromRadians(double radians)
        {
            ValidateRange("radians", radians, -Math.PI, Math.PI);
            return new Longitude(radians);
        }

        /// <summary>
        /// Formats the value of the current instance using the specified format.
        /// </summary>
        /// <param name="format">
        /// The format to use (see remarks) or null to use the default format.
        /// </param>
        /// <param name="formatProvider">
        /// The provider to use to format the value or null to use the format
        /// information from the current locale setting of the operating system.
        /// </param>
        /// <returns>
        /// The value of the current instance in the specified format.
        /// </returns>
        /// <exception cref="ArgumentException">format is unknown.</exception>
        /// <remarks>
        /// Valid format strings are those for
        /// <see cref="Angle.ToString(string, IFormatProvider)"/> plus "ISO"
        /// (without any precision specifier), which returns the angle in
        /// ISO 6709 compatible format.
        /// </remarks>
        public override string ToString(string format, IFormatProvider formatProvider)
        {
            if (format == "ISO")
            {
                char sign = this.Radians < 0 ? '-' : '+';
                return string.Format(
                    CultureInfo.InvariantCulture, // ISO defines the punctuation
                    "{0}{1:000}{2:00}{3:00.####}",
                    sign,
                    Math.Abs(this.Degrees),
                    Math.Abs(this.Minutes),
                    Math.Abs(this.Seconds));
            }

            string formatted = base.ToString(format, formatProvider);

            // We're going to remove the negative sign, but find out what a
            // negative sign is in the current format provider
            var numberFormat = NumberFormatInfo.GetInstance(formatProvider);
            string negativeSign = numberFormat.NegativeSign;
            if (formatted.StartsWith(negativeSign, StringComparison.Ordinal))
            {
                formatted = formatted.Substring(negativeSign.Length);
            }

            return formatted + " " + this.Direction;
        }
    }

/// <summary>
    /// Defines the Geographic Position type a.k.a. <see href="http://geojson.org/geojson-spec.html#positions">Geographic Coordinate Reference System</see>.
    /// </summary>
    /// <summary>Represents a Latitude/Longitude/Altitude coordinate. (on Earth)</summary>
    public sealed partial class Location : Position, IEquatable<Location>, IFormattable
    {
        /// <summary>
        /// 
        /// </summary>
        /// <param name="p"></param>
        /// <returns></returns>
        public static Location FromPoint(Point p) => new Location(new Latitude(new Angle(p.Y)), new Longitude(new Angle(p.X)), p.Z);

        // Equatorial = 6378137, polar = 6356752.
        internal const int EarthRadius = 6366710; // The mean radius of the Earth in meters.

        private Latitude latitude;
        private Longitude longitude;
        private double? altitude;
        //TODO, could be Body{Radius, Flattening, InverseFlattening, UniversalLocation}
        private double radii = EarthRadius;

        /// <summary>Initializes a new instance of the Location class.</summary>
        /// <param name="latitude">The latitude of the coordinate.</param>
        /// <param name="longitude">The longitude of the coordinate.</param>
        /// <exception cref="ArgumentNullException">
        /// latitude/longitude is null.
        /// </exception>
        public Location(Latitude latitude, Longitude longitude)
        {
            if (latitude == null)
            {
                throw new ArgumentNullException("latitude");
            }
            if (longitude == null)
            {
                throw new ArgumentNullException("longitude");
            }

            this.latitude = latitude;
            this.longitude = longitude;
        }

        /// <summary>Initializes a new instance of the Location class.</summary>
        /// <param name="latitude">The latitude of the coordinate.</param>
        /// <param name="longitude">The longitude of the coordinate.</param>
        /// <param name="altitude">
        /// The altitude, specifed in meters, of the coordinate.
        /// </param>
        /// <exception cref="ArgumentNullException">
        /// latitude/longitude is null.
        /// </exception>
        public Location(Latitude latitude, Longitude longitude, double altitude)
            : this(latitude, longitude) => this.altitude = altitude;

        // XmlSerializer requires a parameterless constructor
        private Location()
        {
        }

        /// <summary>
        /// The Radii of the body in which this instance is plotted
        /// </summary>
        public double Radii => radii;

        /// <summary>
        /// Gets the altitude of the coordinate, or null if the coordinate
        /// does not contain altitude information.
        /// </summary>
        public double? Altitude => altitude;

        /// <summary>Gets the latitude of the coordinate.</summary>
        public Latitude Latitude => latitude;

        /// <summary>Gets the longitude of the coordinate.</summary>
        public Longitude Longitude => longitude;

        /// <summary>
        /// Would be better as operator?
        /// </summary>
        public Coordinate Coordinate => new Coordinate(Longitude.Degrees, Latitude.Degrees) { Z = Altitude.GetValueOrDefault(), M = Radii };

        /// <summary>
        /// Would be better as operator?
        /// </summary>
        public Point Point => new Point(Longitude.Degrees, Latitude.Degrees) { Z = Altitude.GetValueOrDefault(), M = Radii };

        /// <summary>
        /// <see cref="Latitude"/>, <see cref="Longitude"/>, <see cref="Altitude"/>, <see cref="Radii"/>
        /// </summary>
        public double[] Coordinates => new double[] { Latitude.Degrees, Longitude.Degrees, Altitude.GetValueOrDefault(), Radii };

        /// <summary>
        /// Determines whether two specified Locations have different values.
        /// </summary>
        /// <param name="locationA">The first Location to compare, or null.</param>
        /// <param name="locationB">The second Location to compare, or null.</param>
        /// <returns>
        /// true if the value of locationA is different from the value of
        /// locationB; otherwise, false.
        /// </returns>
        public static bool operator !=(Location locationA, Location locationB) => !(locationA == locationB);
        /// <summary>
        /// Determines whether two specified Locations have the same value.
        /// </summary>
        /// <param name="locationA">The first Location to compare, or null.</param>
        /// <param name="locationB">The second Location to compare, or null.</param>
        /// <returns>
        /// true if the value of locationA is the same as the value of locationB;
        /// otherwise, false.
        /// </returns>
        public static bool operator ==(Location locationA, Location locationB)
        {
            if (object.ReferenceEquals(locationA, null))
            {
                return object.ReferenceEquals(locationB, null);
            }
            return locationA.Equals(locationB);
        }

        /// <summary>
        /// Determines whether this instance and a specified object, which must
        /// also be a Location, have the same value.
        /// </summary>
        /// <param name="obj">The Location to compare to this instance.</param>
        /// <returns>
        /// true if obj is a Location and its value is the same as this instance;
        /// otherwise, false.
        /// </returns>
        public override bool Equals(object obj) => Equals(obj as Location);

        /// <summary>
        /// Determines whether this instance and another specified Location object
        /// have the same value.
        /// </summary>
        /// <param name="other">The Location to compare to this instance.</param>
        /// <returns>
        /// true if the value of the value parameter is the same as this instance;
        /// otherwise, false.
        /// </returns>
        public bool Equals(Location other)
        {
            if (object.ReferenceEquals(other, null))
            {
                return false;
            }

            return (this.altitude == other.altitude) &&
                   (this.latitude == other.latitude) &&
                   (this.longitude == other.longitude);
        }

        /// <summary>Returns the hash code for this instance.</summary>
        /// <returns>A 32-bit signed integer hash code.</returns>
        public override int GetHashCode() => HashCode.Combine(latitude, longitude, altitude.GetValueOrDefault());

        /// <summary>
        /// Returns a string that represents the current Location in degrees,
        /// minutes and seconds form.
        /// </summary>
        /// <returns>A string that represents the current instance.</returns>
        public override string ToString() => this.ToString(null, null);

        /// <summary>
        /// Formats the value of the current instance using the specified format.
        /// </summary>
        /// <param name="format">
        /// The format to use or null to use the default format (see
        /// <see cref="Angle.ToString(string, IFormatProvider)"/>).
        /// </param>
        /// <param name="formatProvider">
        /// The provider to use to format the value or null to use the format
        /// information from the current locale setting of the operating system.
        /// </param>
        /// <returns>
        /// The value of the current instance in the specified format.
        /// </returns>
        /// <exception cref="ArgumentException">format is unknown.</exception>
        public string ToString(string format, IFormatProvider formatProvider)
        {
            if (string.IsNullOrEmpty(format))
            {
                format = "DMS";
            }

            if (formatProvider == null)
            {
                formatProvider = CultureInfo.CurrentCulture;
            }

            StringBuilder builder = new StringBuilder();
            if (format.Equals("ISO", StringComparison.OrdinalIgnoreCase))
            {
                builder.Append(this.latitude.ToString("ISO", null));
                builder.Append(this.longitude.ToString("ISO", null));
                if (this.altitude != null)
                {
                    builder.AppendFormat(CultureInfo.InvariantCulture, "{0:+0.###;-0.###}", this.altitude.Value);
                }
                builder.Append('/');
            }
            else
            {
                var parsed = Angle.ParseFormatString(format);

                builder.Append(this.latitude.ToString(format, formatProvider));
                builder.Append(' ');
                builder.Append(this.longitude.ToString(format, formatProvider));
                if (this.altitude != null)
                {
                    builder.Append(' ');
                    builder.Append(Angle.GetString(this.altitude.Value, 1, parsed.Item2, formatProvider));
                    builder.Append('m');
                }
            }
            return builder.ToString();
        }

        // These functions are based on the Aviation Formulary V1.45
        // by Ed Williams (http://williams.best.vwh.net/avform.htm)

        const double TwoPi = 2 * Math.PI;

        /// <summary>
        /// Calculates the initial course (or azimuth; the angle measured
        /// clockwise from true north) from this instance to the specified
        /// value.
        /// </summary>
        /// <param name="point">The location of the other point.</param>
        /// <returns>
        /// The initial course from this instance to the specified point.
        /// </returns>
        /// <example>
        /// The azimuth from 0,0 to 1,0 is 0 degrees. From 0,0 to 0,1 is 90
        /// degrees (due east).
        /// </example>
        public Angle Course(Location point)
        {
            double lat1 = this.latitude.Radians;
            double lon1 = this.longitude.Radians;
            double lat2 = point.latitude.Radians;
            double lon2 = point.longitude.Radians;

            double x = (Math.Cos(lat1) * Math.Sin(lat2)) -
                       (Math.Sin(lat1) * Math.Cos(lat2) * Math.Cos(lon2 - lon1));
            double tan = Math.Atan2(Math.Sin(lon2 - lon1) * Math.Cos(lat2), x);

            return Angle.FromRadians(tan % TwoPi);
        }

        /// <summary>
        /// Calculates the great circle distance, in meters, between this instance
        /// and the specified value.
        /// </summary>
        /// <param name="point">The location of the other point.</param>
        /// <returns>The great circle distance, in meters.</returns>
        /// <remarks>The antimeridian was not considered.</remarks>
        /// <exception cref="ArgumentNullException">point is null.</exception>
        public double Distance(Location point, double radius = EarthRadius)
        {
            if (point == null)
            {
                throw new ArgumentNullException(nameof(point));
            }

            double lat1 = this.latitude.Radians;
            double lon1 = this.longitude.Radians;
            double lat2 = point.latitude.Radians;
            double lon2 = point.longitude.Radians;

            double latitudeSqrd = Math.Pow(Math.Sin((lat1 - lat2) / 2), 2);
            double longitudeSqrd = Math.Pow(Math.Sin((lon1 - lon2) / 2), 2);
            double sqrt = Math.Sqrt(latitudeSqrd + (Math.Cos(lat1) * Math.Cos(lat2) * longitudeSqrd));
            double distance = 2 * Math.Asin(sqrt) / this.radii;

            if (this.altitude.HasValue && point.altitude.HasValue)
            {
                double altitudeDelta = point.altitude.Value - this.altitude.Value;
                return Math.Sqrt(Math.Pow(distance, 2) + Math.Pow(altitudeDelta, 2));
            }

            return distance;
        }

        /// <summary>
        /// Calculates a point at the specified distance along the specified
        /// radial from this instance.
        /// </summary>
        /// <param name="distance">The distance, in meters.</param>
        /// <param name="radial">
        /// The course radial from this instance, measured clockwise from north.
        /// </param>
        /// <param name="radius">The default is <see cref="EarthRadius"/></param>
        /// <returns>A Location containing the calculated point.</returns>
        /// <exception cref="ArgumentNullException">radial is null.</exception>
        /// <remarks>The antemeridian is not considered.</remarks>
        public Location GetPoint(double distance, Angle radial, double radius = EarthRadius)
        {
            if (radial == null)
            {
                throw new ArgumentNullException(nameof(radial));
            }

            double lat = this.latitude.Radians;
            double lon = this.longitude.Radians;
            distance /= radius;

            double latDist = Math.Cos(lat) * Math.Sin(distance);
            double radialLat = Math.Asin((Math.Sin(lat) * Math.Cos(distance)) +
                                         (latDist * Math.Cos(radial.Radians)));

            double y = Math.Sin(radial.Radians) * latDist;
            double x = Math.Cos(distance) - (Math.Sin(lat) * Math.Sin(radialLat));
            double atan = Math.Atan2(y, x);

            double radialLon = ((lon + atan + Math.PI) % TwoPi) - Math.PI;

            return new Location(
                Latitude.FromRadians(radialLat),
                Longitude.FromRadians(radialLon));
        }
    }

/// <summary>
    /// 
    /// </summary>
    public enum CompassDirection
    {
        N = 0,
        NNE = 22,
        NE = 45,
        ENE = 67,
        E = 90,
        ESE = 112,
        SE = 135,
        SSE = 157,
        S = 180,
        SSW = 202,
        SW = 225,
        WSW = 247,
        W = 270,
        WNW = 290,
        NW = 315,
        NNW = 337,
    }

Alternatively I find this definition of Location slightly more compact and probably more useful:

/// <summary>Represents a Latitude/Longitude/Altitude coordinate. (on Earth)</summary>
    public sealed partial class Location : Position, IEquatable<Location>, IFormattable
    {
        /// <summary>
        /// 
        /// </summary>
        /// <param name="p"></param>
        /// <returns></returns>
        public static Location FromPoint(Point p) => new Location(new Latitude(new Angle(p.Y)), new Longitude(new Angle(p.X)), p.Z);

        // Equatorial = 6378137, polar = 6356752.
        internal const int EarthRadius = 6366710; // The mean radius of the Earth in meters.

        private System.Numerics.Vector4 _data;

        /// <summary>Initializes a new instance of the Location class.</summary>
        /// <param name="latitude">The latitude of the coordinate.</param>
        /// <param name="longitude">The longitude of the coordinate.</param>
        /// <param name="altitude">
        /// The altitude, specifed in meters, of the coordinate.
        /// </param>
        /// <param name="radii"></param>
        /// <exception cref="ArgumentNullException">
        /// latitude/longitude is null.
        /// </exception>
        public Location(Latitude latitude, Longitude longitude, double altitude = 0, double radii = EarthRadius)
        {
            if (latitude == null)
            {
                throw new ArgumentNullException(nameof(latitude));
            }
            if (longitude == null)
            {
                throw new ArgumentNullException(nameof(longitude));
            }

            _data = new((float)longitude.Radians, (float)latitude.Radians, (float)altitude, (float)radii);
        }

        // XmlSerializer requires a parameterless constructor
        private Location()
        {
        }

        /// <summary>
        /// Gets the altitude of the coordinate, or null if the coordinate
        /// does not contain altitude information.
        /// </summary>
        public double? Altitude => _data.Z;

        /// <summary>
        /// The Radii of the body in which this instance is plotted
        /// </summary>
        public double Radii => _data.W;

        /// <summary>Gets the latitude of the coordinate.</summary>
        public Latitude Latitude => new(new (_data.Y));

        /// <summary>Gets the longitude of the coordinate.</summary>
        public Longitude Longitude => new(new(_data.X));

        /// <summary>
        /// Would be better as operator?
        /// </summary>
        public Coordinate Coordinate => new Coordinate(Longitude.Degrees, Latitude.Degrees) { Z = Altitude.GetValueOrDefault(), M = Radii };

        /// <summary>
        /// Would be better as operator?
        /// </summary>
        public Point Point => new Point(Longitude.Degrees, Latitude.Degrees) { Z = Altitude.GetValueOrDefault(), M = Radii };

        /// <summary>
        /// <see cref="Latitude"/>, <see cref="Longitude"/>, <see cref="Altitude"/>, <see cref="Radii"/>
        /// </summary>
        public double[] Coordinates => new double[] {  Latitude.Degrees, Longitude.Degrees, Altitude.GetValueOrDefault(), Radii };

        /// <summary>
        /// Determines whether two specified Locations have different values.
        /// </summary>
        /// <param name="locationA">The first Location to compare, or null.</param>
        /// <param name="locationB">The second Location to compare, or null.</param>
        /// <returns>
        /// true if the value of locationA is different from the value of
        /// locationB; otherwise, false.
        /// </returns>
        public static bool operator !=(Location locationA, Location locationB) => !(locationA == locationB);
        /// <summary>
        /// Determines whether two specified Locations have the same value.
        /// </summary>
        /// <param name="locationA">The first Location to compare, or null.</param>
        /// <param name="locationB">The second Location to compare, or null.</param>
        /// <returns>
        /// true if the value of locationA is the same as the value of locationB;
        /// otherwise, false.
        /// </returns>
        public static bool operator ==(Location locationA, Location locationB)
        {
            if (object.ReferenceEquals(locationA, null))
            {
                return object.ReferenceEquals(locationB, null);
            }
            return locationA.Equals(locationB);
        }

        /// <summary>
        /// Determines whether this instance and a specified object, which must
        /// also be a Location, have the same value.
        /// </summary>
        /// <param name="obj">The Location to compare to this instance.</param>
        /// <returns>
        /// true if obj is a Location and its value is the same as this instance;
        /// otherwise, false.
        /// </returns>
        public override bool Equals(object obj) => Equals(obj as Location);

        /// <summary>
        /// Determines whether this instance and another specified Location object
        /// have the same value.
        /// </summary>
        /// <param name="other">The Location to compare to this instance.</param>
        /// <returns>
        /// true if the value of the value parameter is the same as this instance;
        /// otherwise, false.
        /// </returns>
        public bool Equals(Location other)
        {
            if (object.ReferenceEquals(other, null))
            {
                return false;
            }

            return _data.Equals(other._data);
        }

        /// <summary>Returns the hash code for this instance.</summary>
        /// <returns>A 32-bit signed integer hash code.</returns>
        public override int GetHashCode() => _data.GetHashCode();

        /// <summary>
        /// Returns a string that represents the current Location in degrees,
        /// minutes and seconds form.
        /// </summary>
        /// <returns>A string that represents the current instance.</returns>
        public override string ToString() => this.ToString(null, null);

        /// <summary>
        /// Formats the value of the current instance using the specified format.
        /// </summary>
        /// <param name="format">
        /// The format to use or null to use the default format (see
        /// <see cref="Angle.ToString(string, IFormatProvider)"/>).
        /// </param>
        /// <param name="formatProvider">
        /// The provider to use to format the value or null to use the format
        /// information from the current locale setting of the operating system.
        /// </param>
        /// <returns>
        /// The value of the current instance in the specified format.
        /// </returns>
        /// <exception cref="ArgumentException">format is unknown.</exception>
        public string ToString(string format, IFormatProvider formatProvider)
        {
            if (string.IsNullOrEmpty(format))
            {
                format = "DMS";
            }

            if (formatProvider == null)
            {
                formatProvider = CultureInfo.CurrentCulture;
            }

            StringBuilder builder = new StringBuilder();
            if (format.Equals("ISO", StringComparison.OrdinalIgnoreCase))
            {
                builder.Append(this.Latitude.ToString("ISO", null));
                builder.Append(this.Longitude.ToString("ISO", null));
                builder.AppendFormat(CultureInfo.InvariantCulture, "{0:+0.###;-0.###}", this.Altitude.Value);
                builder.Append('/');
            }
            else
            {
                var parsed = Angle.ParseFormatString(format);

                builder.Append(this.Latitude.ToString(format, formatProvider));
                builder.Append(' ');
                builder.Append(this.Longitude.ToString(format, formatProvider));
                builder.Append(' ');
                builder.Append(Angle.GetString(this.Altitude.Value, 1, parsed.Item2, formatProvider));
                builder.Append('m');
            }
            return builder.ToString();
        }

        // These functions are based on the Aviation Formulary V1.45
        // by Ed Williams (http://williams.best.vwh.net/avform.htm)

        const double TwoPi = 2 * Math.PI;

        /// <summary>
        /// Calculates the initial course (or azimuth; the angle measured
        /// clockwise from true north) from this instance to the specified
        /// value.
        /// </summary>
        /// <param name="point">The location of the other point.</param>
        /// <returns>
        /// The initial course from this instance to the specified point.
        /// </returns>
        /// <example>
        /// The azimuth from 0,0 to 1,0 is 0 degrees. From 0,0 to 0,1 is 90
        /// degrees (due east).
        /// </example>
        public Angle Course(Location point)
        {
            double lat1 = this._data.Y;
            double lon1 = this._data.X;
            double lat2 = point._data.Y;
            double lon2 = point._data.X;

            double x = (Math.Cos(lat1) * Math.Sin(lat2)) -
                       (Math.Sin(lat1) * Math.Cos(lat2) * Math.Cos(lon2 - lon1));
            double tan = Math.Atan2(Math.Sin(lon2 - lon1) * Math.Cos(lat2), x);

            return Angle.FromRadians(tan % TwoPi);
        }

        /// <summary>
        /// Calculates the great circle distance, in meters, between this instance
        /// and the specified value.
        /// </summary>
        /// <param name="point">The location of the other point.</param>
        /// <returns>The great circle distance, in meters.</returns>
        /// <remarks>The antimeridian was not considered.</remarks>
        /// <exception cref="ArgumentNullException">point is null.</exception>
        public double Distance(Location point, double radius = EarthRadius)
        {
            if (point == null)
            {
                throw new ArgumentNullException(nameof(point));
            }

            //Radians
            //double lat1 = this._data.Y;
            //double lon1 = this._data.X;
            //double lat2 = point._data.Y;
            //double lon2 = point._data.X;

            //double latitudeSqrd = Math.Pow(Math.Sin((lat1 - lat2) / 2), 2);
            //double longitudeSqrd = Math.Pow(Math.Sin((lon1 - lon2) / 2), 2);
            //double sqrt = Math.Sqrt(latitudeSqrd + (Math.Cos(lat1) * Math.Cos(lat2) * longitudeSqrd));
            //double distance = 2 * Math.Asin(sqrt) / EarthRadius;

            var distance = GeometryHelper.GetDistanceInMeters(_data.X, _data.Y, point._data.X, point._data.Y, radius);

            //double altitudeDelta = point._data.Z - this._data.Z;
            //return Math.Sqrt(Math.Pow(distance, 2) + Math.Pow(altitudeDelta, 2));

            return distance;
        }

        /// <summary>
        /// Calculates a point at the specified distance along the specified
        /// radial from this instance.
        /// </summary>
        /// <param name="distance">The distance, in meters.</param>
        /// <param name="radial">
        /// The course radial from this instance, measured clockwise from north.
        /// </param>
        /// <param name="radius">The default is <see cref="EarthRadius"/></param>
        /// <returns>A Location containing the calculated point.</returns>
        /// <exception cref="ArgumentNullException">radial is null.</exception>
        /// <remarks>The antemeridian is not considered.</remarks>
        public Location GetPoint(double distance, Angle radial, double radius = EarthRadius)
        {
            if (radial == null)
            {
                throw new ArgumentNullException(nameof(radial));
            }

            double lat = this._data.Y;
            double lon = this._data.X;
            distance /= radius;

            double latDist = Math.Cos(lat) * Math.Sin(distance);
            double radialLat = Math.Asin((Math.Sin(lat) * Math.Cos(distance)) +
                                         (latDist * Math.Cos(radial.Radians)));

            double y = Math.Sin(radial.Radians) * latDist;
            double x = Math.Cos(distance) - (Math.Sin(lat) * Math.Sin(radialLat));
            double atan = Math.Atan2(y, x);

            double radialLon = ((lon + atan + Math.PI) % TwoPi) - Math.PI;

            return new Location(
                Latitude.FromRadians(radialLat),
                Longitude.FromRadians(radialLon));
        }
    }
juliusfriedman commented 3 years ago

I also modify GeodeticCalculator, GeodecticCurve and GeodecticMeasurement

public sealed class GeodeticCalculator
    {
        private const double TwoPi = Math.PI + Math.PI;
        private const double StandardTolerance = 1e-13;

        /// <summary>
        /// Calculate the destination after traveling a specified distance, and a
        /// specified starting bearing, for an initial location. This is the
        /// solution to the direct geodetic problem.
        /// </summary>
        /// <param name="ellipsoid">reference ellipsoid to use</param>
        /// <param name="start">starting location</param>
        /// <param name="startBearing">starting bearing</param>
        /// <param name="distanceMeters">distance to travel (meters)</param>
        /// <returns></returns>
        public Location CalculateEndingGlobalCoordinates(Manifold ellipsoid, Location start, Angle startBearing, double distanceMeters) => CalculateEndingGlobalCoordinates(ellipsoid, start, startBearing, distanceMeters, StandardTolerance, out var _);

        /// <summary>
        /// Calculate the destination and final bearing after traveling a specified
        /// distance, and a specified starting bearing, for an initial location.
        /// This is the solution to the direct geodetic problem.
        /// </summary>
        /// <param name="manifold">reference ellipsoid to use</param>
        /// <param name="start">starting location</param>
        /// <param name="startBearing">starting bearing</param>
        /// <param name="distanceMeters">distance to travel (meters)</param>
        /// <param name="endBearing">bearing at destination</param>
        /// <returns></returns>
        public Location CalculateEndingGlobalCoordinates(Manifold manifold, Location start, Angle startBearing, double distanceMeters, out Angle endBearing) => this.CalculateEndingGlobalCoordinates(manifold, start, startBearing, distanceMeters, StandardTolerance, out endBearing);
        internal Location CalculateEndingGlobalCoordinates(Manifold manifold, Location start, Angle startBearing, double distanceMeters, double tolerance, out Angle endBearing)
        {
            double a = manifold.SemiMajorAxisMeters;
            double b = manifold.SemiMinorAxisMeters;
            double aSquared = a * a;
            double bSquared = b * b;
            double f = manifold.Flattening;
            double phi1 = start.Latitude.Radians;
            double alpha1 = startBearing.Radians;
            double cosAlpha1 = Math.Cos(alpha1);
            double sinAlpha1 = Math.Sin(alpha1);
            double s = distanceMeters;
            double tanU1 = (1.0 - f) * Math.Tan(phi1);
            double cosU1 = 1.0 / Math.Sqrt(1.0 + tanU1 * tanU1);
            double sinU1 = tanU1 * cosU1;

            // eq. 1
            double sigma1 = Math.Atan2(tanU1, cosAlpha1);

            // eq. 2
            double sinAlpha = cosU1 * sinAlpha1;

            double sin2Alpha = sinAlpha * sinAlpha;
            double cos2Alpha = 1 - sin2Alpha;
            double uSquared = cos2Alpha * (aSquared - bSquared) / bSquared;

            // eq. 3
            double A = 1 + (uSquared / 16384) * (4096 + uSquared * (-768 + uSquared * (320 - 175 * uSquared)));

            // eq. 4
            double B = (uSquared / 1024) * (256 + uSquared * (-128 + uSquared * (74 - 47 * uSquared)));

            // iterate until there is a negligible change in sigma
            double deltaSigma;
            double sOverbA = s / (b * A);
            double sigma = sOverbA;
            double sinSigma;
            double prevSigma = sOverbA;
            double sigmaM2;
            double cosSigmaM2;
            double cos2SigmaM2;

            for (; ; )
            {
                // eq. 5
                sigmaM2 = 2.0 * sigma1 + sigma;
                cosSigmaM2 = Math.Cos(sigmaM2);
                cos2SigmaM2 = cosSigmaM2 * cosSigmaM2;
                sinSigma = Math.Sin(sigma);
                double cosSignma = Math.Cos(sigma);

                // eq. 6
                deltaSigma = B * sinSigma * (cosSigmaM2 + (B / 4.0) * (cosSignma * (-1 + 2 * cos2SigmaM2)
                    - (B / 6.0) * cosSigmaM2 * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM2)));

                // eq. 7
                sigma = sOverbA + deltaSigma;

                // break after converging to tolerance
                if (Math.Abs(sigma - prevSigma) < tolerance) break;

                prevSigma = sigma;
            }

            sigmaM2 = 2.0 * sigma1 + sigma;
            cosSigmaM2 = Math.Cos(sigmaM2);
            cos2SigmaM2 = cosSigmaM2 * cosSigmaM2;

            double cosSigma = Math.Cos(sigma);
            sinSigma = Math.Sin(sigma);

            // eq. 8
            double sinU1sinSigma_cosU1cosSigmacosAlpha1 = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1;
            double phi2 = Math.Atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,
                                     (1.0 - f) * Math.Sqrt(sin2Alpha + (sinU1sinSigma_cosU1cosSigmacosAlpha1 * sinU1sinSigma_cosU1cosSigmacosAlpha1)));

            // eq. 9
            // This fixes the pole crossing defect spotted by Matt Feemster.  When a path
            // passes a pole and essentially crosses a line of latitude twice - once in
            // each direction - the longitude calculation got messed up.  Using Atan2
            // instead of Atan fixes the defect.  The change is in the next 3 lines.
            //double tanLambda = sinSigma * sinAlpha1 / (cosU1 * cosSigma - sinU1*sinSigma*cosAlpha1);
            //double lambda = Math.Atan(tanLambda);
            double lambda = Math.Atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1);

            // eq. 10
            double C = (f / 16) * cos2Alpha * (4 + f * (4 - 3 * cos2Alpha));

            // eq. 11
            double L = lambda - (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cosSigmaM2 + C * cosSigma * (-1 + 2 * cos2SigmaM2)));

            // eq. 12
            double alpha2 = Math.Atan2(sinAlpha, -sinU1 * sinSigma + cosU1 * cosSigma * cosAlpha1);

            // build result
            Angle latitude = Angle.FromRadians(phi2);
            Angle longitude = Angle.FromRadians(start.Longitude.Radians + L);
            endBearing = Angle.FromRadians(alpha2);

            return new Location(new Latitude(latitude), new Longitude(longitude));
        }

        /// <summary>
        /// Calculate the geodetic curve between two points on a specified reference ellipsoid.
        /// This is the solution to the inverse geodetic problem.
        /// </summary>
        /// <param name="ellipsoid">reference ellipsoid to use</param>
        /// <param name="start">starting coordinates</param>
        /// <param name="end">ending coordinates </param>
        /// <returns></returns>
        public GeodeticCurve CalculateGeodeticCurve(Manifold ellipsoid, Location start, Location end) => this.CalculateGeodeticCurve(ellipsoid, start, end, StandardTolerance);
        internal GeodeticCurve CalculateGeodeticCurve(Manifold ellipsoid, Location start, Location end, double tolerance)
        {
            //
            // All equation numbers refer back to Vincenty's publication:
            // See http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
            //

            // get constants
            double a = ellipsoid.SemiMajorAxisMeters;
            double b = ellipsoid.SemiMinorAxisMeters;
            double f = ellipsoid.Flattening;

            // get parameters as radians
            double phi1 = start.Latitude.Radians;
            double lambda1 = start.Longitude.Radians;
            double phi2 = end.Latitude.Radians;
            double lambda2 = end.Longitude.Radians;

            // calculations
            double a2 = a * a;
            double b2 = b * b;
            double a2b2b2 = (a2 - b2) / b2;

            double omega = lambda2 - lambda1;

            double tanphi1 = Math.Tan(phi1);
            double tanU1 = (1.0 - f) * tanphi1;
            double U1 = Math.Atan(tanU1);
            double sinU1 = Math.Sin(U1);
            double cosU1 = Math.Cos(U1);

            double tanphi2 = Math.Tan(phi2);
            double tanU2 = (1.0 - f) * tanphi2;
            double U2 = Math.Atan(tanU2);
            double sinU2 = Math.Sin(U2);
            double cosU2 = Math.Cos(U2);

            double sinU1sinU2 = sinU1 * sinU2;
            double cosU1sinU2 = cosU1 * sinU2;
            double sinU1cosU2 = sinU1 * cosU2;
            double cosU1cosU2 = cosU1 * cosU2;

            // eq. 13
            double lambda = omega;

            // intermediates we'll need to compute 's'
            double A = 0.0;
            double B = 0.0;
            double sigma = 0.0;
            double deltasigma = 0.0;
            double lambda0;
            bool converged = false;

            for (int i = 0; i < 20; i++)
            {
                lambda0 = lambda;

                double sinlambda = Math.Sin(lambda);
                double coslambda = Math.Cos(lambda);

                // eq. 14
                double cosU1sinU2_sinU2cosU2coslambda = cosU1sinU2 - sinU1cosU2 * coslambda;
                double sin2sigma = (cosU2 * sinlambda * cosU2 * sinlambda) + (cosU1sinU2_sinU2cosU2coslambda * cosU1sinU2_sinU2cosU2coslambda);
                double sinsigma = Math.Sqrt(sin2sigma);

                // eq. 15
                double cossigma = sinU1sinU2 + (cosU1cosU2 * coslambda);

                // eq. 16
                sigma = Math.Atan2(sinsigma, cossigma);

                // eq. 17    Careful!  sin2sigma might be almost 0!
                double sinalpha = (sin2sigma == 0) ? 0.0 : cosU1cosU2 * sinlambda / sinsigma;
                double alpha = Math.Asin(sinalpha);
                double cosalpha = Math.Cos(alpha);
                double cos2alpha = cosalpha * cosalpha;

                // eq. 18    Careful!  cos2alpha might be almost 0!
                double cos2sigmam = cos2alpha == 0.0 ? 0.0 : cossigma - 2 * sinU1sinU2 / cos2alpha;
                double u2 = cos2alpha * a2b2b2;

                double cos2sigmam2 = cos2sigmam * cos2sigmam;

                // eq. 3
                A = 1.0 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)));

                // eq. 4
                B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)));

                // eq. 6
                deltasigma = B * sinsigma * (cos2sigmam + B / 4 * (cossigma * (-1 + 2 * cos2sigmam2) - B / 6 * cos2sigmam * (-3 + 4 * sin2sigma) * (-3 + 4 * cos2sigmam2)));

                // eq. 10
                double C = f / 16 * cos2alpha * (4 + f * (4 - 3 * cos2alpha));

                // eq. 11 (modified)
                lambda = omega + (1 - C) * f * sinalpha * (sigma + C * sinsigma * (cos2sigmam + C * cossigma * (-1 + 2 * cos2sigmam2)));

                if (i < 2)
                {
                    continue;
                }

                // see how much improvement we got
                double change = Math.Abs((lambda - lambda0) / lambda);

                if (change < tolerance)
                {
                    converged = true;
                    break;
                }
            }

            // eq. 19
            double s = b * A * (sigma - deltasigma);
            Angle alpha1;
            Angle alpha2;

            // didn't converge?  must be N/S
            if (!converged)
            {
                if (phi1 > phi2)
                {
                    alpha1 = Angle.Angle180;
                    alpha2 = Angle.Zero;
                }
                else if (phi1 < phi2)
                {
                    alpha1 = Angle.Zero;
                    alpha2 = Angle.Angle180;
                }
                else
                {
                    alpha1 = Angle.NaN;
                    alpha2 = Angle.NaN;
                }
            }
            else
            {
                double radians;

                // eq. 20
                radians = Math.Atan2(cosU2 * Math.Sin(lambda), (cosU1sinU2 - sinU1cosU2 * Math.Cos(lambda)));
                if (radians < 0.0) radians += TwoPi;
                alpha1 = Angle.FromRadians(radians);

                // eq. 21
                radians = Math.Atan2(cosU1 * Math.Sin(lambda), (-sinU1cosU2 + cosU1sinU2 * Math.Cos(lambda))) + Math.PI;
                if (radians < 0.0) radians += TwoPi;
                alpha2 = Angle.FromRadians(radians);
            }

            if (alpha1.Radians >= TwoPi) alpha1 = Angle.FromRadians(alpha1.Radians - TwoPi);
            if (alpha2.Radians >= TwoPi) alpha2 = Angle.FromRadians(alpha2.Radians - TwoPi);

            return new GeodeticCurve(s, alpha1, alpha2);
        }

        /// <summary>
        /// Calculate the three dimensional geodetic measurement between two positions
        /// measured in reference to a specified ellipsoid.
        /// 
        /// This calculation is performed by first computing a new ellipsoid by expanding or contracting
        /// the reference ellipsoid such that the new ellipsoid passes through the average elevation
        /// of the two positions.  A geodetic curve across the new ellisoid is calculated.  The
        /// point-to-point distance is calculated as the hypotenuse of a right triangle where the length
        /// of one side is the ellipsoidal distance and the other is the difference in elevation.
        /// </summary>
        /// <param name="manifold">reference ellipsoid to use</param>
        /// <param name="start">starting position</param>
        /// <param name="end">ending position</param>
        /// <returns></returns>
        public GeodeticMeasurement CalculateGeodeticMeasurement(Manifold manifold, Location start, Location end) => this.CalculateGeodeticMeasurement(manifold, start, end, StandardTolerance);
        internal GeodeticMeasurement CalculateGeodeticMeasurement(Manifold manifold, Location start, Location end, double tolerance)
        {
            // calculate elevation differences
            double elev1 = start.Altitude.GetValueOrDefault();
            double elev2 = end.Altitude.GetValueOrDefault();
            double elev12 = (elev1 + elev2) / 2.0;

            // calculate latitude differences
            double phi1 = start.Latitude.Radians;
            double phi2 = end.Latitude.Radians;
            double phi12 = (phi1 + phi2) / 2.0;

            // calculate a new ellipsoid to accommodate average elevation
            double refA = manifold.SemiMajorAxisMeters;
            double f = manifold.Flattening;
            double a = refA + elev12 * (1.0 + f * Math.Sin(phi12));
            Manifold ellipsoid = Manifold.FromAAndF(a, f);

            // calculate the curve at the average elevation
            GeodeticCurve averageCurve = CalculateGeodeticCurve(ellipsoid, start, end, tolerance);

            // return the measurement
            return new GeodeticMeasurement(averageCurve, elev2 - elev1);
        }
    }

Serializable]
    public struct GeodeticCurve : IEquatable<GeodeticCurve>, ISerializable
    {
        /// <summary>
        /// Create a new GeodeticCurve.
        /// </summary>
        /// <param name="ellipsoidalDistanceMeters">ellipsoidal distance in meters</param>
        /// <param name="azimuth">azimuth in degrees</param>
        /// <param name="reverseAzimuth">reverse azimuth in degrees</param>
        public GeodeticCurve(double ellipsoidalDistanceMeters, Angle azimuth, Angle reverseAzimuth)
        {
            EllipsoidalDistanceMeters = ellipsoidalDistanceMeters;
            Azimuth = azimuth;
            ReverseAzimuth = reverseAzimuth;
        }

        /// <summary>Ellipsoidal distance (in meters).</summary>
        public double EllipsoidalDistanceMeters { get; }

        /// <summary>
        /// Get the azimuth.  This is angle from north from start to end.
        /// </summary>
        public Angle Azimuth { get; }

        /// <summary>
        /// Get the reverse azimuth.  This is angle from north from end to start.
        /// </summary>
        public Angle ReverseAzimuth { get; }

        // TODO: consider just leaving it at ellipsoidal distance...
        public static int GetHashCode(GeodeticCurve value) => HashCode.Combine(value.EllipsoidalDistanceMeters, value.Azimuth, value.ReverseAzimuth);

        public static bool Equals(GeodeticCurve first, GeodeticCurve second) => first.EllipsoidalDistanceMeters == second.EllipsoidalDistanceMeters &&
                                                                                first.Azimuth == second.Azimuth &&
                                                                                first.ReverseAzimuth == second.ReverseAzimuth;
        /// <summary>
        /// 
        /// </summary>
        /// <param name="value"></param>
        /// <returns></returns>
        public static string ToString(GeodeticCurve value) => $"GeodeticCurve[EllipsoidalDistanceMeters={value.EllipsoidalDistanceMeters}, Azimuth={value.Azimuth}, ReverseAzimuth={value.ReverseAzimuth}]";

        /// <summary>
        /// 
        /// </summary>
        /// <returns></returns>
        public override int GetHashCode() => GetHashCode(this);

        /// <summary>
        /// 
        /// </summary>
        /// <param name="obj"></param>
        /// <returns></returns>
        public override bool Equals(object obj) => obj is GeodeticCurve && Equals(this, (GeodeticCurve)obj);

        /// <summary>
        /// 
        /// </summary>
        /// <param name="other"></param>
        /// <returns></returns>
        public bool Equals(GeodeticCurve other) => Equals(this, other);

        /// <summary>
        /// Get curve as a string.
        /// </summary>
        /// <returns></returns>
        public override string ToString() => ToString(this);

        #region Serialization / Deserialization

        private GeodeticCurve(SerializationInfo info, StreamingContext context)
        {
            this.EllipsoidalDistanceMeters = info.GetDouble("ellipsoidalDistanceMeters");

            double azimuthRadians = info.GetDouble("azimuthRadians");
            double reverseAzimuthRadians = info.GetDouble("reverseAzimuthRadians");

            this.Azimuth = Angle.FromRadians(azimuthRadians);
            this.ReverseAzimuth = Angle.FromRadians(reverseAzimuthRadians);
        }

        void ISerializable.GetObjectData(SerializationInfo info, StreamingContext context)
        {
            info.AddValue("ellipsoidalDistanceMeters", this.EllipsoidalDistanceMeters);

            info.AddValue("azimuthRadians", this.Azimuth.Radians);
            info.AddValue("reverseAzimuthRadians", this.ReverseAzimuth.Radians);
        }

        #endregion

        #region Operators

        public static bool operator ==(GeodeticCurve lhs, GeodeticCurve rhs) => Equals(lhs, rhs);
        public static bool operator !=(GeodeticCurve lhs, GeodeticCurve rhs) => !Equals(lhs, rhs);

        #endregion
    }

[Serializable]
    public struct GeodeticMeasurement : IEquatable<GeodeticMeasurement>, ISerializable
    {
        /// <summary>
        /// Creates a new instance of GeodeticMeasurement.
        /// </summary>
        /// <param name="averageCurve">the geodetic curve as measured at the average elevation between two points</param>
        /// <param name="elevationChangeMeters">the change in elevation, in meters, going from the starting point to the ending point</param>
        public GeodeticMeasurement(GeodeticCurve averageCurve, double elevationChangeMeters)
        {
            double ellipsoidalDistanceMeters = averageCurve.EllipsoidalDistanceMeters;

            this.AverageCurve = averageCurve;
            this.ElevationChangeMeters = elevationChangeMeters;
            this.PointToPointDistanceMeters = Math.Sqrt((ellipsoidalDistanceMeters * ellipsoidalDistanceMeters) + (elevationChangeMeters * elevationChangeMeters));
        }

        /// <summary>
        /// Get the average geodetic curve.  This is the geodetic curve as measured
        /// at the average elevation between two points.
        /// </summary>
        public GeodeticCurve AverageCurve { get; }

        /// <summary>
        /// Get the ellipsoidal distance (in meters).  This is the length of the average geodetic
        /// curve.  For actual point-to-point distance, use PointToPointDistance property.
        /// </summary>
        public double EllipsoidalDistanceMeters => this.AverageCurve.EllipsoidalDistanceMeters;

        /// <summary>
        /// Get the azimuth.  This is angle from north from start to end.
        /// </summary>
        public Angle Azimuth => this.AverageCurve.Azimuth;

        /// <summary>
        /// Get the reverse azimuth.  This is angle from north from end to start.
        /// </summary>
        public Angle ReverseAzimuth => this.AverageCurve.ReverseAzimuth;

        /// <summary>
        /// Get the elevation change, in meters, going from the starting to the ending point.
        /// </summary>
        public double ElevationChangeMeters { get; }

        /// <summary>
        /// Get the distance travelled, in meters, going from one point to the next.
        /// </summary>
        public double PointToPointDistanceMeters { get; }

        // p2p is a derived metric, no need to test.
        public static int GetHashCode(GeodeticMeasurement value) => HashCode.Combine(value.AverageCurve, value.ElevationChangeMeters);

        public static bool Equals(GeodeticMeasurement first, GeodeticMeasurement second) => first.AverageCurve == second.AverageCurve &&
                                                                                            first.ElevationChangeMeters == second.ElevationChangeMeters;
        public static string ToString(GeodeticMeasurement value) => $"GeodeticMeasurement[AverageCurve={value.AverageCurve}, ElevationChangeMeters={value.ElevationChangeMeters}, PointToPointDistanceMeters={value.PointToPointDistanceMeters}]";

        public override int GetHashCode() => GetHashCode(this);
        public override bool Equals(object obj) => obj is GeodeticMeasurement && Equals(this, (GeodeticMeasurement)obj);
        public bool Equals(GeodeticMeasurement other) => Equals(this, other);

        /// <summary>
        /// Get the GeodeticMeasurement as a string
        /// </summary>
        /// <returns></returns>
        public override string ToString() => ToString(this);

        #region Serialization / Deserialization

        private GeodeticMeasurement(SerializationInfo info, StreamingContext context)
        {
            double elevationChangeMeters = this.ElevationChangeMeters = info.GetDouble("elevationChangeMeters");

            double ellipsoidalDistanceMeters = info.GetDouble("averageCurveEllipsoidalDistanceMeters");
            double azimuthRadians = info.GetDouble("averageCurveAzimuthRadians");
            double reverseAzimuthRadians = info.GetDouble("averageCurveReverseAzimuthRadians");

            this.AverageCurve = new GeodeticCurve(ellipsoidalDistanceMeters, Angle.FromRadians(azimuthRadians), Angle.FromRadians(reverseAzimuthRadians));
            this.PointToPointDistanceMeters = Math.Sqrt((ellipsoidalDistanceMeters * ellipsoidalDistanceMeters) + (elevationChangeMeters * elevationChangeMeters));
        }

        void ISerializable.GetObjectData(SerializationInfo info, StreamingContext context)
        {
            info.AddValue("elevationChangeMeters", this.ElevationChangeMeters);

            info.AddValue("averageCurveEllipsoidalDistanceMeters", this.AverageCurve.EllipsoidalDistanceMeters);
            info.AddValue("averageCurveAzimuthRadians", this.AverageCurve.Azimuth.Radians);
            info.AddValue("averageCurveReverseAzimuthRadians", this.AverageCurve.ReverseAzimuth.Radians);
        }

        #endregion

        #region Operators

        public static bool operator ==(GeodeticMeasurement lhs, GeodeticMeasurement rhs) => Equals(lhs, rhs);
        public static bool operator !=(GeodeticMeasurement lhs, GeodeticMeasurement rhs) => !Equals(lhs, rhs);

        #endregion
    }