高斯投影正反算 c#代码

高斯投影正反算程序设计

一．程序设计流程

本程序的设计思路如下：

（1），程序采用VS08版本作为开发平台，并采用C#语言作为开发语言，设计为WindowsForm 窗体程序形式。

（2），程序主要的算法来自于教材。但是本程序为了更加实用，添加了更多的解算基准，包括：WGS-84，国际椭球1975，克氏椭球，和2000国家大地坐标系。

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

二．代码

using System;

using Syst

using ponentModel;

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 (pareTo("克氏椭球")==0)

{

a = 6378245.00;

}

if (pareTo("WGS-84") == 0)

{

a = 6378137.00;

}

if (pareTo("1975国际椭球") == 0)

{

a = 6378140.00;

ee =

}

if (pareTo("2000国家大地坐标系") == 0)

{

a = 6378137.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 (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;