coordtrans.cpp

坐标转换程序

// CoordTrans.cpp: implementation of the CCoordTrans class.
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "CASMCoordTrans.h"
#include "CoordTrans.h"
#include <math.h>

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

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CCoordTrans::CCoordTrans(double BX,double LY,double HZ)
{
	XB=BX;
	YL=LY;
	ZH=HZ;
	PI =3.141592653589793;

	t1=97.5;                       // 三参数法转换参数
	t2=58.8;
	t3=6.3;

	arf84 =1.0/298.257223563;     // WGS84椭球参数
	a84 =6378137.0;
	arf80 =1.0/298.257;           // Xi'an80椭球参数
	a80 =6378140.0;
}

CCoordTrans::~CCoordTrans()
{
	
}
//////////////////////////////////////////////////////////////////////
// Function 1-4
//////////////////////////////////////////////////////////////////////

double CCoordTrans::jdhhd(double b)   //degree to radian
{
	double z,b1,b2,b3;
	z = 3.1415926535897932;
	if (b>1.0e-12)
	{
		b1 = b-int(b) ; 
		b2 = int(b1*100)*1.0/100;
		b3 = b1-b2;
		b = int(b)+b2*100.0/60+b3*10000.0/3600;
		b= b*z/180;
	} 
	else
	{
		b = fabs(b);
		b1 = b- int(b) ; 
		b2 = int(b1*100)*1.0/100;
		b3 = b1-b2;
		b = int(b)+b2*100.0/60+b3*10000.0/3600;
		b = b*z/180;
		b = -b;
	}
	return b;
}

double CCoordTrans::hdhjd(double b)   //radian to degree
{
	double z,b1,b2,b3,b4;
	z = 3.1415926535897932;
	if (b>1.0e-12)
	{
		b =b*180.0/z;
		b1 =b-int(b) ; 
		b2 =b1*60.0/100;
		b3 =int(b2*100)*1.0/100;
		b4 =(b2-b3)*60.0/100;
		b =int(b)+b3+b4;
	} 
	else
	{
		b=fabs(b);
		b=b*180.0/z;
		b1=b- int(b) ; 
		b2=b1*60.0/100;
		b3=int(b2*100)*1.0/100;
		b4=(b2-b3)*60.0/100;
		b=int(b)+b3+b4;
		b=-b;
	}
	return b;
}

void CCoordTrans::blhxyz(double b84, double l84, double h84)  //(B84,L84,H84) to (X84,Y84,Z84)
{
	double e2,a2,b2,dn;
	e2 =arf84*(2-arf84);
	a2 =a84*a84;
	b2 =a2*(1-e2);
	dn =a2/sqrt(a2*cos(b84)*cos(b84)+b2*sin(b84)*sin(b84));
	x84 =(dn+h84)*cos(b84)*cos(l84);
	y84 =(dn+h84)*cos(b84)*sin(l84);
	z84 =(dn*(1-e2)+h84)*sin(b84);
} 

void CCoordTrans::trans(double x84, double y84,double z84)   // (X84,Y84,Z84) to (x80,y80,z80)
{
	x80 =x84+t1;
	y80 =y84+t2;
	z80 =z84+t3;
}

void CCoordTrans::xyzblh(double x80, double y80, double z80)  // (x80,y80,z80) to (B80,L80,H80)
{
	double e2,ep2,r,sinfi2,sin2fi2;
	double p,q,pn,tgb,atanb;
	double xy;            // y/x
	double atanxy;
	e2 = arf80*(2-arf80);
	ep2 = e2/(1-e2);
	r = sqrt(x80*x80+y80*y80+z80*z80);
	sinfi2 =(z80*z80)/(r*r);
    sin2fi2 =4*(sinfi2-sinfi2*sinfi2);
	q =e2/(2-e2);
	p =a80/sqrt(1+ep2*sinfi2);
	h80 =(r-p)*(1-(1.0/2)*q*q*sin2fi2);
	pn =a80/sqrt((1-e2*(sinfi2+q*sin2fi2)));
	tgb =z80/sqrt(x80*x80+y80*y80)*(1+e2/(1-e2+h80/pn));
	atanb =atan(fabs(tgb));
	double pi=3.1415926535897932;
	if(tgb<0)
	{ 
		b80=pi-atanb;
	}
	else
	{
		b80=atanb;
	}
	xy=y80/x80;
	atanxy=atan(fabs(xy));
	if (xy<0)
	{
		l80=pi-atanxy;
	} 
	else
	{
		l80=atanxy;
	}
}