csharp源码.txt

根据地球上两点经纬度（以度分秒为单位）

       /// <summary>
        /// 贝赛尔大地问题正解
        /// </summary>
        /// <param name="B1">已知点纬度</param>
        /// <param name="L1">已知点经度</param>
        /// <param name="S">大地线长</param>
        /// <param name="A1">大地方位角</param>
        /// <param name="B2">待求点纬度</param>
        /// <param name="L2">待求点经度</param>
        /// <param name="A2">大地反方位角</param>
        /// <param name="a">参考椭球长半轴</param>
        /// <param name="f">参考椭球扁率倒数</param>
        public static void Bessel_PntSA_Pnt(double B1, double L1, double S, double A1, out double B2, out double L2, out double A2, double a, double f)
        {
            double ee = (2 * f - 1) / f / f;       //第一偏心率的平方
            double ee2 = ee / (1 - ee);            //第二偏心率的平方

            double u1 = Math.Atan(Math.Sqrt(1 - ee) * Math.Tan(B1 * Math.PI / 180));
            double m = Math.Cos(u1) * Math.Sin(A1);
            m = Math.Atan(m / Math.Sqrt(1 - m * m));

            double M_ = Math.Atan(Math.Tan(u1) / Math.Cos(A1));
            if (M_ < 0) M_ += Math.PI;

            double KK = ee2 * Math.Pow(Math.Cos(m), 2);
            double alpha = Math.Sqrt(1 + ee2) * (1 - KK / 4 + 7 * KK * KK / 64 - 15 * Math.Pow(KK, 3) / 256) / a;
            double beta = KK / 4 - KK * KK / 8 + 37 * Math.Pow(KK, 3) / 512;
            double gama = KK * KK * (1 - KK) / 128;

            double sigma, temp;
            sigma = alpha * S;
            temp = 0;
            while (Math.Abs(temp - sigma) > 0.0000000001)
            {
                temp = sigma;
                sigma = alpha * S + beta * Math.Sin(sigma) * Math.Cos(2 * M_ + sigma) + gama * Math.Sin(2 * sigma) * Math.Cos(4 * M_ + 2 * sigma);
            }

            A2 = Math.Atan(Math.Tan(m) / Math.Cos(M_ + sigma));
            if (Math.Cos(M_) > 0) A2 += Math.PI;
            if (m < 0 && Math.Cos(M_) < 0) A2 += 2 * Math.PI;

            double u2 = Math.Atan(-Math.Cos(A2) * Math.Tan(M_ + sigma));

            double lamda1;
            lamda1 = Math.Atan(Math.Sin(u2) * Math.Tan(A2));
            if (Math.Cos(M_) < 0) lamda1 += Math.PI;
            if (m < 0 && Math.Cos(M_) > 0) lamda1 += 2 * Math.PI;

            double lamda2;
            lamda2 = Math.Atan(Math.Sin(u1) * Math.Tan(A1));
            if (m > 0)
            {
                if (Math.Cos(M_ + sigma) < 0)
                    lamda2 += Math.PI;
                else if (Math.Sin(M_ + sigma) < 0)
                    lamda2 += 2 * Math.PI;
            }
            else
            {
                if (Math.Cos(2 * Math.PI - M_ - sigma) < 0)
                    lamda2 += Math.PI;
                else if (Math.Sin(M_ + sigma) < 0)
                    lamda2 += 2 * Math.PI;
            }

            B2 = Math.Atan(Math.Sqrt(1 + ee2) * Math.Tan(u2)) * 180 / Math.PI;

            KK = ee * Math.Pow(Math.Cos(m), 2);
            alpha = ee / 2 + ee * ee / 8 + Math.Pow(ee, 3) / 16 - ee * (1 + ee) * KK / 16 + 3 * ee * KK * KK / 128;
            gama = ee * KK * KK / 256;

            L2 = lamda1 + lamda2 - Math.Sin(m) * (alpha * sigma + beta * Math.Sin(sigma) * Math.Cos(2 * M_ + sigma) + gama * Math.Sin(2 * sigma) * Math.Cos(4 * M_ + 2 * sigma));
            L2 = L1 + L2 * 180 / Math.PI;
        }

        /// <summary>
        /// 由两点大地坐标求解大地方位角
        /// </summary>
        /// <param name="B1">起点纬度</param>
        /// <param name="L1">起点经度</param>
        /// <param name="B2">末点纬度</param>
        /// <param name="L2">末点经度</param>
        /// <param name="a">参考椭球长半轴</param>
        /// <param name="f">参考椭球扁率倒数</param>
        /// <returns>两点间的大地方位角</returns>
        public static double Bessel_BL_A(double B1, double L1, double B2, double L2, double a, double f)
        {
            if (L1 == L2) return 0;

            double ee = (2 * f - 1) / f / f;       //第一偏心率的平方

            double u1 = Math.Atan(Math.Sqrt(1 - ee) * Math.Tan(B1 * Math.PI / 180));
            double u2 = Math.Atan(Math.Sqrt(1 - ee) * Math.Tan(B2 * Math.PI / 180));
            double dL = (L2 - L1) * Math.PI / 180;

            double sigma = Math.Sin(u1) * Math.Sin(u2) + Math.Cos(u1) * Math.Cos(u2) * Math.Cos(dL);
            sigma = Math.Atan(Math.Sqrt(1 - sigma * sigma) / sigma);
            if (sigma <= 0) sigma += Math.PI;

            double m = Math.Cos(u1) * Math.Cos(u2) * Math.Sin(dL) / Math.Sin(sigma);
            m = Math.Atan(m / Math.Sqrt(1 - m * m));

            double KK = ee * Math.Pow(Math.Cos(m), 2);
            double alpha = ee / 2 + ee * ee / 8 + Math.Pow(ee, 3) / 16 - ee * (1 + ee) * KK / 16 + 3 * ee * KK * KK / 128;
            double lamda = dL + alpha * sigma * Math.Sin(m);

            sigma += Math.Sin(m) * (alpha * sigma * Math.Sin(m));

            m = Math.Cos(u1) * Math.Cos(u2) * Math.Sin(lamda) / Math.Sin(sigma);
            m = Math.Atan(m / Math.Sqrt(1 - m * m));

            double A = Math.Atan(Math.Sin(lamda) / (Math.Cos(u1) * Math.Tan(u2) - Math.Sin(u1) * Math.Cos(lamda)));
            if (A <= 0) A += Math.PI;
            if (m <= 0) A += Math.PI;

            double M_ = Math.Atan(Math.Sin(u1) * Math.Tan(A) / Math.Sin(m));
            if (M_ <= 0) M_ += Math.PI;

            KK = ee * Math.Pow(Math.Cos(m), 2);
            alpha = ee / 2 + ee * ee / 8 + Math.Pow(ee, 3) / 16 - ee * (1 + ee) * KK / 16 + 3 * ee * KK * KK / 128;
            double beta = ee * (1 + ee) * KK / 16 - ee * KK * KK / 32;
            lamda = dL + Math.Sin(m) * (alpha * sigma + beta * Math.Sin(sigma) * Math.Cos(2 * M_ + sigma));

            A = Math.Atan(Math.Sin(lamda) / (Math.Cos(u1) * Math.Tan(u2) - Math.Sin(u1) * Math.Cos(lamda)));
            if (A <= 0) A += Math.PI;
            if (m <= 0) A += Math.PI;

            return A * 180 / Math.PI;
        }

        /// <summary>
        /// 由大地坐标计算大地线长
        /// </summary>
        /// <param name="B1">起点纬度</param>
        /// <param name="L1">起点经度</param>
        /// <param name="B2">末点纬度</param>
        /// <param name="L2">末点经度</param>
        /// <param name="a">参考椭球长半轴</param>
        /// <param name="f">参考椭球扁率倒数</param>
        /// <returns>两点间的大地线长</returns>
        public static double Bessel_BL_S(double B1, double L1, double B2, double L2, double a, double f)
        {
            if (L1 == L2) return Math.Abs(MeridianLength(B1, a, f) - MeridianLength(B2, a, f));

            double ee = (2 * f - 1) / f / f;       //第一偏心率的平方
            double ee2 = ee / (1 - ee);            //第二偏心率的平方

            double u1 = Math.Atan(Math.Sqrt(1 - ee) * Math.Tan(B1 * Math.PI / 180));
            double u2 = Math.Atan(Math.Sqrt(1 - ee) * Math.Tan(B2 * Math.PI / 180));
            double dL = (L2 - L1) * Math.PI / 180;

            double sigma = Math.Sin(u1) * Math.Sin(u2) + Math.Cos(u1) * Math.Cos(u2) * Math.Cos(dL);
            sigma = Math.Atan(Math.Sqrt(1 - sigma * sigma) / sigma);
            if (sigma <= 0) sigma += Math.PI;

            double m = Math.Cos(u1) * Math.Cos(u2) * Math.Sin(dL) / Math.Sin(sigma);
            m = Math.Atan(m / Math.Sqrt(1 - m * m));

            double KK = ee * Math.Pow(Math.Cos(m), 2);
            double alpha = ee / 2 + ee * ee / 8 + Math.Pow(ee, 3) / 16 - ee * (1 + ee) * KK / 16 + 3 * ee * KK * KK / 128;
            double lamda = dL + alpha * sigma * Math.Sin(m);

            sigma += Math.Sin(m) * (alpha * sigma * Math.Sin(m));

            m = Math.Cos(u1) * Math.Cos(u2) * Math.Sin(lamda) / Math.Sin(sigma);
            m = Math.Atan(m / Math.Sqrt(1 - m * m));

            double A = Math.Atan(Math.Sin(lamda) / (Math.Cos(u1) * Math.Tan(u2) - Math.Sin(u1) * Math.Cos(lamda)));
            if (A <= 0) A += Math.PI;
            if (m <= 0) A += Math.PI;

            double M_ = Math.Atan(Math.Sin(u1) * Math.Tan(A) / Math.Sin(m));
            if (M_ <= 0) M_ += Math.PI;

            KK = ee * Math.Pow(Math.Cos(m), 2);
            alpha = ee / 2 + ee * ee / 8 + Math.Pow(ee, 3) / 16 - ee * (1 + ee) * KK / 16 + 3 * ee * KK * KK / 128;
            double beta = ee * (1 + ee) * KK / 16 - ee * KK * KK / 32;

            lamda = dL + Math.Sin(m) * (alpha * sigma + beta * Math.Sin(sigma) * Math.Cos(2 * M_ + sigma));
            sigma = Math.Sin(u1) * Math.Sin(u2) + Math.Cos(u1) * Math.Cos(u2) * Math.Cos(lamda);
            sigma = Math.Atan(Math.Sqrt(1 - sigma * sigma) / sigma);
            if (sigma < 0) sigma += Math.PI;

            KK = ee2 * Math.Pow(Math.Cos(m), 2);
            alpha = Math.Sqrt(1 + ee2) / a * (1 - KK / 4 + 7 * KK * KK / 64 - 15 * Math.Pow(KK, 3) / 256);
            beta = KK / 4 - KK * KK / 8 + 37 * Math.Pow(KK, 3) / 512;
            double gama = KK * KK * (1 - KK) / 128;

            double S = (sigma - beta * Math.Sin(sigma) * Math.Cos(2 * M_ + sigma) - gama * Math.Sin(2 * gama) * Math.Cos(4 * M_ + 2 * sigma)) / alpha;

            return S;
        }