高斯投影C++代码

高斯投影C++代码

CCoordinate CEarthPlaneDialog::EarthToPlane(const CEarth & earth) {

double B = earth.m_dB.m_realRadians; double L = earth.m_dL.m_realRadians; L = SetRad_PIPI(L);

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

double L0 = SetL0(L); 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);

double 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);

double 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; }

//将大地坐标系转换为平面直角坐标系

CEarth CEarthPlaneDialog::PlaneToEarth(const CCoordinate & coor, UINT Sign) {

double x = coor.m_dX; double y = coor.m_dY; CCoordinate coor(x,y); return coor;

double L0 = (6.0 * Sign - 3.0) / 180 * PI;

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.0000000000000001); B0 = B0 * PI / 180.0; 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;

double 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;

double 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;

CEarth earth(B,L,0,'r'); }

return earth;

//将平面直角坐标系转换为大地坐标系

void CEarthPlaneDialog::OnRADIOPlane1954() {

// TODO: Add your control notification handler code here a = 6378245;

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

f = 298.3;

A1 = 111134.8611; A2 = -16036.4803; A3 = 16.8281; A4 = -0.0220; }

void CEarthPlaneDialog::OnRADIOPlane1980() {

// TODO: Add your control notification handler code here a = 6378140; f = 298.257;

e2 = 1 - ((f - 1) / f) * ((f - 1) / f); e12 = (f / (f - 1)) * (f / (f - 1)) - 1; A1 = 111133.0047; A2 = -16038.5282; A3 = 16.8326; A4 = -0.0220; }

void CEarthPlaneDialog::OnRADIOPlaneWGS84() {

// TODO: Add your control notification handler code here a = 6378137;

m_iReferenceFrame = 1; m_iReferenceFrame = 0;

f = 298.257223563;

e2 = 1 - ((f - 1) / f) * ((f - 1) / f); e12 = (f / (f - 1)) * (f / (f - 1)) - 1; A1 = 111133.0047; A2 = -16038.5282; A3 = 16.8326; A4 = -0.0220; }

double CEarthPlaneDialog::SetL0(double L) {

else

else if( L < 0 ) { L = -L;

for(m_iCincture = 1 ; !(( L > (6 * m_iCincture - 6 )) & ( L <= ( 6 * {}

L0 = m_iCincture * 6 - 363; }

if( L > 0 ) {

for(m_iCincture = 1 ; !(( L >= (6 * m_iCincture - 6 )) & ( L < ( 6 * {}

L0 = 6 * m_iCincture - 3; }

L = L / PI * 180;

m_iReferenceFrame = 2;

m_iCincture ))); m_iCincture++)

m_iCincture ))); m_iCincture++) m_iCincture = 61 - m_iCincture;

}

{ L0 = 3; }

double temp = ((double)L0) / 180 * PI; return temp;

m_iCincture = 1;

double CEarthPlaneDialog::SetRad_PIPI(double Rad) {

if(Rad < 0) }

if(Rad > 0) {

double i = floor(Rad / PI / 2); }

Rad -= 2 * PI * i;

{

Rad = -Rad;

double i = floor(Rad / PI / 2); Rad = -Rad + PI * 2; }

if(Rad > PI)

return Rad;

Rad -= 2 * PI;

Rad -= 2 * PI * i;