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

高斯投影正反算程序设计
一．程序设计流程
本程序的设计思路如下：
（1），程序采用VS08版本作为开发平台，并采用C#语言作为开发语言，设计为WindowsForm 窗体程序形式。

（2），程序主要的算法来自于教材。

但是本程序为了更加实用，添加了更多的解算基准，包括：WGS-84，国际椭球1975，克氏椭球，和2000国家大地坐标系。

（3），程序为了更方便的读取数据和输出数据，故需要自己定义了固定的数据输入格式和数据输出格式或形式，请老师注意查看。

二．代码
using System;
using System.Collections.Generic;
using ponentModel;
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.Items[boBox1.SelectedIndex].ToString(); }
catch
{
MessageBox.Show("Gauss Inverse: Choose datum error!");
return;
}
if (pareTo("克氏椭球")==0)
{
a = 6378245.00;
ee = Math.Sqrt(0.006693421622);
}
if (pareTo("WGS-84") == 0)
{
a = 6378137.00;
ee = Math.Sqrt(0.00669437999013);
}
if (pareTo("1975国际椭球") == 0)
{
a = 6378140.00;
ee = Math.Sqrt(0.006694384999588);
}
if (pareTo("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 double[3];
for (int i = 0; i < 3; i++)
{
tempradius[i] = Convert.ToDouble(temp[i]);
}
pcrdGeo.dLatitude = tempradius[0] / 180.0 * pai + tempradius[1] / 180.0 / 60.0 * pai + tempradius[2] / 180 / 60.0 / 60 * pai;
//求经度
temp = textBox2.Text.Split(' ');
for (int i = 0; i < 3; i++)
{
tempradius[i] = Convert.ToDouble(temp[i]);
}
pcrdGeo.dLongitude = tempradius[0] / 180.0 * pai + tempradius[1] / 180.0 / 60.0 * pai + tempradius[2] / 180 / 60.0 / 60 * pai;
int deglon = Convert.ToInt32(pcrdGeo.dLongitude * 180 / pai);
//求中央经度
int num; //带号 midlong = 0; //默认值，需要制定分带
try
{
tt = boBox3.Items[boBox3.SelectedIndex].ToString();
}
catch
{
MessageBox.Show("Choose 3/6 error!");
return;
}
if (pareTo("3度带") == 0)
{
num = Convert.ToInt32(deglon / 6 + 1);
midlong = (6 * num - 3) / 180.0 * pai;
}
if (pareTo("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 < 0)
pcrdCar.y += 500000;
richTextBox1.Text = "Results:\nX:\t" + Convert.ToString(pcrdCar.x) +"\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.Items[boBox2.SelectedIndex].ToString(); }
catch
{
MessageBox.Show("Gauss Inverse: Choose datum error!");
return;
}
if (pareTo("克氏椭球") == 0)
{
a = 6378245.00;
ee = Math.Sqrt(0.006693421622);
}
if (pareTo("WGS-84") == 0)
{
a = 6378137.00;
ee = Math.Sqrt(0.00669437999013);
}
if (pareTo("1975国际椭球") == 0)
{
a = 6378140.00;
ee = Math.Sqrt(0.006694384999588);
}
if (pareTo("2000国家大地坐标系") == 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，Y和带号
pcrdCar.x = Convert.ToDouble(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.Items[boBox4.SelectedIndex].ToString(); }
catch
{
MessageBox.Show("Choose 3/6 error!");
return;
}
if (pareTo("3度带") == 0)
{
midlong = num * 3 * pai / 180;
}
if (pareTo("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 = 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 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 - 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(pcrdGeo.dLatitude) + "\nLongtitude: " +
Convert.ToString(pcrdGeo.dLongitude);
}
private void label13_Click(object sender, EventArgs e)
{
}
}
}
三．程序运行结果分析
通过选取书上的具体实例进行测试，本程序的精度大体满足要求，一般正算的精度在0.01米和0.001米之间，反算的精度在0.0001秒左右。

以下是程序运行的截图。