C#白塞尔大地正算

白塞尔大地正算是地理测量中的重要计算方法，用于计算两个地理坐标之间的距离和方位角。在地图应用、导航系统和测绘工程中广泛应用。

椭球面点的大地经度 L、大地纬度 B，两点间的大地线长度 S 及其正、反大地方位角 A1、 A2 ,统称为大地元素,如图 1 所示。如果知道某些大地元素推求另外一些大地元素，这样的计算就叫做大地主题解算。

已知：大地线起点P1 的纬度B1，经度L1，大地方位角A1，起点P1 到终点P2的大地线长度 S；

计算：大地线终点 P2 的纬度B2，经度L2及大地方位角 A2。

namespace GeoDir
{
    public class BesselDirect
    {
        private Ellipsoid Ell;

        public BesselDirect(Ellipsoid Ell)
        {
            this.Ell = Ell;
        }

        private void CalReduceLat(double B1,ref double sinu1,ref double cosu1)
        {
            double e1 = Ell.e1;
            double W1 = GeodesyCal.GetW(e1, B1);
            sinu1 = Math.Sqrt(1 - e1 * e1) / W1 * Math.Sin(B1);
            cosu1 = Math.Cos(B1) / W1;
        }

        private void CalABC_AlphaBeta(double sinA0,double[] ABC,ref double alpha,ref double beta,ref double gama)
        {
            double cos2_A0 = 1 - sinA0 * sinA0;
            double k_2 = 0;
            double e1 = Ell.e1;
            double e2 = Ell.e2;
            double b = Ell.b;
            k_2 = GeodesyCal.Getk_2(e2, cos2_A0);
            GeodesyCal.GetABC(b, k_2, ABC);
            alpha = GeodesyCal.GetAlpha(e1, cos2_A0);
            beta = GeodesyCal.GetBeta(e1, cos2_A0);
            gama = GeodesyCal.GetGama(e1, cos2_A0);
        }

        private void CalGeodesicLength(double sigma1,double[] ABC,double S,ref double sigma)
        {
            double _sigma = 0;
            do
            {
                _sigma = sigma;
                sigma = ABC[0] * S + ABC[1] * Math.Sin(sigma) * Math.Cos(2 * sigma1 + sigma)
                    + ABC[2] * Math.Sin(2 * sigma) * Math.Cos(4 * sigma1 + 2 * sigma);
            } while (Math.Abs(_sigma - sigma) > 0.000000001);
        }

        public void DirectSolution(GeodesicInfo geodesic)
        {
            double B1 = GeodesyCal.DMS2RAD(geodesic.P1.B);
            double L1 = GeodesyCal.DMS2RAD(geodesic.P1.L);
            double A1 = GeodesyCal.DMS2RAD(geodesic.A1);
            //double A1 = geodesic.A1;
            double S = geodesic.S;
            //椭球参数
            double e1, e2, b;
            e1 = Ell.e1;
            e2 = Ell.e2;
            b = Ell.b;
            //计算归化纬度
            double sinu1 = 0, cosu1 = 0;
            CalReduceLat(B1, ref sinu1, ref cosu1);
            //计算辅助函数
            double sinA0 = cosu1 * Math.Sin(A1);
            double cot_sigma1 = cosu1 * Math.Cos(A1) / sinu1;
            double sigma1 = Math.Atan(1.0 / cot_sigma1);
            //计算ABC及α与β
            double[] ABC = new double[3];
            double alpha = 0;
            double beta = 0;
            double gama = 0;
            CalABC_AlphaBeta(sinA0, ABC, ref alpha, ref beta, ref gama);

            //计算球面长度           
            double sigma = ABC[0] * S;
            CalGeodesicLength(sigma1, ABC, S, ref sigma);

            //计算经差改正数
            double lamda_L = sinA0 * (alpha * sigma + beta * Math.Sin(sigma) * Math.Cos(2 * sigma1 + sigma)
                + gama * Math.Sin(2 * sigma) * Math.Cos(4 * sigma1 + 2 * sigma));

            //计算终点大地坐标及大地方位角
            double sinu2 = sinu1 * Math.Cos(sigma) + cosu1 * Math.Cos(A1) * Math.Sin(sigma);
            double B2 = Math.Atan(1 / Math.Sqrt(1 - e1 * e1) * sinu2 / Math.Sqrt(1 - sinu2 * sinu2));
            double lamba = Math.Atan(Math.Sin(sigma) * Math.Sin(A1) / (cosu1 * Math.Cos(sigma) - sinu1 * Math.Sin(sigma)
                * Math.Cos(A1)));
            double sinA1 = Math.Sin(A1);
            lamba = GeodesyCal.DirJudgeLamba(sinA1, lamba);

            double L2 = L1 + lamba - lamda_L;

            double A2 = Math.Atan(cosu1 * Math.Sin(A1) / (cosu1 * Math.Cos(sigma)
                * Math.Cos(A1) - sinu1 * Math.Sin(sigma)));
            //double sinA1 = Math.Sin(A1);
            A2 = GeodesyCal.DirJudgeA2(sinA1, A2);
            if (A2 > 2 * Math.PI) A2 -= 2 * Math.PI;
            if (A2 < 0) A2 += 2 * Math.PI;
            if (A1 >= Math.PI && A2 >= Math.PI) A2 = A2 - Math.PI;
            if (A1 < Math.PI && A2 < Math.PI) A2 = A2 + Math.PI;

            //角度转换
            geodesic.P2.B = GeodesyCal.RAD2DMS(B2);
            geodesic.P2.L = GeodesyCal.RAD2DMS(L2);
            geodesic.A2 = GeodesyCal.RAD2DMS(A2);
        }

        public void CalPro(List<GeodesicInfo> geodesics)
        {
            for(int i = 0; i < geodesics.Count; i++)
            {
                DirectSolution(geodesics[i]);
            }
        }
        
    }
}

主要计算部分实现代码，这是其中计算类的部分。

具体代码实现以及程序开发文档可以在我的资源中获取。

(44条消息) C#白塞尔大地主题正算资源-CSDN文库