大地经纬度和国家标准坐标转换

      平常的工作中就需要用到坐标转换的问题，特别是我做的webgis，由于需要在web页面上输入经纬度，而在实际上中采用国家标准的大地坐标来作图，所以就需要考虑坐标转换，否则人家输入到经纬度，想定位，结果坐标无法匹配，而最终无法实现功能，客户肯定会发牢骚的。这个是我开发来转换坐标的部分代码，结果很准确的，项目中应用过了的。

        /// <summary>
        /// 高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres)
        /// </summary>
        /// <param name="longitude">经度,以度为单位</param>
        /// <param name="latitude">纬度,以度为单位</param>
        /// <returns></returns>
        public Point3D GaussProjCal(double longitude, double latitude)
        {
            Point3D pnt = new Point3D();
            int ProjNo = 0; int ZoneWide; ////带宽
            double longitude1, latitude1, longitude0,  X0, Y0, xval, yval;
            double e2, ee, NN, T, C, A, M;
            ZoneWide = paraDegree;
            ProjNo = (int)(longitude / ZoneWide);
            longitude0 = originLongitude;
            longitude1 = Angle.DmsToRad(longitude);//经度转换为弧度
            latitude1 = Angle.DmsToRad(latitude);//纬度转换为弧度
            e2 = 2 * base.flat - flat * flat;
            ee = e2 * (1.0 - e2);
            double sinB = Math.Sin(latitude1);
            double tanB = Math.Tan(latitude1);
            double cosB = Math.Cos(latitude1);
            NN = longRadius / Math.Sqrt(1.0 - e2 * sinB * sinB);
            T = tanB * tanB;
            C = ee * cosB * cosB;
            A = (longitude1 - longitude0) * cosB;
            double a0 = 1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256;
            double a2 = 3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024;
            double a4 = 15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024;
            double a6 = 35 * e2 * e2 * e2 / 3072;
            M = longRadius * (a0 * latitude1 - a2 * Math.Sin(2 * latitude1) + a4 * Math.Sin(4 * latitude1) - a6 * Math.Sin(6 * latitude1));
            xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120);
            yval = M + NN * Math.Tan(latitude1) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24
            + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720);

            //X0 = 1000000L * ProjNo + 500000L;
            X0 = 500000L;
            Y0 = 0;
            xval = xval + X0;
            yval = yval + Y0;
            pnt.X = xval;
            pnt.Y = yval;

            return pnt;
        }

        /// <summary>
        /// 高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres)
        /// </summary>
        /// <param name="lola"></param>
        /// <returns></returns>
        public Point3D GaussProjCal(LoLatude lola)
        {
            double lo = lola.ToLongitude();
            double la = lola.ToLatitude();
            return GaussProjCal(lo, la);
        }

        /// <summary>
        /// 高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres)
        /// </summary>
        /// <param name="longitude">经度,以度为单位</param>
        /// <param name="latitude">纬度,以度为单位</param>
        /// <param name="x">X坐标</param>
        /// <param name="y">Y坐标</param>
        public void GaussProjCal(double longitude, double latitude,out double x,out double y)
        {
            Point3D pnt = GaussProjCal(longitude, latitude);
            x = pnt.X;
            y = pnt.Y;
        }

        /// <summary>
        /// 高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)
        /// </summary>
        /// <param name="X">X坐标</param>
        /// <param name="Y">Y坐标</param>
        /// <param name="longitude">经度</param>
        /// <param name="latitude">纬度</param>
        public void GaussProjInvCal(double X, double Y, out double longitude, out double latitude)
        {
            int ProjNo; int ZoneWide; ////带宽
            double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
            double e1, e2,ee, NN, T, C, M, D, R, u, fai;
            ZoneWide = paraDegree; ////6度带宽
            ProjNo = (int)(X / 1000000L); //查找带号
            longitude0 = originLongitude; //中央经线
            //X0 = ProjNo * 1000000L + 500000L;
            X0 = 500000L;
            Y0 = 0;
            xval = X - X0;
            yval = Y - Y0; //带内大地坐标
            e2 = 2 * flat - flat * flat;
            double sqrtE2 = Math.Sqrt(1 - e2);
            e1 = (1.0 - sqrtE2) / (1.0 + sqrtE2);
            ee = e2 / (1 - e2);
            M = yval;
            double q0 = 1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256;
            double q2 = 3 * e1 / 2 - 27 * e1 * e1 * e1 / 32;
            double q4 = 21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32;
            double q6 = 151 * e1 * e1 * e1 / 96;
            double q8 = 1097 * e1 * e1 * e1 * e1 / 512;
            u = M / (longRadius * q0);
            fai = u + q2 * Math.Sin(2 * u) + q4 * Math.Sin(4 * u)
            + q6 * Math.Sin(6 * u) + q8 * Math.Sin(8 * u);
            double cosB = Math.Cos(fai);
            double tanB = Math.Tan(fai);
            double sinB = Math.Sin(fai);
            C = ee * cosB * cosB;
            T = tanB * tanB;
            NN = longRadius / Math.Sqrt(1.0 - e2 * sinB * sinB);
            R = longRadius * (1 - e2) / Math.Sqrt((1 - e2 * sinB * sinB) * (1 - e2 * sinB * sinB) * (1 - e2 * sinB * sinB));
            D = xval / NN;
            //计算经度(Longitude) 纬度(Latitude)
            longitude1 = longitude0 + (D - (1 + 2 * T + C) * D * D * D / 6 + (5 - 2 * C + 28 * T - 3 * C * C + 8 * ee + 24 * T * T) * D
            * D * D * D * D / 120) / cosB;
            latitude1 = fai - (NN * tanB / R) * (D * D / 2 - (5 + 3 * T + 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24
            + (61 + 90 * T + 298 * C + 45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720);
            //转换为度 DD
            longitude = longitude1 / Pi2;
            latitude = latitude1 / Pi2;
        }

        /// <summary>
        /// 大地坐标反算经纬度坐标
        /// </summary>
        /// <param name="x">X坐标</param>
        /// <param name="y">Y坐标</param>
        /// <returns></returns>
        public LoLatude GaussProjInvCal(double x,double y)
        {
            double lo, la;
            GaussProjInvCal(x, y, out lo, out la);
            lo = lo * Pi2;
            la = la * Pi2;
            LoLatude lola = new LoLatude();
            lola.Longitude = Angle.RadToDms(lo);
            lola.Latitude = Angle.RadToDms(la);
            return lola;
        }

        /// <summary>
        /// 大地坐标反算经纬度坐标
        /// </summary>
        /// <param name="pnt">坐标点</param>
        /// <returns></returns>
        public LoLatude GaussProjInvCal(Point3D pnt)
        {
            double lo, la;
            GaussProjInvCal(pnt.X, pnt.Y, out lo, out la);
            lo = lo * Pi2;
            la = la * Pi2;
            LoLatude lola = new LoLatude();
            lola.Longitude = Angle.RadToDms(lo);
            lola.Latitude = Angle.RadToDms(la);
            return lola;
        }