C# 高斯正算 高斯反算 高斯换带等

首先，你要确定椭球参数：

a = 6378140;  //西安80椭球 IGA75
e2 = 0.006694384999588;
m0 = a * (1 - e2);
m2 = 3.0 / 2 * e2 * m0;
m4 = 5.0 / 4 * e2 * m2;
m6 = 7.0 / 6 * e2 * m4;
m8 = 9.0 / 8 * e2 * m6;
a0 = m0 + m2 / 2 + (3.0 / 8.0) * m4 + (5.0 / 16.0) * m6 + (35.0 / 128.0) * m8;
a2 = m2 / 2 + m4 / 2 + 15.0 / 32 * m6 + 7.0 / 16 * m8;
a4 = m4 / 8 + 3.0 / 16 * m6 + 7.0 / 32 * m8;
a6 = m6 / 32 + m8 / 16;
a8 = m8 / 128;
xx = 0;
yy = 0;
_x = 0;
_y = 0;
BB = 0;
LL = 0;

 下面才开始正题：

 高斯正算：把经纬度坐标转换为平面坐标

 

        void GaussPositive(double B, double L, double L0)
        {
            double X, t, N, h2, l, m, Bmiao, Lmiao;
            int Bdu, Bfen, Ldu, Lfen;
            Bdu = (int)B;
            Bfen = (int)(B * 100) % 100;
            Bmiao = (B - Bdu - Bfen * 0.01) * 10000.0;
            B = Bdu * PI / 180.0 + (Bfen / 60.0) * PI / 180.0 + Bmiao / 3600.0 * PI / 180.0;
            Ldu = (int)L;
            Lfen = (int)(L * 100) % 100;
            Lmiao = (L - Ldu - Lfen * 0.01) * 10000.0;
            L = Ldu * PI / 180.0 + (Lfen / 60.0) * PI / 180 + Lmiao / 3600.0 * PI / 180.0;
            l = L - L0 * PI / 180;
            X = a0 * B - Math.Sin(B) * Math.Cos(B) * ((a2 - a4 + a6) + (2 * a4 - 16.0 / 3.0 * a6) * Math.Sin(B) * Math.Sin(B) + 16.0 / 3.0 * a6 * Math.Pow(Math.Sin(B), 4)) + a8 / 8.0 * Math.Sin(8 * B);
            t = Math.Tan(B);
            h2 = e2 / (1 - e2) * Math.Cos(B) * Math.Cos(B);
            N = a / Math.Sqrt(1 - e2 * Math.Sin(B) * Math.Sin(B));
            m = Math.Cos(B) * l;
            xx = X + N * t * ((0.5 + (1.0 / 24.0 * (5 - t * t + 9 * h2 + 4 * h2 * h2) + 1.0 / 720.0 * (61 - 58 * t * t + Math.Pow(t, 4)) * m * m) * m * m) * m * m);
            yy = N * ((1 + (1.0 / 6.0 * (1 - t * t + h2) + 1.0 / 120.0 * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * h2 - 58 * h2 * t * t) * m * m) * m * m) * m);
            yy = yy + 500000;
        }

    高斯反算：把平面坐标转换为经纬度坐标：

   

         void GaussNegative(double x, double y, double L0)
        {
            double Bf, Vf, l, tf, hf2, Nf, Bmiao, Lmiao;
            int Bdu, Bfen, Ldu, Lfen;
            y = y - 500000;
            Bf = hcfansuan(x);
            Vf = Math.Sqrt(1 + e2 / (1 - e2) * Math.Cos(Bf) * Math.Cos(Bf));
            tf = Math.Tan(Bf);
            hf2 = e2 / (1 - e2) * Math.Cos(Bf) * Math.Cos(Bf);
            Nf = a / Math.Sqrt(1 - e2 * Math.Sin(Bf) * Math.Sin(Bf));
            BB = (Bf - 0.5 * Vf * Vf * tf * (Math.Pow(y / Nf, 2) - 1.0 / 12 * (5 + 3 * tf * tf + hf2 - 9 * hf2 * tf * tf) * Math.Pow(y / Nf, 4) + 1.0 / 360 * (61 + 90 * tf * tf + 45 * tf * tf) * Math.Pow(y / Nf, 6))) * 180 / PI;
            Bdu = (int)BB;
            Bfen = (int)((BB - Bdu) * 60);
            Bmiao = ((BB - Bdu) * 60 - Bfen) * 60;
            BB = Bdu + 0.01 * Bfen + 0.0001 * Bmiao;
            l = 1.0 / Math.Cos(Bf) * (y / Nf - 1.0 / 6.0 * (1 + 2 * tf * tf + hf2) * Math.Pow(y / Nf, 3) + 1.0 / 120.0 * (5 + 28 * tf * tf + 24 * Math.Pow(tf, 4) + 6 * hf2 + 8 * hf2 * tf * tf) * Math.Pow(y / Nf, 5)) * 180.0 / PI;
            LL = L0 + l;
            Ldu = (int)LL;
            Lfen = (int)((LL - Ldu) * 60);
            Lmiao = ((LL - Ldu) * 60 - Lfen) * 60;
            LL = Ldu + 0.01 * Lfen + 0.0001 * Lmiao;
        }

   里面涉及到的弧长反算：

  

        double hcfansuan(double pX)
        {
            double Bf0 = pX / a0;
            double Bf1, Bf2;
            Bf1 = Bf0;
            Bf2 = (pX - F(Bf1)) / a0;
            while ((Bf2 - Bf1) > 1.0E-11)
            {
                Bf1 = Bf2;
                Bf2 = (pX - F(Bf1)) / a0;
            }
            return Bf1;
        }

    高斯换带就比较简单了：

   

        void GaussZone(double x, double y, double L0, double L0new)
        {
            GaussNegative(x, y, L0);
            GaussPositive(BB, LL, L0new);
        }