besselearth.cs

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

                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 W2 = Math.Sqrt(1 - earthParam.e1 * earthParam.e1 * Math.Sin(m_dEndLatitude * Math.PI / 180) * Math.Sin(m_dEndLatitude * Math.PI / 180));
            double sinu1 = Math.Sin(m_dStartLatitude * Math.PI / 180) * Math.Sqrt(1 - earthParam.e1 * earthParam.e1) / W1;
            double sinu2 = Math.Sin(m_dEndLatitude * Math.PI / 180) * Math.Sqrt(1 - earthParam.e1 * earthParam.e1) / W2;
            double cosu1 = Math.Cos(m_dStartLatitude * Math.PI / 180) / W1;
            double cosu2 = Math.Cos(m_dEndLatitude * Math.PI / 180) / W2;
            double l = m_dEndLongitude - m_dStartLongitude;
            double a1 = sinu1 * sinu2;
            double a2 = cosu1 * cosu2;
            double b1 = cosu1 * sinu2;
            double b2 = sinu1 * cosu2;
            // 用逐次趋近法同时计算起点大地方位角、球面长度及经差。
            double theta = 0;
            double theta0, lamuda, sinA0, cos2A0, x, zita;
            double p, q;
            do
            {
                theta0 = theta;
                lamuda = l + theta;
                p = cosu2 * Math.Sin(lamuda * Math.PI / 180);
                q = b1 - b2 * Math.Cos(lamuda * Math.PI / 180);
                m_dA1 = Math.Atan2(p, q) * 180 / Math.PI;
                if (p > 0)
                {
                    if (q > 0)
                    {
                        m_dA1 = Math.Abs(m_dA1);
                    }
                    else
                    {
                        m_dA1 = 180 - Math.Abs(m_dA1);
                    }
                }
                else
                {
                    if (q > 0)
                    {
                        m_dA1 = 360 - Math.Abs(m_dA1);
                    }
                    else
                    {
                        m_dA1 = 180 + Math.Abs(m_dA1);
                    }
                }
                double sinzita = p * Math.Sin(m_dA1 * Math.PI / 180) + q * Math.Cos(m_dA1 * Math.PI / 180);
                double coszita = a1 + a2 * Math.Cos(lamuda * Math.PI / 180);
                zita = Math.Atan2(sinzita, coszita) * 180 / Math.PI;
                if (coszita > 0)
                {
                    zita = Math.Abs(zita);
                }
                else
                {
                    zita = 180 - Math.Abs(zita);
                }
                sinA0 = cosu1 * Math.Sin(m_dA1 * Math.PI / 180);
                cos2A0 = 1 - sinA0 * sinA0;
                x = 2 * a1 - cos2A0 * Math.Cos(zita * Math.PI / 180);
                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 BetaP = 2 * Beta;
                theta = (Alpha * (zita * Math.PI / 180) - BetaP * x * Math.Sin(zita * Math.PI / 180)) * sinA0 * 180 / Math.PI;
            }
            while (theta == theta0 || Math.Abs(theta - theta0) < dThreshold);
            lamuda = l + theta;
            // 计算系数A, B'', C''
            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 BPP = 2 * B / cos2A0;
            double CPP = 2 * C / cos2A0 / cos2A0;
            // 计算大地线长度
            double y = (cos2A0 * cos2A0 - 2 * x * x) * Math.Cos(zita * Math.PI / 180);
            m_dLength = A * zita * Math.PI / 180 + (BPP * x + CPP * y) * Math.Sin(zita * Math.PI / 180);
            // 计算反方位角
            double p2 = cosu1 * Math.Sin(lamuda * Math.PI / 180);
            double q2 = b1 * Math.Cos(lamuda * Math.PI / 180) - b2;
            m_dA2 = Math.Atan2(p2, q2) * 180 / Math.PI;
            if (p2 > 0)
            {
                if (q2 > 0)
                {
                    m_dA2 = 180 + Math.Abs(m_dA2);
                }
                else
                {
                    m_dA2 = 360 - Math.Abs(m_dA2);
                }
            }
            else
            {
                if (q2 > 0)
                {
                    m_dA2 = Math.Abs(m_dA2);
                }
                else
                {
                    m_dA2 = 180 - Math.Abs(m_dA2);
                }
            }
        }
    }

    /// <remarks>地球椭球参数</remarks>
    public class EarthParameter
    {
        /// <remarks>椭圆的长半轴（米）</remarks>
        private double m_dA;
        /// <remarks>椭圆的短半轴（米）</remarks>
        private double m_dB;

        /// <remarks>输入椭球长短半轴</remarks>
        public EarthParameter(double dA, double dB)
        {
            m_dA = dA;
            m_dB = dB;
        }

        /// <remarks>输入椭球类型</remarks>
        public EarthParameter(EnumEarthParameter typeEarthParameter)
        {
            switch (typeEarthParameter)
            {
                case EnumEarthParameter.International1975:
                    {
                        // 1975国际椭球体参数
                        m_dA = 6378140.0;
                        m_dB = 6356755.2881575287;
                        break;
                    }
                case EnumEarthParameter.KeLaSuoFuSiJi:
                    {
                        // 克拉索夫斯基椭球体参数
                        m_dA = 6378245.0;
                        m_dB = 6356863.0187730473;
                        break;
                    }
                case EnumEarthParameter.WGS84:
                    {
                        // WGS-84椭球体参数
                        m_dA = 6378137.0;
                        m_dB = 6356752.3142;
                        break;
                    }
            }
        }

        /// <remarks>获取/设置椭圆的长半轴（米）</remarks>
        public double a
        {
            get
            {
                return m_dA;
            }
            set
            {
                m_dA = value;
            }
        }

        /// <remarks>获取/设置椭圆的短半轴（米）</remarks>
        public double b
        {
            get
            {
                return m_dB;
            }
            set
            {
                m_dB = value;
            }
        }

        /// <remarks>获取引入的参数</remarks>
        public double c
        {
            get
            {
                return m_dA * m_dA / m_dB;
            }
        }

        /// <remarks>获取椭圆的第一偏心率</remarks>
        public double e1
        {
            get
            {
                return Math.Sqrt(m_dA * m_dA - m_dB * m_dB) / m_dA;
            }
        }

        /// <remarks>获取椭圆的第二偏心率</remarks>
        public double e2
        {
            get
            {
                return Math.Sqrt(m_dA * m_dA - m_dB * m_dB) / m_dB;
            }
        }

        /// <remarks>获取椭圆的扁率</remarks>
        public double alpha
        {
            get
            {
                return (m_dA - m_dB) / m_dA;
            }
        }
    }

    /// <remarks>椭球体类型</remarks>
    public enum EnumEarthParameter
    {
        /// <remarks>克拉索夫斯基椭球体</remarks>
        KeLaSuoFuSiJi,
        /// <remarks>1975国际椭球体</remarks>
        International1975,
        /// <remarks>WGS-84椭球体</remarks>
        WGS84,
    }
}