测绘程序设计-参心大地坐标转换为参心空间直角坐标

用C#写的关于坐标的转换，由BLH转换为XYZ。

在项目中添加类（class）的方法来写的程序，然后运行程序。

废话不多说，直接上代码：

internal class CoorTrans
    {
        private const double a = 6378137;
        private const double f = 1 / 298.257222101;

        public void BLH2XYZ(double B, double L, double H, 
            out double X, out double Y, out double Z)
        {
            double e2 = ComputeE2();
            double N = ComputeN(B);

            double b = Deg2Rad(B);
            double l = Deg2Rad(L);
            X = (N + H) * Math.Cos(b) * Math.Cos(l);
            Y = (N + H) * Math.Cos(b) * Math.Sin(l);
            Z = (N * (1 - e2) + H) * Math.Sin(b);
        }
        public void XYZ2BLH(double X,double Y,double Z,
            out double B,out double L,out double H)
        {
            L= ComputeL(X,Y);
            double B0 = 0;
            double dateB = 90;
            do
            {
                H = ComputeH(X, Y, B0);
                B = ComputeB(B0, X, Y, Z, H);
                dateB = B - B0;
                B0 = B;
            }
            while (Math.Abs(dateB) < 1e-10);
        }
        //B is not 90 deg. 
        private double ComputeH(double X,double Y,double B)
        {
            double N = ComputeN(B);
            double H = Math.Sqrt(X * X + Y * Y) / Math.Cos(Deg2Rad(B)) - N;
            return H;
        }
        private double ComputeB(double B0,double X,double Y,double Z,double H)
        {
            double N = ComputeN(B0);
            double e2 = ComputeE2();
            double upper = Z * (N + H);
            double under = Math.Sqrt(X * X + Y * Y) * (N * (1 - e2) + H);
            double B = Math.Atan(upper / under);
            B = Rad2Deg(B);
            return B;
        }
        private double ComputeL(double X,double Y)
        {
            double L = 0;
            if (Math.Abs(X) < 1e-15 && Y > 0)
                L = 0;
            else if (Math.Abs(X) < 1e-15 && Y < 0)
                L = 180;
            else
            {
                L = Math.Atan(Y / X);
                L = Rad2Deg(L);
                if (Y < 0 && X > 0)
                    L = L + 360;
                if (Y > 0 && X < 0)
                    L = 180 + L;
                if (Y < 0 && X < 0)
                    L = 180 + L;
                
            }
            return L;
        }
        private double ComputeE2()
        {
            double e2 = 2 * f - f * f;
            return e2;
        }
        private double ComputeN(double B)
        {
            double e2 = ComputeE2();
            double b = Deg2Rad(B);
            double sinb = Math.Sin(b);
            double W = Math.Sqrt(1 - e2 * sinb *sinb);
            double N = a / W;
            return N;
        }
        private double Deg2Rad(double deg)
        {
            double red = deg * Math.PI / 180.0;
            return red;
        }
        private double Rad2Deg(double rad)
        {
            double deg = rad * Math.PI;
            return deg;
        }
    }

主代码如下：

internal class Program
    {
        static void Main(string[] args)
        {
            double B = 30.30, L = 114.45, H = 50;
            CoorTrans trans = new CoorTrans();
            double X, Y, Z;
            trans.BLH2XYZ(B, L, H, out X, out Y, out Z);
            //trans.XYZ2BLH(X, Y, Z, out B, out L, out H);
            Console.WriteLine("{0} {1} {2}",X,Y,Z);
            Console.ReadKey();
        }
    }