// 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;//结果应该大致是:3736714.783, 627497.303
double OutMyY;
OutMyX =
MyPrj.x;//正算结果: 北坐标x
OutMyY =
MyPrj.y;//结果：东坐标y
//反算
double InputMyX = 3736714.783;//如果是独立计算，应该给出中央子午线L0
double InputMyY = 627497.303;
MyPrj.Setxy(InputMyX, InputMyY);
MyPrj.GetBL(&MyPrj.B, &MyPrj.L);
//把计算出的BL的弧度值换算为dms形式
double OutputMyB;
double OutputMyL;
OutputMyB = MyPrj.B;//反算结果：B
OutputMyL = MyPrj.L;//反算结果：L
//分析表明，此程序的结果和Coord4.2的转换结果是一样的，只差到毫米级
//原程序有几个问题，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的计算公式见孔祥元等主编武汉大学出版社2002年出版的《控制测量学》的第72页
//书的的括号有问题， ( 和 [ 应该交换
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()
//在类外定义构造成员函数，要加上类名和域限定符 ::
{
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()//析构函数，释放构造函数占用的内存
{
}
///
// 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; // 中央子午线经度
double B, L; // 大地坐标
double x, y; // 高斯投影平面坐标
public:
bool BL2xy();
bool xy2BL();
protected:
double a, f, e2, e12; // 基本椭球参数
double A1, A2, A3, A4; //
用于计算X的椭球参数
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();
//定义构造函数，用来初始化。(函数名和类名相同)
~PrjPoint_Krasovsky();
};
class PrjPoint_IUGG1975 : virtual public
PrjPoint
{
public:
PrjPoint_IUGG1975();
~PrjPoint_IUGG1975();
};
#endif /* ndef _COORTRANS_H_INCLUDED
*/
高斯投影正反算C语言程序代码 高斯投影正反算-对网络上面流行的C代码的修改(已正确运行)...
