projecttran.cpp

坐标转换工具 WGS84坐标转换为地方坐标系的程序 以及各种三维坐标系之间的转换

// ProjectTran.cpp: implementation of the CProjectTran class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "ProjectTran.h"
#include "math.h"
#include "common.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CProjectTran::CProjectTran()
{

	L0=0;
	m_ntype=CPROJ_TRAN;
}

CProjectTran::~CProjectTran()
{

}

//高斯坐标正算;
void CProjectTran::GussCoordTranNormal()
{
//	double C0,C1,C2,C3;

	//临时变量
	double A,B,C,D,E,F;		
	double t1,t2,t3,t4,t5,t6;
	double e4,e6,e8,e10,e12;
	double e2t;		//1-e2

	sax=lax-af*lax;
	e2=2*af-af*af;



	e4=e2*e2; e6=e4*e2; e8=e6*e2; e10=e8*e2; e12=e10*e2;
		
	A=1+3*e2/4+45*e4/64+175*e6/256+11025*e8/16384+43659*e10/65536;
	B=3*e2/4+15*e4/16+525*e6/512+2205*e8/2048+72765*e10/65536;
	C=15*e4/64+105*e6/256+2205*e8/4096+10395*e10/16384;
	D=35*e6/512+315*e8/2048+31185*e10/131072;
	E=315*e8/16384+3465*e10/65536;
	F=693*e10/131072;

	double A1,B1,C1,D1,E1,F1,G1;
	A1=A+693693*e12/1048576;
	B1=3*e2/8+15*e4/32+525*e6/1024+2205*e8/4096+72765*e10/131072+297297*e12/524288;
	C1=15*e4/256+105*e6/1024+2205*e8/16384+10395*e10/65536+1486485*e12/8388608;
	D1=35*e6/3072+105*e8/4096+10395*e10/262144+55055*e12/1048576;
	E1=315*e8/131072+3465*e10/524288+99099*e12/8388608;
	F1=693*e10/131072+9009*e12/5242880;
	G1=1001*e12/8388608;

	e2t=1-e2;
/*	t1=A*lax*e2t;	t2=-1*B*lax*e2t/2;
	t3=C*lax*e2t/4;	t4=-1*D*lax*e2t/6;
	t5=E*lax*e2t/8; t6=-1*F*lax*e2t/10;

	C0=t1;	C1=2*t2+4*t3+6*t4;  
	C2=-1*(8*t3+32*t4); 
	C3=32*t4; 
*/
	


	double l,t,m0,N;
	double t7;        //
	double BB,LL;
	xyz temp;

	m_xyz.RemoveAll();

	for(int i=0;i<m_blh.GetSize();i++)
	{
		BB=m_blh.GetAt(i).b;
		LL=m_blh.GetAt(i).l;

		temp.x=lax*e2t*(A1*BB-B1*sin(2*BB)+C1*sin(4*BB)-D1*sin(6*BB)+E1*sin(8*BB)-F1*sin(10*BB)+G1*sin(12*BB));
		
		l=LL-L0;
		t=tan(BB);
		m0=l*cos(BB);
		t7=e2*(1+cos(2*BB))/(2*e2t);
		N=lax/sqrt(1-e2*(1-cos(2*BB))/2);

		//temp.x=C0*BB+cos(BB)*(C1*sin(BB)+C2*(3*sin(BB)-sin(3*BB))/4+C3*(sin(5*BB)-5*sin(3*BB)+10*sin(BB))/16);
		temp.x=temp.x+N*t*m0*m0/2+(5-t*t+9*t7+4*t7*t7)*N*t*m0*m0*m0*m0/24;
		temp.x=temp.x+(61-58*t*t+t*t*t*t)*N*t*m0*m0*m0*m0*m0*m0/720;

		temp.y=N*m0+(1-t*t+t7)*N*m0*m0*m0/6;
		temp.y=temp.y+(5-18*t*t+t*t*t*t+14*t7-58*t7*t*t)*N*m0*m0*m0*m0*m0/120;
		temp.y+=500000;

		temp.z=0.00;


		m_xyz.Add(temp);	


	}

}

