高斯投影正反算c代码

1、高斯投影正反算程序设计一. 程序设计流程本程序的设计思路如下：(1) ,程序采用VS08版本作为开发平台，并采用C雄言作为开发语言，设计为WindowsForm 窗体程序形式。(2) ,程序主要的算法来自于教材。但是本程序为了更加实用，添加了更多的解算基准，包括：WGS-84国际椭球1975,克氏椭球，和 2000国家大地坐标系。(3) ,程序为了更方便的读取数据和输出数据，故需要自己定义了固定的数据输入格式和 数据输出格式或形式，请老师注意查看。二. 代码using System;using System.Collections.Generic;using System.ComponentM

2、odel;using System.Data;using System.Drawing;using System.Text;using System.Windows.Forms;namespace Gausspublic partial class Form1 : Form/大地坐标/Geodetic Coordinatepublic struct CRDGEODETICpublic double dLongitude;public double dLatitude;public double dHeight;/笛卡尔坐标/Cartesian Coordinatepublic struct C

3、RDCARTESIAN(public double x;public double y;public double z;)public Form1()(InitializeComponent();)private void button1_Click(object sender, EventArgs e)(double ee = 0;double a = 0;string tt;try(tt= boBox1.ItemsboBox1.SelectedIndex.ToString();)catch(MessageBox.Show("Gauss Inverse: Choose datum

4、error!");return;)if (tt.CompareTo(克氏椭球")=0)a = 6378245.00;ee = Math.Sqrt(0.006693421622);)if (tt.CompareTo("WGS-84") = 0)(a = 6378137.00;ee = Math.Sqrt(0.00669437999013);)if (tt.CompareTo("1975国际椭球")=0)(a = 6378140.00;ee = Math.Sqrt(0.006694384999588);)if (tt.CompareTo(

5、"2000国家大地坐标系 ")=0)(a = 6378137.0;ee =Math.Sqrt(0.0066943802290);)const double pai = 3.1415926;double b = Math.Sqrt(a * a * (1 - ee * ee);double c = a * a / b;double epp = Math.Sqrt(a * a - b * b) / b / b);CRDGEODETIC pcrdGeo;CRDCARTESIAN pcrdCar;double midlong;/ 求纬度string temp;temp = textB

6、ox1.Text.Split('');double tempradius = new double3;for (int i = 0; i < 3; i+)(tempradiusi = Convert.ToDouble(tempi);)pcrdGeo.dLatitude = tempradius0 / 180.0 * pai + tempradius1 / 180.0/ 60.0 * pai + tempradius2 / 180 / 60.0 / 60 * pai;/ 求经度temp = textBox2.Text.Split('');for (int i

7、 = 0; i < 3; i+)(tempradiusi = Convert.ToDouble(tempi);)pcrdGeo.dLongitude = tempradius0 / 180.0 * pai + tempradius1 / 180.0/ 60.0 * pai + tempradius2 / 180 / 60.0 / 60 * pai;int deglon = Convert.ToInt32(pcrdGeo.dLongitude * 180 / pai);/ 求中央经度int num;/带号midlong= 0;/默认值，需要制定分带try(tt= boBox3.Itemsb

8、oBox3.SelectedIndex.ToString();)catch(MessageBox.Show("Choose 3/6 error!");return;)if (tt.CompareTo("3度带")=0)num = Convert.ToInt32(deglon / 6 + 1);midlong = (6 * num - 3) / 180.0 * pai;)if (tt.CompareTo("6度带")=0)(num = Convert.ToInt32(deglon + 1.5) / 3);midlong = num *

