// GaussBL2xy.cpp : Defines the entry point

for the console application.

//

#include "stdafx.h"

#include "math.h"

#include "CoorTrans.h"

#include

using namespace std;

void main(int argc, char* argv[])

{

double MyL0 = 108; //中央子午線

double MyB = 33.44556666; //33 du 44 fen 55.6666

miao

double MyL = 109.22334444; //3度帶，109 d 22 m

33.4444 s

PrjPoint_Krasovsky MyPrj;

MyPrj.SetL0(MyL0);

MyPrj.SetBL(MyB, MyL);

double

OutMyX;?//結(jié)果應(yīng)該大致是:3736714.783, 627497.303

double OutMyY;

OutMyX =

MyPrj.x;?//正算結(jié)果: 北坐標(biāo)x

OutMyY =

MyPrj.y;?//結(jié)果：東坐標(biāo)y

//反算

double InputMyX = 3736714.783;?//如果是獨立計算，應(yīng)該給出中央子午線L0

double InputMyY = 627497.303;

MyPrj.Setxy(InputMyX, InputMyY);

MyPrj.GetBL(&MyPrj.B, &MyPrj.L);

//把計算出的BL的弧度值換算為dms形式

double OutputMyB;

double OutputMyL;

OutputMyB = MyPrj.B;?//反算結(jié)果：B

OutputMyL = MyPrj.L;?//反算結(jié)果：L

//分析表明，此程序的結(jié)果和Coord4.2的轉(zhuǎn)換結(jié)果是一樣的，只差到毫米級

//原程序有幾個問題，1.Pi的值不對。2.SetBL中多了兩行錯誤代碼

}

double Dms2Rad(double Dms)

{

double Degree, Miniute;

double Second;

int Sign;

double Rad;

if(Dms >= 0)

Sign = 1;

else

Sign = -1;

Dms = fabs(Dms);

Degree = floor(Dms);

Miniute = floor(fmod(Dms * 100.0, 100.0));

Second = fmod(Dms * 10000.0, 100.0);

Rad = Sign * (Degree + Miniute / 60.0 + Second /

3600.0) * PI / 180.0;

return Rad;

}

double Rad2Dms(double Rad)

{

double Degree, Miniute;

double Second;

int Sign;

double Dms;

if(Rad >= 0)

Sign = 1;

else

Sign = -1;

Rad = fabs(Rad * 180.0 / PI);

Degree = floor(Rad);

Miniute = floor(fmod(Rad * 60.0, 60.0));

Second = fmod(Rad * 3600.0, 60.0);

Dms = Sign * (Degree + Miniute / 100.0 + Second /

10000.0);

return Dms;

}

///

// Definition of PrjPoint

///

bool PrjPoint::BL2xy()

{

double X, N, t, t2, m, m2, ng2;

double sinB, cosB;

X = A1 * B * 180.0 / PI + A2 * sin(2 * B) + A3 *

sin(4 * B) + A4 * sin(6 * B);

sinB = sin(B);

cosB = cos(B);

t = tan(B);

t2 = t * t;

N = a / sqrt(1 - e2 * sinB * sinB);

m = cosB * (L - L0);

m2 = m * m;

ng2 = cosB * cosB * e2 / (1 - e2);

//x,y的計算公式見孔祥元等主編武漢大學(xué)出版社2002年出版的《控制測量學(xué)》的第72頁

//書的的括號有問題， ( 和 [ 應(yīng)該交換

x = X + N * t * ((0.5 + ((5 - t2 + 9 * ng2 + 4 *

ng2 * ng2) / 24.0 + (61 -?58 * t2 + t2 * t2) * m2 /

720.0) * m2) * m2);

y = N * m * ( 1 + m2 * ( (1 - t2 + ng2) / 6.0 +

m2 * ( 5 - 18 * t2 + t2 * t2

+ 14 * ng2 - 58 * ng2 * t2 ) /

120.0));

y += 500000;

return true;

}

bool PrjPoint::xy2BL()

