高斯投影正反算c代码

高斯投影正反算程序设计
一．程序设计流程
本程序(de)设计思路如下：
（1）,程序采用VS08版本作为开发平台,并采用C语言作为开发语言,设计为WindowsForm窗体程序形式.
（2）,程序主要(de)算法来自于教材.但是本程序为了更加实用,添加了更多(de)解算基准,包括：WGS-84,国际椭球1975,克氏椭球,和2000国家大地坐标系.
（3）,程序为了更方便(de)读取数据和输出数据,故需要自己定义了固定(de)数据输入格式和数据输出格式或形式,请老师注意查看.
二．代码
using System;
using Syst
using SystemponentModel;
using System.Data;
using System.Drawing;
using System.Text;
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
{
}
catch
{
MessageBox.Show("Gauss Inverse: Choose datum error");
return;
}
if (ttpareTo("克氏椭球")==0)
{
a = 6378245.00;
}
if (ttpareTo("WGS-84") == 0)
{
a = 6378.00;
}
if (ttpareTo("1975国际椭球") == 0)
{
a = 6378140.00;
ee =
}
if (ttpareTo("2000国家大地坐标系") == 0)
{
a = 6378.0;
}
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;
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;
//求经度
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
{
}
catch
{
MessageBox.Show("Choose 3/6 error");
return;
}
if (ttpareTo("3度带") == 0)
{
num = Convert.ToInt32(deglon / 6 + 1);
midlong = (6 num - 3) / 180.0 pai;
}
if (ttpareTo("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
{
}
catch
{
MessageBox.Show("Gauss Inverse: Choose datum
error");
return;
}
if (ttpareTo("克氏椭球") == 0)
{
a = 6378245.00;
}
if (ttpareTo("WGS-84") == 0)
{
a = 6378.00;
}
if (ttpareTo("1975国际椭球") == 0)
{
a = 6378140.00;
}
if (ttpareTo("2000国家大地坐标系") == 0)
{
a = 6378.0;
}
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
{
}
catch
{
MessageBox.Show("Choose 3/6 error");
return;
}
if (ttpareTo("3度带") == 0)
{
midlong = num 3 pai / 180;
}
if (ttpareTo("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) {
}
}
}
三．程序运行结果分析
通过选取书上(de)具体实例进行测试,本程序(de)精度大体满足要求,一般正算(de)精度在0.01米和0.001米之间,反算(de)精度在0.0001秒左右.以下是程序运行(de)截图.。