利用C#编写一个高斯正反算程序

一、代码界面展示

整个界面控件为tabControl，groupBox，label，textbox，comboBox，button，richTextBook。

二、代码运算结果展示

数据结果采用国家统一坐标（横坐标上加上500000m，在坐标前面冠以带号）。
内置了4种椭球参数（CGCS200,WGS84,西安80，克拉索夫斯基椭球体）。
可以计算三度带和六度带。
计算曲率半径功能是我用来做大地测量学作业用的，顺便写在这里的。

2.1 数据相互转换验证

相互转换的结果肯定是一样的（我自己测试了10几次不同的数据，是没有问题的）。

2.2 数据结果正确性验证

2.2.1 高斯正算

数据结果的正确性需要使用其他的软件进行验证。
这里我选取的是https://blog.csdn.net/starmings/article/details/123034189计算的数据进行验证。

软件结果如下：
采用西安80椭球参数，高斯正算的结果与软件的结果是一致的。

2.2.2 高斯反算

软件结果如下：

采用西安80椭球参数，高斯反算的结果与软件的结果是一致的。（我只保留了0.01，已经可以了。）

三、源代码

3.1 变量设置

        //大地经度 度分秒public double L;public double Ld;public double Lf;public double Lm;//大地纬度 度分秒public double B;public double Bd;public double Bf;public double Bm;double L0;//中央子午线度数double N1;//带号double a;//长半轴double b;//短半轴double e1;//第一偏心率平方double e2; //第二偏心率平方double m0,m2,m4,m6,m8;double a0, a2, a4, a6, a8;int i ;

3.2 高斯正算

             //m0-m8的值m0 = a * (1 - e1);m2 = 1.5 * m0 * e1;m4 = 1.25 * m2 * e1;m6 = (7.0 / 6) * m4 * e1;m8 = (9.0 / 8) * m6 * e1;//a0-a8的值a0 = m0 + m2 / 2.0 + (3.0 / 8) * m4 + (5.0 / 16) * m6 + (35.0 / 128) * m8;a2 = (1.0 / 2) * (m2 + m4) + (15.0 / 32) * m6 + (7.0 / 16) * m8;a4 = (1.0 / 8) * m4 + (3.0 / 16) * m6 + (7.0 / 32) * m8;a6 = (1.0 / 32) * m6 + (1.0 / 16) * m8;a8 = (1.0 / 128) * m8;//将大地经度转化为度数Ld = double.Parse(textBox1.Text);Lf = double.Parse(textBox2.Text);Lm = double.Parse(textBox3.Text);L = Ld + Lf / 60 + Lm / 3600;//将大地纬度转化为度数Bd = double.Parse(textBox4.Text);Bf = double.Parse(textBox5.Text);Bm = double.Parse(textBox6.Text);B = Bd + Bf / 60 + Bm / 3600;if (comboBox2.Text == "三度带"){N1 = Math.Truncate((L + 1.5) / 3);//带号L0 = 3 * N1;}if (comboBox2.Text == "六度带"){N1 = Math.Truncate(L / 6) + 1;L0 = 6 * N1 - 3;}//求椭球面上某点到赤道的子午线弧长，精度至0.001mdouble radB = (Math.PI / 180) * B; //将纬度转换为弧度double Δ1, Δ2, Δ3, Δ4, Δ5;Δ1 = a0 * radB;Δ2 = (a2 / 2) * Math.Sin(2 * radB);Δ3 = (a4 / 4) * Math.Sin(4 * radB);Δ4 = (a6 / 6) * Math.Sin(6 * radB);Δ5 = (a8 / 8) * Math.Sin(8 * radB);double X =Math.Round(Δ1 - Δ2 + Δ3 - Δ4 + Δ5, 8);//卯酉圈曲率半径Ndouble N = a / Math.Sqrt(1 - e1 * Math.Sin(radB) * Math.Sin(radB));// tan(B)的值double t = Math.Tan(radB);// η的值double η =Math.Sqrt(e2) * Math.Cos(radB) ;//l的值double ρm = 180 * 60 * 60 / Math.PI;double l = ((L - L0) * 3600) / ρm;//高斯坐标中的x值和y值double x = Math.Round( X + (N / 2) * t * Math.Cos(radB) * Math.Cos(radB) * l * l + (N / 24) * t * (5 - t * t + 9 * η* η + 4 * Math.Pow(η,4)) * Math.Pow(Math.Cos(radB), 4) * Math.Pow(l, 4) + (N / 720) * t * (61 - 58 * t * t + Math.Pow(t, 4)) * Math.Pow(l, 6) * (Math.Pow(Math.Cos(radB), 6)),4);double y = Math.Round( N * Math.Cos(radB) * l + (N / 6.0) * (1 - t * t + η* η) * Math.Pow(Math.Cos(radB), 3) * Math.Pow(l, 3) + (N / 120.0) * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * η * η - 58 * η * η * t * t) * Math.Pow(Math.Cos(radB), 5) * Math.Pow(l, 5),4);richTextBox1.AppendText(Environment.NewLine);richTextBox1.Text = "大地经度为L:" + Ld + "°" + Lf + "′" + Lm + "″" + "        " + "大地经度为B:" + Bd + "°" + Bf + "′" + Bm+"″" + "\r\n";richTextBox1.Text += "大地坐标转高斯平面坐标后：" + "\r\n";richTextBox1.Text += "高斯X坐标：" + x + "米" + "        " + "高斯Y坐标：" + +N1 + "" + (y+500000)+  "米" + "\r\n";