{

double sinB, cosB, t, t2, N ,ng2, V, yN;

double preB0, B0;

double eta;

y -= 500000;

B0 = x / A1;

do

{

preB0 = B0;

B0 = B0 * PI / 180.0;

B0 = (x - (A2 * sin(2 * B0) +

A3 * sin(4 * B0) + A4 * sin(6 * B0))) / A1;

eta = fabs(B0 - preB0);

}while(eta > 0.000000001);

B0 = B0 * PI / 180.0;

B = Rad2Dms(B0);

sinB = sin(B0);

cosB = cos(B0);

t = tan(B0);

t2 = t * t;

N = a / sqrt(1 - e2 * sinB * sinB);

ng2 = cosB * cosB * e2 / (1 - e2);

V = sqrt(1 + ng2);

yN = y / N;

B = B0 - (yN * yN - (5 + 3 * t2 + ng2 - 9 * ng2 *

t2) * yN * yN * yN * yN /

12.0 + (61 + 90 * t2 + 45 * t2

* t2) * yN * yN * yN * yN * yN * yN / 360.0)

* V * V * t / 2;

L = L0 + (yN - (1 + 2 * t2 + ng2) * yN * yN * yN

/ 6.0 + (5 + 28 * t2 + 24

* t2 * t2 + 6 * ng2 + 8 * ng2 *

t2) * yN * yN * yN * yN * yN / 120.0) / cosB;

return true;

}

bool PrjPoint::SetL0(double dL0)

{

L0 = Dms2Rad(dL0);

return true;

}

bool PrjPoint::SetBL(double dB, double dL)

{

B = Dms2Rad(dB);

L = Dms2Rad(dL);

//B =

dB;?//我靠，I

wana say fuck

//L = dL; //del it !

BL2xy();

return true;

}

bool PrjPoint::GetBL(double *dB, double *dL)

{

*dB = Rad2Dms(B);

*dL = Rad2Dms(L);

return true;

}

bool PrjPoint::Setxy(double dx, double dy)

{

x = dx;

y = dy;

xy2BL();

return true;

}

bool PrjPoint::Getxy(double *dx, double *dy)

{

*dx = x;

*dy = y;

return true;

}

///

// Definition of PrjPoint_IUGG1975

///

PrjPoint_IUGG1975::PrjPoint_IUGG1975()

//在類外定義構(gòu)造成員函數(shù)，要加上類名和域限定符 ::

{

a = 6378140;

f = 298.257;

e2 = 1 - ((f - 1) / f) * ((f - 1) / f);

e12 = (f / (f - 1)) * (f / (f - 1)) - 1;

A1 = 111133.0047;?//這幾個A是什么意思？

A2 = -16038.5282;

A3 = 16.8326;

A4 = -0.0220;

}

PrjPoint_IUGG1975::~PrjPoint_IUGG1975()?//析構(gòu)函數(shù)，釋放構(gòu)造函數(shù)占用的內(nèi)存

{

}

///

// Definition of PrjPoint_Krasovsky

///

PrjPoint_Krasovsky::PrjPoint_Krasovsky()

{

a = 6378245;

f = 298.3;

e2 = 1 - ((f - 1) / f) * ((f - 1) / f);

e12 = (f / (f - 1)) * (f / (f - 1)) - 1;

A1 = 111134.8611;

A2 = -16036.4803;

A3 = 16.8281;

A4 = -0.0220;

}

PrjPoint_Krasovsky::~PrjPoint_Krasovsky()

{

}

//CoorTrans.h 文件

#ifndef _COORTRANS_H_INCLUDED

#define _COORTRANS_H_INCLUDED

#include "stdafx.h"

#include "math.h"

const double PI =

3.141592653589793;

double Dms2Rad(double Dms);

double Rad2Dms(double Rad);

class PrjPoint

{

public:

double L0; // 中央子午線經(jīng)度

double B, L; // 大地坐標(biāo)

double x, y; // 高斯投影平面坐標(biāo)

public:

bool BL2xy();

bool xy2BL();

protected:

double a, f, e2, e12; // 基本橢球參數(shù)

double A1, A2, A3, A4; //

用于計算X的橢球參數(shù)

public:

bool SetL0(double dL0);

bool SetBL(double dB, double dL);

bool GetBL(double *dB, double *dL);

bool Setxy(double dx, double dy);

bool Getxy(double *dx, double *dy);

};

class PrjPoint_Krasovsky : virtual public

PrjPoint //類的繼承，聲明基類是 PrjPoint，虛基類

{

public:

PrjPoint_Krasovsky();

//定義構(gòu)造函數(shù)，用來初始化。(函數(shù)名和類名相同)

~PrjPoint_Krasovsky();

};

class PrjPoint_IUGG1975 : virtual public

PrjPoint

{

public:

PrjPoint_IUGG1975();

~PrjPoint_IUGG1975();

};

#endif /* ndef _COORTRANS_H_INCLUDED

*/

總結(jié)

以上是生活随笔為你收集整理的高斯投影正反算C语言程序代码,高斯投影正反算-对网络上面流行的C代码的修改(已正确运行)...的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。