src/cs/Geodetic/Geodetic.cs

using System;
using static System.Math;

namespace Prelude.Geodetic {
    public static class Datum {
        public const double SemiMajorAxis = 6378137.0; // a (in meters)
        public const double SemiMinorAxis = 6356752.31424518; // b (in meters)
        public const double FlatteningFactor = 298.257223563; // 1/f
        public const double LinearEccentricity = 521854.00842339; // E (in meters)
        public const double Eccentricity = 0.0818191908426215; // e
        public const double EccentricitySquared = 0.006694379990141; // e2
        public const double Radius = 6371001; // mean radius (in meters)
        public const double RadiusAuthalic = 6371007.1810; // radius constant surface area (in meters)
    }
    public class Coordinate {
        private double _Latitude;
        private double _Longitude;
        private string[] _Hemisphere = new string[] { "N", "E" };
        public string[] Hemisphere {
            get {
                return _Hemisphere;
            }
            private set {
                _Hemisphere = value;
            }
        }
        public double Latitude {
            get {
                return _Latitude;
            }
            set {
                Hemisphere[0] = value < 0 ? "S" : "N";
                _Latitude = value;
            }
        }
        public double Longitude {
            get {
                return _Longitude;
            }
            set {
                Hemisphere[1] = value < 0 ? "W" : "E";
                _Longitude = value;
            }
        }
        public double Height;
        private static double ToDegree(double value) => value * (180 / PI);
        private static double ToRadian(double value) => value * (PI / 180);
        public static double[] ToSexagesimal(double value) {
            double fractionalPart = Abs(value - Truncate(value));
            double degree = Truncate(value);
            double minute = Truncate(fractionalPart * 60);
            double second = Round(((fractionalPart * 60) - minute) * 60, 2);
            return new double[] { degree, minute, second };
        }
        public Coordinate() {
            Latitude = 0.0;
            Longitude = 0.0;
            Height = 0.0;
        }
        public Coordinate(double latitude, double longitude, double height = 0.0) {
            Latitude = latitude;
            Longitude = longitude;
            Height = height;
        }
        public static Coordinate FromCartesian(double x, double y, double z) {
            double[] geodetic = ToGeodetic(x, y, z);
            double latitude = geodetic[0];
            double longitude = geodetic[1];
            double height = geodetic[2];
            return new Coordinate(latitude, longitude, height);
        }
        public static double[] ToGeodetic(double x, double y, double z) {
            double a = Datum.SemiMajorAxis;
            double b = Datum.SemiMinorAxis;
            double E = Datum.LinearEccentricity;
            double E2 = Pow(E, 2);
            double x2 = Pow(x, 2), y2 = Pow(y, 2), z2 = Pow(z, 2);
            double r2 = x2 + y2 + z2;
            double Q = Sqrt(x2 + y2);
            double u = Sqrt(0.5 * ((r2 - E2) + Sqrt(Pow(r2 - E2, 2) + (4 * E2 * z2))));
            double u2 = Pow(u, 2);
            double beta = Atan(Sqrt(u2 + E2) * z / (u * Q));
            double latitude = Atan(a / b * Tan(beta));
            double longitude = Atan2(y, x);
            double height = Sqrt(Pow(z - (b * Sin(beta)), 2) + Pow(Q - (a * Cos(beta)), 2));
            return new double[] { ToDegree(latitude), ToDegree(longitude), height };
        }
        public static double[] ToCartesian(double latitude, double longitude, double height = 0) {
            double a = Datum.SemiMajorAxis;
            double e2 = Datum.EccentricitySquared;
            double lat = ToRadian(latitude);
            double lon = ToRadian(longitude);
            double h = height;
            double v = a / Sqrt(1 - (e2 * Pow(Sin(lat), 2)));
            double x = Cos(lat) * Cos(lon) * (v + h);
            double y = Cos(lat) * Sin(lon) * (v + h);
            double z = Sin(lat) * ((v * (1 - e2)) + h);
            return new double[] { x, y, z };
        }
        public override string ToString() {
            double latitude = Abs(Latitude);
            double longitude = Abs(Longitude);
            string NS = Hemisphere[0];
            string WE = Hemisphere[1];
            string[] lat = Array.ConvertAll(ToSexagesimal(latitude), Convert.ToString);
            string[] lon = Array.ConvertAll(ToSexagesimal(longitude), Convert.ToString);
            return $"{lat[0]}°{lat[1]}'{lat[2]}\"{NS} {lon[0]}°{lon[1]}'{lon[2]}\"{WE}";
        }
    }
}