大地测量学高斯投影正反算

namespace 大地测量学2
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }

        private void button1_Click(object sender, EventArgs e)
        {
            double ρ = 206264.806247096355, 
                 a = 6378245, b = 6356863.0187730473, c = 6399698.9017827110, α = 1 / 298.3, e2 = 0.006693421622966, e3 = 0.006738525414683;//先BLA是角度，BmLmAm是弧度直接进行运算，可以直接进行三角函数运算，感觉在最上边声明MATH后就不用再一个一个的代用


            double B = 51.0 + 38.0 / 60.0 + 43.9023 / 3600.0, L = 111.0 + 2.0 / 60.0 + 13.136 / 3600.0;//角度
            B = B / 180 * Math.PI;
            L = L / 180 * Math.PI;//弧度

            double W = Math.Sqrt(1 - e2 * Math.Sin(B) * Math.Sin(B));
            double V = Math.Sqrt(1 +e3 * Math.Cos(B) * Math.Cos(B));

           // double N = a / W;

            double M = a * (1 - e2) / (W * W * W);
            double t = Math.Tan(B);
            double η = e3 * Math.Cos(B);

            B = 51.0 + 38.0 / 60.0 + 43.9023 / 3600.0;
            L = 111 + 2.0 / 60.0+ 13.136 / 3600;//重新化为角度 B= 51.012195083333332

            int n6,n3;



            n6 = (int)(L / 6) + 1;//投影带编号19
            n3 = (int)(L / 3 + 0.5);
            double L06 = 6 * n6 - 3;//L06即使这样存在，是否有需要强制转换 此时为111
            double l2 = L - L06;//此时是角度
            double zc = 2 * Math.PI * b + 4 * (a - b);//椭圆的周长40026876.244215891

            double X6=zc/60/2*(2*(n6-1)+1);//12341620.1752999

            B = B / 180 * Math.PI;
            L = L / 180 * Math.PI;//弧度
            l2 = l2 / 180 * Math.PI;//没卵用

            //double x = X6 + N / 2 / ρ / ρ * Math.Sin(B) * Math.Cos(B) * l2 * l2 + N / 24 / ρ / ρ / ρ / ρ * Math.Sin(B) * Math.Cos(B) * Math.Cos(B) * Math.Cos(B) * (5 - t * t + 9 * η * η + 4 * η * η * η * η)*l2*l2
            //  * l2 * l2 + N / 720 / ρ / ρ / ρ / ρ / ρ / ρ * Math.Sin(B) * Math.Cos(B) * Math.Cos(B) * Math.Cos(B) * Math.Cos(B) * Math.Cos(B) * (61 - 58 * t * t + t * t * t * t) * l2 * l2 * l2 * l2 * l2 * l2;
            //double y = N / ρ * Math.Cos(B) * l2 + N / 6 / ρ / ρ / ρ * Math.Cos(B) * Math.Cos(B) * Math.Cos(B) * (1 - t * t + η * η) * l2 * l2 * l2 + N / 120 / ρ / ρ / ρ / ρ / ρ * Math.Cos(B) * Math.Cos(B)
            //* Math.Cos(B) * Math.Cos(B) * Math.Cos(B) * (5 - 18 * t * t + t * t * t * t + 14 * η * η - 58 * η * η * t * t) * l2 * l2 * l2 * l2 * l2;
            double N=6399698.902-(21562.267-(108.973-0.62*Math.Cos(B) * Math.Cos(B)) * Math.Cos(B)* Math.Cos(B)) * Math.Cos(B)* Math.Cos(B) ;
            double B2=(51.0 + 38.0 / 60.0 + 43.9023 / 3600.0)*3600;//这个为秒的格式
            double a0 = (32140.404 - (135.3302 - (0.7092 - 0.004 * Math.Cos(B) * Math.Cos(B)) * Math.Cos(B) * Math.Cos(B)) * Math.Cos(B) * Math.Cos(B))/*a0*/;
            double a4 = (0.25 + 0.00252 * Math.Cos(B) * Math.Cos(B)) * Math.Cos(B) * Math.Cos(B) - 0.04166;
            double a6 = (0.166*Math.Cos(B) * Math.Cos(B) - 0.084) * Math.Cos(B) * Math.Cos(B);
            double a3 = (0.3333333 + 0.001123 * Math.Cos(B) * Math.Cos(B)) * Math.Cos(B) * Math.Cos(B) - 0.1666667;
            double a5 = 0.0083 - (0.1667 - (0.1968 + 0.004 * Math.Cos(B) * Math.Cos(B)) * Math.Cos(B) * Math.Cos(B)) * Math.Cos(B) * Math.Cos(B);
            double L0 = 6 * 19 - 3;//已知带号为19
            double l = L + L0;//角度
            l=l*3600/ρ;
            double x=6367558.4969*B2/ρ-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*Math.Sin (B)*Math .Cos (B);
            double y = (1 + (a3 + a5 * l * l) * l * l) * N * l * Math.Cos(B);

            Console.WriteLine("{0} {1} ", x, y);
            textBox1.Text = Convert.ToString(x);
            textBox2.Text = Convert.ToString(y);

        }

        private void button2_Click(object sender, EventArgs e)
        {
            textBox1.Text = "";
            textBox2.Text = "";

        }

        private void groupBox1_Enter(object sender, EventArgs e)
        {

        }

        private void button3_Click(object sender, EventArgs e)
        {
            double ρ = 206264.806247096355,
                a = 6378245, b = 6356863.0187730473, c = 6399698.9017827110, α = 1 / 298.3, e2 = 0.006693421622966, e3 = 0.006738525414683;//先BLA是角度，BmLmAm是弧度直接进行运算，可以直接进行三角函数运算，感觉在最上边声明MATH后就不用再一个一个的代用
         //计算子午面的弧长
            double Round = 2 * Math.PI * b + 4 * (a - b);
            double X=5724004.82211011;
            double Y=502559.918404555;
            //以度°为单位的bf
            double Bf = 27.11115372595 + 9.02468257083 * (X / 1000000 - 3) - 0.00579740442 * (X / 1000000 - 3) * (X / 1000000 - 3) - 0.00043532572 * (X / 1000000 - 3) * (X / 1000000 - 3) * (X / 1000000 - 3) + 0.00004857285 * (X / 1000000 - 3) * (X / 1000000 - 3) * (X / 1000000 - 3) * (X / 1000000 - 3);

            double W = Math.Sqrt(1 - e2 * Math.Sin(Bf) * Math.Sin(Bf));
            double V = Math.Sqrt(1 + e3 * Math.Cos(Bf) * Math.Cos(Bf));

            double N = a / W;

            double M = a * (1 - e2) / (W * W * W);
            double t = Math.Tan(Bf);
            double η = e3 * Math.Cos(Bf);

            double B = Bf - t / 2 / M / N * Y * Y + t / 24 / M / N / N / N * (5 + 3 * t * t + η * η - 9 * η * η * t * t) * Y * Y * Y * Y;
            Bf=Bf/180*Math .PI;
            double l = 1 / N / Math.Cos(Bf) * Y - 1 / 6 / N / N / N / Math.Cos(Bf) * (1 + 2 * t * t + η * η) * Y * Y * Y + 1 / 120 / N / N / N / N / N / Math.Cos(Bf) * (5 + 28 * t * t + 24 * t * t * t * t) * Y * Y * Y * Y * Y;


            double L0 = 6 * 19 - 3;//已知带号为19
            double y = Y - N - 500000;
            double L = L0 + l;//角度
            double D=(int)(L);
            double F=(int)((L-(int)(L))*60);
            double MM=(L-D-F/60)*3600;


            textBox3.Text = Convert.ToString(B);
            textBox4.Text = Convert.ToString(L);

            Console.WriteLine("{0} {1} {2}", D, F, MM);//111° 7分 35.6402392184452秒
        }

        private void button4_Click(object sender, EventArgs e)
        {
            textBox3.Text = "";
            textBox4.Text = "";
        }
    }
}

实验截图如上，在这里为了简单只用了一个固定数据，界面好不人性化。。。需要的读者可以自己改换为相应的数据，需要完整程序的可以联系我。
感谢阅读