3.3 高斯反算

 //输入坐标x,ydouble x = double.Parse(textBox7.Text);double y1 = double.Parse(textBox8.Text);double y = y1-Math.Truncate(y1 / Math.Pow(10, 6)) * Math.Pow(10, 6)-500000;double[] Bf1 = new double[20];double[] F = new double[20];Bf1[0] = x / a0;//底点纬度F[0] = (-a2 / 2.0) * Math.Sin(2 * Bf1[0]) + a4 / 4.0 * Math.Sin(4 * Bf1[0]) - a6 / 6.0 * Math.Sin(6 * Bf1[0]);Bf1[1] = (x - F[0]) / a0;double Δd =Bf1[1] - Bf1[0];for ( i = 0;i<10 ; i++){F[i] = (-a2 / 2.0) * Math.Sin(2 * Bf1[i]) + a4 / 4.0 * Math.Sin(4 * Bf1[i]) - a6 / 6.0 * Math.Sin(6 * Bf1[i]);Bf1[i + 1] = (x - F[i]) / a0;Δd = Bf1[i + 1] - Bf1[i];if (Δd < Math.Pow(10, -10)){i = i + 1;break;}}// tan(Bf)的值double tf = Math.Tan(Bf1[i]);// ηf的值double ηf = e2 * Math.Cos(Bf1[i]) * Math.Cos(Bf1[i]);//计算Mf的值double Mf = a * (1 - e1) / Math.Pow(1 - e1 * Math.Sin(Bf1[i]) * Math.Sin(Bf1[i]), 1.5);//计算Nf的值double Nf = a / Math.Sqrt(1 - e1 * Math.Sin(Bf1[i]) * Math.Sin(Bf1[i]));//计算y坐标的前2位int aa = Convert.ToInt32(Math.Truncate(y1 / Math.Pow(10, 6)));if (comboBox2.Text == "三度带"){L0 = 3 * aa;}if (comboBox2.Text == "六度带"){L0 = 6 * aa - 3;}//求大地经度LL = (Math.PI / 180) * L0 + y / (Nf * Math.Cos(Bf1[i])) - (1 + 2 * tf * tf + ηf ) * Math.Pow(y, 3) / (6 * Math.Pow(Nf, 3) * Math.Cos(Bf1[i])) + (5 + 28 * tf * tf + 24 * Math.Pow(tf, 4) + 6 * ηf + 8 * ηf  * Math.Pow(tf,4)) * Math.Pow(y, 5) / (120 * Math.Pow(Nf, 5) * Math.Cos(Bf1[i]));//转换经度double degrees = (180 / Math.PI) * L;//将弧度转化为度int d1 = Convert.ToInt32(Math.Truncate(degrees));//度int f1 = Convert.ToInt32(Math.Truncate((degrees - d1) * 60));//分double m1 = Math.Round(((degrees - d1) * 60 - f1) * 60,2);//秒if (m1 == 60){m1 = 0;f1 += 1;}//求大地纬度BB = Bf1[i] - (tf * y * y) / (2 * Mf * Nf)+ tf * (5 + 3 * tf * tf + ηf  - 9 * ηf  * tf * tf) * Math.Pow(y, 4) / (24 * Mf * Math.Pow(Nf, 3)) + tf * (61 + 90 * tf * tf + 45 * Math.Pow(tf, 4)) * Math.Pow(y, 6) / (720 * Mf * Math.Pow(Nf, 5));double degrees1 = (180 / Math.PI) * B;//将弧度转化为度int d = Convert.ToInt32(Math.Truncate(degrees1));//度int f = Convert.ToInt32(Math.Truncate((degrees1 - d) * 60));//分double m = Math.Round(((degrees1 - d) * 60 - f) * 60,2);//秒if (m == 60){m = 0;f += 1;}richTextBox1.AppendText(Environment.NewLine);richTextBox1.Text += "高斯平面坐标转大地坐标后：" + "\r\n";richTextBox1.Text +=  "大地坐标L：" + d1+"°   "+f1+"′  "+m1+"″   "+ "大地坐标B：" + d + "°   " + f + "′  " + m + "″   "+"\r\n";

四、结语

最开始写这个高斯正反算程序的时候，公式完全按照大地测量学实习指导书上的进行编写，写完过后进行数据验证，发现相互转换就不能实现，分位就出现了偏差，说明正算公式或者反算公式出现了问题，经过检查，正算公式修改成了与https://blog.csdn.net/du12300/article/details/109307386所展示的公式。反算公式是对的，主要底点纬度的精度，我没在CSDN上找到相关的高斯反算文章，需要你自己发力了。
这个程序没有很困难的地方，就是公式的编写，可能就是需要注意double的精度丢失对数据结果带来的影响。
如果你需要批量处理，其实也简单，把我那个水准或者导线的程序里面的批量导入部分copy过来，设置几个变量，再写几个for语句就完事。
我的其他文章：

利用Python编写一个高斯正反算程序：https://blog.csdn.net/Zj1638/article/details/125740379
利用C#编写一个附和闭合导线简易平差程序：https://blog.csdn.net/Zj1638/article/details/125639541
利用C#编写一个水准测量近似平差程序：https://blog.csdn.net/Zj1638/article/details/119303957
利用C#编写一个GPS高程拟合（二次曲面拟合模型）程序：https://blog.csdn.net/Zj1638/article/details/125752243