C#编程练习（02）：大地坐标系（LBH）向空间直角坐标系（XYZ）的转换及其逆转换

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace XYZ2BLH
{
    class Program
    {
        //输出为度格式的字符串 （度分秒.秒的小数部分）
        //例如：35度23分32.9秒    表述为352332.9
        static string[] XYZ2LBH(double x,double y,double z)   
        {
            double[] blh = {0.0, 0.0, 0.0}; 
            double a = 6378137;
            double b = 6356752.3142;
            double e2 = 0.0066943799013;
            double ePie2 = 0.00673949674227;
            double xita = Math.Atan(z * a / (Math.Sqrt(x * x + y * y) * b));
            double xitaSin3 = Math.Sin(xita) * Math.Sin(xita) * Math.Sin(xita);
            double xitaCos3 = Math.Cos(xita) * Math.Cos(xita) * Math.Cos(xita);
            double B = Math.Atan((z + ePie2 * b * xitaSin3) / ((Math.Sqrt(x * x + y * y) - e2 * a *  xitaCos3)));
            double L = Math.Atan(y / x);
            double N = a / Math.Sqrt(1 - e2 * Math.Sin(B) * Math.Sin(B)); ;
            double H = Math.Sqrt(x * x + y * y) / Math.Cos(B) - N;
            //转化为角度
            B = 180 * B / Math.PI;
            if (B < 0.0) { B = B + 90; }                                 
            L = 180 * L / Math.PI;
            if (L < 0.0) { L = L + 180; }
            return new string[] {rad2dms(B), rad2dms(L), H.ToString()};  
        }
        //获取一个小数的整数部分和小数部分
        static double[] modf(double t)
        {
            t += 1.0e-14;//避免角度转换过程产生无效近似
            double t_integer = Math.Floor(t);
            return new double[] { t_integer, t - t_integer };      
        }
        //角度转弧度
        static double turn1(double t)
        {
            double[] t1 = modf(t);
            double[] t2 = modf(t1[1] * 100);
            double x1 = t1[0];
            double x2 = t2[0];
            double x3 = t2[1] * 100;        
            //Console.Write("{0:N} \t {1:N}\t {2:N}\r\n", x1, x2, x3);
            return (x1 + x2 / 60 + x3 / 3600) / 180 * Math.PI;
        }
        //弧度转角度
        static string rad2dms(double rad)
        {
            int du, fen;
            double miao;
            du = (int)rad;
            fen = (int)((rad - (double)du) * 60);
            miao = ((rad - (double)du) * 60 - (double)fen) * 60;
            return du.ToString() + fen.ToString().PadLeft(2,'0') + miao.ToString();        
        }
        //经度、纬度和大地高转向空间直角坐标
        static double[] LBH2XYZ(double[] lbh)
        {
            double B = turn1(lbh[1]);
            double L = turn1(lbh[0]);
            double H = lbh[2];
            double A = 1 / 298.257223563;
            double a = 6378137.0000000000;
            double E = 2 * A - A * A;
            double W = Math.Sqrt(1 - E * Math.Sin(B) * Math.Sin(B));
            double N = a / W;
            double[] xyz = { 0.0, 0.0, 0.0 };
            xyz[0] = (N + H) * Math.Cos(B) * Math.Cos(L);
            xyz[1] = (N + H) * Math.Cos(B) * Math.Sin(L);
            xyz[2] = (N * (1 - E) + H) * Math.Sin(B);
            return xyz;
        }
        static void Main(string[] args)
        {
            double[] lbh = { 108.5743, 32.2352, 100.59 };
            double[] xyz = LBH2XYZ(lbh);
            Console.WriteLine("{0:N6} {1:N6} {2:N6}", xyz[0], xyz[1], xyz[2]);
            string[] lbhTmp = XYZ2LBH(xyz[0], xyz[1], xyz[2]);
            Console.WriteLine("{0:N6} {1:N6} {2:N6}", double.Parse(lbhTmp[0]), double.Parse(lbhTmp[1]), double.Parse(lbhTmp[2]));         
        }
    }
}