besselearth.cs

白塞尔大地主题解算 包括正算与反算。可由经纬度计算距离

using System;
using System.Collections.Generic;
using System.Text;

namespace BesselEarth
{
    /// <summary>
    /// 白塞尔大地主体解算
    /// </summary>
    /// <remarks>
    /// 白塞尔大地主体解算
    /// 作者：Polaris
    /// 来源：开源盛世-源代码下载网
    ///       http://www.vscodes.com
    /// 版权所有，转载，引用代码请保留以上信息！
    /// </remarks>
    class BesselEarth
    {
        /// <remarks>大地起点经度</remarks>
        private double m_dStartLongitude = 0;
        /// <remarks>大地终点经度</remarks>
        private double m_dEndLongitude = 0;
        /// <remarks>大地起点纬度</remarks>
        private double m_dStartLatitude = 0;
        /// <remarks>大地终点纬度</remarks>
        private double m_dEndLatitude = 0;
        /// <remarks>大地方位角A1</remarks>
        private double m_dA1 = 0;
        /// <remarks>大地方位角A2</remarks>
        private double m_dA2 = 0;
        /// <remarks>大地线长度（米）</remarks>
        private double m_dLength = 0;
        /// <remarks>克拉索夫斯基椭球体</remarks>
        private EarthParameter m_earthKeLaSuoFuSiJi = new EarthParameter(EnumEarthParameter.KeLaSuoFuSiJi);
        /// <remarks>1975国际椭球体</remarks>
        private EarthParameter m_earthInternational1975 = new EarthParameter(EnumEarthParameter.International1975);
        /// <remarks>WGS-84椭球体</remarks>
        private EarthParameter m_earthWGS84 = new EarthParameter(EnumEarthParameter.WGS84);

        /// <remarks>获取/设置大地方位角A1</remarks>
        public double A1
        {
            get
            {
                return m_dA1;
            }
            set
            {
                m_dA1 = value;
            }
        }

        /// <remarks>获取/设置大地方位角A2</remarks>
        public double A2
        {
            get
            {
                return m_dA2;
            }
            set
            {
                m_dA2 = value;
            }
        }

        /// <remarks>获取/设置大地线长度（米）</remarks>
        public double Length
        {
            get
            {
                return m_dLength;
            }
            set
            {
                m_dLength = value;
            }
        }

        /// <remarks>获取/设置大地起点纬度</remarks>
        public double StartLatitude
        {
            get
            {
                return m_dStartLatitude;
            }
            set
            {
                m_dStartLatitude = value;
            }
        }

        /// <remarks>获取/设置大地起点经度</remarks>
        public double StartLongitude
        {
            get
            {
                return m_dStartLongitude;
            }
            set
            {
                m_dStartLongitude = value;
            }
        }

        /// <remarks>获取/设置大地终点纬度</remarks>
        public double EndLatitude
        {
            get
            {
                return m_dEndLatitude;
            }
            set
            {
                m_dEndLatitude = value;
            }
        }