void CProjectTran::GussCoordTranRevers()
{
	

	//临时变量
	double A,B,C,D,E,F;		
	double t1,t11,t2,t3,t4;
	double e4,e6,e8,e10,e12;
	double e2t;		//1-e2

	sax=lax-af*lax;
	e2=2*af-af*af;


	e4=e2*e2; e6=e4*e2; e8=e6*e2; e10=e8*e2; e12=e10*e2;
		
	A=1+3*e2/4+45*e4/64+175*e6/256+11025*e8/16384+43659*e10/65536;
	B=3*e2/4+15*e4/16+525*e6/512+2205*e8/2048+72765*e10/65536;
	C=15*e4/64+105*e6/256+2205*e8/4096+10395*e10/16384;
	D=35*e6/512+315*e8/2048+31185*e10/131072;
	E=315*e8/16384+3465*e10/65536;
	F=693*e10/131072;

	double AA,BB,CC,DD,EE,FF,GG;
	AA=A+693693*e12/1048576;
	BB=3*e2/8+15*e4/32+525*e6/1024+2205*e8/4096+72765*e10/131072+297297*e12/524288;
	CC=15*e4/256+105*e6/1024+2205*e8/16384+10395*e10/65536+1486485*e12/8388608;
	DD=35*e6/3072+105*e8/4096+10395*e10/262144+55055*e12/1048576;
	EE=315*e8/131072+3465*e10/524288+99099*e12/8388608;
	FF=693*e10/131072+9009*e12/5242880;
	GG=1001*e12/8388608;

	e2t=1-e2;
	t1=lax*e2t; t11=t1*AA; t1=t1*A;
	t2=-1*B*lax*e2t/2;
	t3=C*lax*e2t/4;	t4=-1*D*lax*e2t/6;
	t2=-t2/t1; t3=-t3/t1;
	t4=-t4/t1; 

	double tm1,tm2,tm3;
	double tm21,tm22,tm23;
	double tm31,tm32,tm33;
	double K1,K2,K3;
	double bf0,bf,vf;

	tm1=t3+t2*t2;
	tm2=t2+t2*t3-3*(t2*t2*t2)/2-2*t2*t3;
	tm3=t4+t2*t3+(t2*t2*t2)/2+2*t3*t2;

	tm21=t3+t2*tm2;
	tm22=t2+t2*tm1-3*(t2*tm2*tm2)/2-2*t3*tm2;
	tm23=t4+t2*tm1+(t2*tm2*tm2)/2+2*t3*tm2;

	tm31=t3+t2*tm22;
	tm32=t2+t2*tm21-3*(t2*tm22*tm22)/2-2*t3*tm22;
	tm33=t4+t2*tm21+(t2*tm22*tm22)/2+2*t3*tm22;

	K1=2*tm32+4*tm31+6*tm33;
	K2=8*tm31+32*tm33;
	K3=32*tm33;

	m_blh.RemoveAll();

	double X,Y;
	double N,t7,t8,tf;
	blh temp;
	double bf01;

	double fb;


	for(int i=0;i<m_xyz.GetSize();i++)
	{
		X=m_xyz.GetAt(i).x;
		Y=m_xyz.GetAt(i).y-500000;

	
		bf0=X/t11;
		fb=lax*e2t*((-1)*BB*sin(2*bf0)+CC*sin(4*bf0)-DD*sin(6*bf0)+EE*sin(8*bf0)-FF*sin(10*bf0)+GG*sin(12*bf0));
		bf01=(X-fb)/t11;
	
		while(fabs(bf0-bf01)>1.0E-16)
		{

			bf0=bf01;	
			fb=lax*e2t*((-1)*BB*sin(2*bf0)+CC*sin(4*bf0)-DD*sin(6*bf0)+EE*sin(8*bf0)-FF*sin(10*bf0)+GG*sin(12*bf0));
			bf01=(X-fb)/t11;

		}
	
		bf=bf0;
	
		N=lax/sqrt(1-e2*(1-cos(2*bf))/2);
		t7=e2*(1+cos(2*bf))/(2*(1-e2));

		vf=sqrt(1+t7);
		tf=tan(bf);

		t8=Y/N;

		temp.b=bf-vf*vf*tf*t8*t8/2;
		temp.b=temp.b+(5+3*tf*tf+t7-9*t7*tf*tf)*t8*t8*t8*t8*vf*vf*tf/24;
		temp.b=temp.b-(61+90*tf*tf+45*tf*tf*tf*tf)*vf*vf*tf*t8*t8*t8*t8*t8*t8/720;

		temp.l=t8/cos(bf)-(1+2*tf*tf+t7)*t8*t8*t8/(6*cos(bf));
		temp.l=temp.l+(5+28*tf*tf+24*tf*tf*tf*tf+6*t7+8*t7*tf*tf)*t8*t8*t8*t8*t8/(120*cos(bf));
		temp.l=temp.l+L0;

		temp.h=0.00;
		
		m_blh.Add(temp);


	}


}


void CProjectTran::SetL0(double m_L0)
{
	L0=m_L0;

}