9、3 * pai / 180;)double lp=pcrdGeo.dLongitude - midlong;double N = c / Math.Sqrt(1 + epp * epp * Math.Cos(pcrdGeo.dLatitude) * Math.Cos(pcrdGeo.dLatitude);double M = c / Math.Pow(1 + epp * epp * Math.Cos(pcrdGeo.dLatitude) * Math.Cos(pcrdGeo.dLatitude), 1.5);double ita = epp * Math.Cos(pcrdGeo.dLatitu

10、de);double t = Math.Tan(pcrdGeo.dLatitude);double Nscnb = N * Math.Sin(pcrdGeo.dLatitude) * Math.Cos(pcrdGeo.dLatitude);double Ncosb = N * Math.Cos(pcrdGeo.dLatitude);double cosb = Math.Cos(pcrdGeo.dLatitude);double X;double m0, m2, m4, m6, m8;double a0, a2, a4, a6, a8;m0 = a * (1 - ee * ee);m2 = 3.

11、0 / 2.0 * m0 * ee * ee;m4 = 5.0 / 4.0 * ee * ee * m2;m6 = 7.0 / 6.0 * ee * ee * m4;m8 = 9.0 / 8.0 * ee * ee * m6;a0 = m0 + m2 / 2.0 + 3.0 / 8.0 * m4 + 5.0 / 16.0 * m6 + 35.0 / 128.0 * m8;a2 = m2 / 2 + m4 / 2 + 15.0 / 32.0 * m6 +7.0/16.0 * m8;a4 = m4 / 8.0 + 3.0 / 16.0 * m6 + 7.0 / 32.0*m8;a6 = m6 /

12、32.0 + m8 / 16.0;a8 = m8 / 128.0;double B = pcrdGeo.dLatitude;double sb = Math.Sin(B);double cb = Math.Cos(B);double s2b = sb * cb * 2;double s4b = s2b * (1 - 2 * sb * sb) * 2;double s6b = s2b * Math.Sqrt(1 - s4b * s4b) + s4b * Math.Sqrt(1 - s2b * s2b);X = a0 * B - a2 / 2.0 * s2b + a4 * s4b / 4.0 -

13、a6 / 6.0 * s6b;pcrdCar.x = Nscnb * lp * lp / 2.0 + Nscnb * cosb * cosb * Math.Pow(lp, 4)* (5 - t * t + 9 * ita * ita + 4 * Math.Pow(ita, 4) / 24.0+ Nscnb * Math.Pow(cosb, 4) * Math.Pow(lp, 6) * (61 - 58 * t * t + Math.Pow(t, 4) / 720.0 + X;pcrdCar.y = Ncosb * Math.Pow(lp, 1) + Ncosb * cosb * cosb *

14、(1 - t * t + ita * ita) / 6.0 * Math.Pow(lp, 3) + Ncosb * Math.Pow(lp, 5)* Math.Pow(cosb, 4) * (5 - 18* t * t+ Math.Pow(t, 4) + 14 * ita * ita - 58 * ita * ita * t * t) / 120.0 ;if (pcrdCar.y < 0)pcrdCar.y += 500000;richTextBox1.Text = "Results:nX:t" + Convert.ToString(pcrdCar.x) +"

15、;nY:t"+ Convert.ToString(pcrdCar.y);)private void button2_Click(object sender, EventArgs e)(double ee = 0;double a = 0;string tt;int num;/带号string ytext;/利用y值求带号和中央经线)try(tt= boBox2.ItemsboBox2.SelectedIndex.ToString();)catch(MessageBox.Show("Gauss Inverse: Choose datum error!");retur

16、n;)if (tt.CompareTo("克氏椭球")=0)(a = 6378245.00;ee = Math.Sqrt(0.006693421622);)if (tt.CompareTo("WGS-84") = 0)(a = 6378137.00;ee = Math.Sqrt(0.00669437999013);)if (tt.CompareTo("1975国际椭球")=0)(a = 6378140.00;ee = Math.Sqrt(0.006694384999588);)if (tt.CompareTo("2000国家

17、大地坐标系 ")=0)(a = 6378137.0;ee =Math.Sqrt(0.0066943802290);const double pai = 3.1415926535898;double b = Math.Sqrt(a * a * (1 - ee * ee);double c = a * a / b;double epp = Math.Sqrt(a * a - b * b) / b / b);CRDGEODETIC pcrdGeo;CRDCARTESIAN pcrdCar;double midlong = 0;/ 求X, YW带号pcrdCar.x = Convert.To

18、Double(textBox4.Text);ytext = textBox5.Text;string temp = ytext.Substring(0, 2);num = Convert.ToInt32(temp);ytext = ytext.Remove(0, 2);pcrdCar.y = Convert.ToDouble(ytext) - 500000;try(tt= boBox4.ItemsboBox4.SelectedIndex.ToString();)catch(MessageBox.Show("Choose 3/6 error!");return;)if (tt

19、.CompareTo("3度带")=0)(midlong = num * 3 * pai / 180;)if (tt.CompareTo("6度带")=0)(midlong = (6 * num - 3) * pai / 180;)b = Math.Sqrt(a * a * (1 - ee * ee);c = a * a / b;epp = Math.Sqrt(a * a - b * b) / b;double m0, m2, m4, m6, m8;double a0, a2, a4, a6, a8;m0 = a * (1 - ee * ee);m2 =

20、 3.0 / 2.0 * m0 * ee * ee;m4 = 5.0 / 4.0 * ee * ee * m2;m6 = 7.0 / 6.0 * ee * ee * m4;m8 = 9.0 / 8.0 * ee * ee * m6;a0 = m0 + m2 / 2.0 + 3.0 / 8.0 * m4 + 5.0 / 16.0 * m6 + 35.0 / 128.0 * m8;a2 = m2 / 2 + m4 / 2 + 15.0 / 32.0 * m6 +7.0 /16.0 * m8;a4 = m4 / 8.0 + 3.0 / 16.0 * m6 + 7.0/ 32.0* m8;a6 = m

21、6 / 32.0 + m8 / 16.0;a8 = m8 / 128.0;double Bf, B;Bf = pcrdCar.x / a0;B = 0.0;while (Math.Abs(Bf - B) > 1E-10)(B = Bf;double sb = Math.Sin(B);double cb = Math.Cos(B);double s2b = sb * cb * 2;double s4b = s2b * (1 - 2 * sb * sb) * 2;double s6b = s2b * Math.Sqrt(1 - s4b * s4b) + s4b * Math.Sqrt(1 -

22、 s2b* s2b);Bf = (pcrdCar.x - (-a2 / 2.0 * s2b + a4 / 4.0 * s4b - a6 / 6.0 * s6b)/ a0;)double itaf, tf, Vf, Nf;itaf = epp * Math.Cos(Bf);tf = Math.Tan(Bf);Vf = Math.Sqrt(1 + epp * epp * Math.Cos(Bf) * Math.Cos(Bf);Nf = c / Vf;double ynf = pcrdCar.y / Nf;pcrdGeo.dLatitude = Bf - 1.0 / 2.0 * Vf * Vf *tf * (ynf * ynf - 1.0 / 12.0* Math.Pow(ynf, 4) * (5 + 3 * tf * tf + itaf * it