        /// <remarks>获取/设置大地终点经度</remarks>
        public double EndLongitude
        {
            get
            {
                return m_dEndLongitude;
            }
            set
            {
                m_dEndLongitude = value;
            }
        }
        /// <remarks>正算</remarks>
        public void Computation(EnumEarthParameter typeEarthParameter)
        {
            // 获取地球椭球参数
            EarthParameter earthParam;
            switch(typeEarthParameter)
            {
                case EnumEarthParameter.International1975:
                    {
                        earthParam = m_earthInternational1975;
                        break;
                    }
                case EnumEarthParameter.KeLaSuoFuSiJi:
                    {
                        earthParam = m_earthKeLaSuoFuSiJi;
                        break;
                    }
                case EnumEarthParameter.WGS84:
                    {
                        earthParam = m_earthWGS84;
                        break;
                    }
                default:
                    {
                        earthParam = m_earthWGS84;
                        break;
                    }
            }
            // 计算起点的归化纬度
            double W1 = Math.Sqrt(1 - earthParam.e1 * earthParam.e1 * Math.Sin(m_dStartLatitude * Math.PI / 180) * Math.Sin(m_dStartLatitude * Math.PI / 180));
            double sinu1 = Math.Sin(m_dStartLatitude * Math.PI / 180) * Math.Sqrt(1 - earthParam.e1 * earthParam.e1) / W1;
            double cosu1 = Math.Cos(m_dStartLatitude * Math.PI / 180) / W1;
            // 计算辅助函数值
            double sinA0 = cosu1 * Math.Sin(m_dA1 * Math.PI / 180);
            double cotq1 = cosu1 * Math.Cos(m_dA1 * Math.PI / 180);
            double sin2q1 = 2 * cotq1 / (Math.Pow(cotq1, 2) + 1);
            double cos2q1 = (Math.Pow(cotq1, 2) - 1) / (Math.Pow(cotq1, 2) + 1);
            // 计算系数A,B,C及Alpha, Beta的值。
            double cos2A0 = 1 - Math.Pow(sinA0, 2);
            double k2 = earthParam.e2 * earthParam.e2 * cos2A0;
            double A = earthParam.b * (1 + k2 / 4 - 3 * k2 * k2 / 64 + 5 * k2 * k2 * k2 / 256);
            double B = earthParam.b * (k2 / 8 - k2 * k2 / 32 + 15 * k2 * k2 * k2 / 1024);
            double C = earthParam.b * (k2 * k2 / 128 - 3 * k2 * k2 * k2 / 512);
            double e2 = earthParam.e1 * earthParam.e1;
            double Alpha = (e2 / 2 + e2 * e2 / 8 + e2 * e2 * e2 / 16) - (e2 * e2 / 16 + e2 * e2 * e2 / 16) * cos2A0 + (3 * e2 * e2 * e2 / 128) * cos2A0 * cos2A0;
            double Beta = (e2 * e2 / 32 + e2 * e2 * e2 / 32) * cos2A0 - (e2 * e2 * e2 / 64) * cos2A0 * cos2A0;
            // 计算球面长度
            double q0 = (m_dLength - (B + C * cos2q1) * sin2q1) / A;
            double sin2q1q0 = sin2q1 * Math.Cos(2 * q0) + cos2q1 * Math.Sin(2 * q0);
            double cos2q1q0 = cos2q1 * Math.Cos(2 * q0) - sin2q1 * Math.Sin(2 * q0);
            double q = q0 + (B + 5 * C * cos2q1q0) * sin2q1q0 / A;
            // 计算经度差改正数
            double theta = (Alpha * q + Beta * (sin2q1q0 - sin2q1)) * sinA0;
            // 计算终点大地坐标及大地方位角
            double sinu2 = sinu1 * Math.Cos(q) + cosu1 * Math.Cos(m_dA1 * Math.PI / 180) * Math.Sin(q);
            m_dEndLatitude = Math.Atan2(sinu2, (Math.Sqrt(1 - earthParam.e1 * earthParam.e1) * Math.Sqrt(1 - sinu2 * sinu2))) * 180 / Math.PI;
            double lamuda = Math.Atan2(Math.Sin(m_dA1 * Math.PI / 180) * Math.Sin(q), (cosu1 * Math.Cos(q) - sinu1 * Math.Sin(q) * Math.Cos(m_dA1 * Math.PI / 180))) * 180 / Math.PI;
            if (Math.Sin(m_dA1 * Math.PI / 180) > 0)
            {
                if (Math.Sin(m_dA1 * Math.PI / 180) * Math.Sin(q) / (cosu1 * Math.Cos(q) - sinu1 * Math.Sin(q) * Math.Cos(m_dA1 * Math.PI / 180)) > 0)
                {
                    lamuda = Math.Abs(lamuda);
                }
                else
                {
                    lamuda = 180 - Math.Abs(lamuda);
                }
            }
            else
            {
                if (Math.Sin(m_dA1 * Math.PI / 180) * Math.Sin(q) / (cosu1 * Math.Cos(q) - sinu1 * Math.Sin(q) * Math.Cos(m_dA1 * Math.PI / 180)) > 0)
                {
                    lamuda = Math.Abs(lamuda) - 180;
                }
                else
                {
                    lamuda = -Math.Abs(lamuda);
                }
            }
            m_dEndLongitude = m_dStartLongitude + lamuda - theta * 180 / Math.PI;
            m_dA2 = Math.Atan2(cosu1 * Math.Sin(m_dA1 * Math.PI / 180), (cosu1 * Math.Cos(q) * Math.Cos(m_dA1 * Math.PI / 180) - sinu1 * Math.Sin(q))) * 180 / Math.PI;
            if (Math.Sin(m_dA1 * Math.PI / 180) > 0)
            {
                if (cosu1 * Math.Sin(m_dA1 * Math.PI / 180) / (cosu1 * Math.Cos(q) * Math.Cos(m_dA1 * Math.PI / 180) - sinu1 * Math.Sin(q)) > 0)
                {
                    m_dA2 = 180 + Math.Abs(m_dA2);
                }
                else
                {
                    m_dA2 = 360 - Math.Abs(m_dA2);
                }
            }
            else
            {
                if (cosu1 * Math.Sin(m_dA1 * Math.PI / 180) / (cosu1 * Math.Cos(q) * Math.Cos(m_dA1 * Math.PI / 180) - sinu1 * Math.Sin(q)) > 0)
                {
                    m_dA2 = Math.Abs(m_dA2);
                }
                else
                {
                    m_dA2 = 180 - Math.Abs(m_dA2);
                }
            }
        }

        /// <remarks>反算</remarks>
        public void InverseComputation(EnumEarthParameter typeEarthParameter)
        {
            InverseComputation(typeEarthParameter, 0.000000000005);
        }

        /// <remarks>反算</remarks>
        public void InverseComputation(EnumEarthParameter typeEarthParameter, double dThreshold)
        {
            // 获取地球椭球参数
            EarthParameter earthParam;
            switch (typeEarthParameter)
            {