高斯投影正反算-c#代码.doc

精品文档高斯投影正反算程序设计一程序设计流程本程序的设计思路如下：（1），程序采用VS08版本作为开发平台，并采用C#语言作为开发语言，设计为WindowsForm窗体程序形式。（2），程序主要的算法来自于教材。但是本程序为了更加实用，添加了更多的解算基准，包括：WGS-84，国际椭球1975，克氏椭球，和2000国家大地坐标系。（3），程序为了更方便的读取数据和输出数据，故需要自己定义了固定的数据输入格式和数据输出格式或形式，请老师注意查看。二代码using System;using System.Collections.Generic;using System.ComponentModel;using System.Data;using System.Drawing;using System.Text;using System.Windows.Forms;namespace Gauss public partial class Form1 : Form /大地坐标 /Geodetic Coordinate public struct CRDGEODETIC public double dLongitude; public double dLatitude; public double dHeight; /笛卡尔坐标 /Cartesian Coordinate public struct CRDCARTESIAN 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.IboBox1.SelectedIndex.ToString(); catch MessageBox.Show(Gauss Inverse: Choose datum 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(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 = textBox1.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 = 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.IboBox3.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 * 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.dLatitude); 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.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 / 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 - 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 * (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 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 - 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 * itaf - 9 * Math.Pow(itaf * tf, 2) + 1.0 / 360.0 * (61 + 90 * tf * tf + 45 * Math.Pow(tf, 4) * Math.Pow(ynf, 6); pcrdGeo.dLongitude = (ynf / Math.Cos(Bf) - (1 + 2 * tf * tf + itaf * itaf) * Math.Pow(ynf, 3) / 6.0 / Math.Cos(Bf) + (5 + 28 * tf * tf + 24 * Math.Pow(tf, 4) + 6 * itaf * itaf + 8 * Math.Pow(itaf * tf, 2) * Math.Pow(ynf, 5) / 120.0 / Math.Cos(Bf); pcrdGeo.dLongitude = pcrdGeo.dLongitude + midlong; /pcrdGeo.dLatitude = pcrdGeo.dLatitude; richTextBox2.Text = Results:nLatitude: + Convert.ToString(pc