gauseprojectview.cpp

大地测量课程实习程序：高斯坐标投影与大地坐标之间的转换

// GauseProjectView.cpp : implementation of the CGauseProjectView class
//

#include "stdafx.h"
#include "GauseProject.h"

#include "GauseProjectDoc.h"
#include "GauseProjectView.h"
/////////////////////////////////
#include "GauseDlg.h"
#include "math.h"
/////////////////////////////////
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

// data
#define  PI       (3.14159265)
#define  rad_dgr  (180.0/PI)
#define  rad_min  (rad_dgr*60)
#define  rad_sec  (rad_min*60)
/*  // kers
#define  _e2    0.006738525
#define  c      6399698.901
*/
//
#define  Maxcnt  5

/////////////////////////////////////////////////////////////////////////////
// CGauseProjectView

IMPLEMENT_DYNCREATE(CGauseProjectView, CView)

BEGIN_MESSAGE_MAP(CGauseProjectView, CView)
	//{{AFX_MSG_MAP(CGauseProjectView)
	ON_COMMAND(ID_GAUSE, OnGause)
	//}}AFX_MSG_MAP
	// Standard printing commands
	ON_COMMAND(ID_FILE_PRINT, CView::OnFilePrint)
	ON_COMMAND(ID_FILE_PRINT_DIRECT, CView::OnFilePrint)
	ON_COMMAND(ID_FILE_PRINT_PREVIEW, CView::OnFilePrintPreview)
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CGauseProjectView construction/destruction

CGauseProjectView::CGauseProjectView()
{
	// TODO: add construction code here

}

CGauseProjectView::~CGauseProjectView()
{
}

BOOL CGauseProjectView::PreCreateWindow(CREATESTRUCT& cs)
{
	// TODO: Modify the Window class or styles here by modifying
	//  the CREATESTRUCT cs

	return CView::PreCreateWindow(cs);
}

/////////////////////////////////////////////////////////////////////////////
// CGauseProjectView drawing

void CGauseProjectView::OnDraw(CDC* pDC)
{
	CGauseProjectDoc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);
	// TODO: add draw code for native data here
}

/////////////////////////////////////////////////////////////////////////////
// CGauseProjectView printing

BOOL CGauseProjectView::OnPreparePrinting(CPrintInfo* pInfo)
{
	// default preparation
	return DoPreparePrinting(pInfo);
}

void CGauseProjectView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO: add extra initialization before printing
}

void CGauseProjectView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO: add cleanup after printing
}

/////////////////////////////////////////////////////////////////////////////
// CGauseProjectView diagnostics

#ifdef _DEBUG
void CGauseProjectView::AssertValid() const
{
	CView::AssertValid();
}

void CGauseProjectView::Dump(CDumpContext& dc) const
{
	CView::Dump(dc);
}

CGauseProjectDoc* CGauseProjectView::GetDocument() // non-debug version is inline
{
	ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CGauseProjectDoc)));
	return (CGauseProjectDoc*)m_pDocument;
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CGauseProjectView message handlers

void CGauseProjectView::OnGause() 
{
	// TODO: Add your command handler code here

	CGauseProjectDoc*  pDoc=GetDocument();
	pDoc->UpdateAllViews(NULL);

	CDC * pDC = (CDC *) GetDC();

	CGauseDlg  dlg;
	if(dlg.DoModal()==IDOK)
	{
		double X,Y;
		float B[4]={0,0,0,0},L[4]={0,0,0,0};  /////////////////////////
		int Daihao,L0;
		double XX;
		int elps=dlg.m_nEllipsoid;
		
		if(dlg.m_nMode==0)	//正算
		{
			B[1] = dlg.m_B1;
			B[2] = dlg.m_B2;
			B[3] = dlg.m_B3;
			L[1] = dlg.m_L1;
			L[2] = dlg.m_L2;
			L[3] = dlg.m_L3;
			Daihao=int( L[1]/6 )+1;
			L0=6*Daihao-3;

			CString strB,strL;
			strB.Format("B = %.0f度%.0f分%.4f秒",B[1],B[2],B[3]);
			strL.Format("L = %.0f度%.0f分%.4f秒",L[1],L[2],L[3]);

			pDC->TextOut(100,100,"正算：");
			pDC->TextOut(100,125,"输入的大地坐标为：");
			pDC->TextOut(100,150,strB);
			pDC->TextOut(100,175,strL);

//			Daihao=int( L[1]/6 )+1;
//			L0=6*Daihao-3;

			Zhengsuan(B,L,L0,X,Y,XX,elps);
//			Y += Daihao*1000000;

			CString strX,strY;
			strX.Format("X = %.3f",X);
			strY.Format("Y = %d %.3f",Daihao,Y);

			pDC->TextOut(100,225,"结果：");
			pDC->TextOut(100,250,strX);
			pDC->TextOut(100,275,strY);
//			pDC->TextOut
		}
		else if(dlg.m_nMode==1)	//反算
		{
			X = dlg.m_X;
			Y = dlg.m_Y;

			CString strX,strY;
			strX.Format("X = %.3f",X);
			strY.Format("Y = %.3f",Y);

			pDC->TextOut(100,100,"反算：");
			pDC->TextOut(100,125,"输入的高斯平面坐标为：");
			pDC->TextOut(100,150,strX);
			pDC->TextOut(100,175,strY);
			
			Fansuan(X,Y,Daihao,B,L,XX,elps);

			CString strB,strL;
			strB.Format("B = %.0f度%.0f分%.4f秒",B[1],B[2],B[3]);
			strL.Format("L = %.0f度%.0f分%.4f秒",L[1],L[2],L[3]);

			pDC->TextOut(100,225,"结果：");
			pDC->TextOut(100,250,strB);
			pDC->TextOut(100,275,strL);
		}
		else if(dlg.m_nMode==2)	//正算后反算
		{
			B[1] = dlg.m_B1;
			B[2] = dlg.m_B2;
			B[3] = dlg.m_B3;
			L[1] = dlg.m_L1;
			L[2] = dlg.m_L2;
			L[3] = dlg.m_L3;
			Daihao=int( L[1]/6 )+1;
			L0=6*Daihao-3;

			CString strB,strL;
			strB.Format("B = %.0f度%.0f分%.4f秒",B[1],B[2],B[3]);
			strL.Format("L = %.0f度%.0f分%.4f秒",L[1],L[2],L[3]);

			pDC->TextOut(100,100,"正算后反算：");
			pDC->TextOut(100,125,"输入的大地坐标为：");
			pDC->TextOut(100,150,strB);
			pDC->TextOut(100,175,strL);

//			Daihao=int( L[1]/6 )+1;
//			L0=6*Daihao-3;

			Zhengsuan(B,L,L0,X,Y,XX,elps);
			Y += Daihao*1000000;
			Fansuan(X,Y,Daihao,B,L,XX,elps);

			strB.Format("B = %.0f度%.0f分%.4f秒",B[1],B[2],B[3]);
			strL.Format("L = %.0f度%.0f分%.4f秒",L[1],L[2],L[3]);
			pDC->TextOut(100,225,"结果：");
			pDC->TextOut(100,250,strB);
			pDC->TextOut(100,275,strL);
		}
	}
}

void CGauseProjectView::Zhengsuan(float *B, float *L, int L0, double &X, double &Y, double &XX,int elps)
{
	double N,V,l2,sBcB,cBcB,tx,t2,etar2,ty,x,y;
	double a,b,_e2,c;
	double B_dgr,B_rad,L_dgr,L_rad,l_rad;
	double bata0, bata2, bata4, bata6, bata8;
	double C0,C2,C4,C6,C8;

	// select ellipsoid
	if(elps==1)			//1975国际椭球体
	{
		a = 6378140.0;
		b = 6356755.2881575287;
	}
	else if(elps==2)	//WGS-84椭球体
	{
		a = 6378137.0;
		b = 6356752.3142;
	}
	else 
	{
		a = 6378245.0;
		b = 6356863.0187330473;
//		_e2 = 0.006738525414683;
//		 c = 6399698.901782711;
	}
	_e2 = (a*a-b*b)/b/b;
	c = a*a/b;
	// degree to rad
	B_dgr= *(B+1) + *(B+2)/60 + *(B+3)/3600 ;
	B_rad= B_dgr / rad_dgr ;
	L_dgr= *(L+1) + *(L+2)/60 + *(L+3)/3600 ;
									//	L_rad= L_dgr / rad_dgr ;	
	l_rad= (L_dgr - L0)/rad_dgr;//rad_sec ;

	// sinB,cosB
	sBcB=sin(B_rad)*cos(B_rad);
	cBcB=cos(B_rad)*cos(B_rad);

	// Calculate XX
	bata0 = 1.0 - 3.0/4.0*_e2 + 45.0/64.0*_e2*_e2 - 175.0/256.0*pow(_e2,3) + 11025.0/16384.0*pow(_e2,4);
	C0 = bata0 * c;
	bata2 = bata0-1.0;
	C2 = bata2 * c;
	bata4 = 15.0/32.0*_e2*_e2 - 175.0/384.0*pow(_e2,3) + 3675.0/8192.0*pow(_e2,4);
	C4 = bata4 * c;
	bata6 = -35./96*pow(_e2,3) + 735./2048*pow(_e2,4);
	C6 = bata6 * c;
	bata8 = 315./1024*pow(_e2,4);
	C8 = bata8 * c;

	XX = C0*B_rad + sBcB*(C2 + C4*cBcB + C6*cBcB*cBcB + C8*pow(cBcB,3) );
//	XX=111134.861*B_dgr-16036.48*sin(2*B_rad)+16.828*sin(4*B_rad)-0.022*sin(6*B_rad);

	V=sqrt(1+_e2*cBcB);
	N=c/V;
	t2=tan(B_rad)*tan(B_rad);
	etar2=_e2*cBcB;

	l2=l_rad*l_rad;

	tx=N*l2*sBcB;
	x=XX+tx/2+tx/24*cBcB*l2*(5-t2+9*etar2+4*etar2*etar2)+tx/720*cBcB*cBcB*l2*l2*(61-58*t2+t2*t2);
	X=x;

	ty=N*l_rad*cos(B_rad);
	y=ty+ty/6*cBcB*l2*(1-t2+etar2)+ty/120*cBcB*cBcB*l2*l2*(5-18*t2+t2*t2+14*etar2-58*etar2*t2);
	Y=500000+y;
}

void CGauseProjectView::Fansuan(double X, double Y, int& Daihao, float *B, float *L, double XX,int elps)
{
	double x,y,cnt,Bf1,Bf,FBf1,Vf,Nf,Mf,tf2,etarf2,tb,tl,Nf2,y2;
	double a,b,e2,_e2,c;
	double bata0,C0;
	double m0,m2,m4,m6,m8,a2,a4,a6,a8;

	// select ellipsoid
	if(elps==1)			//1975国际椭球体
	{
		a = 6378140.0;
		b = 6356755.2881575287;
	}
	else if(elps==2)	//WGS-84椭球体
	{
		a = 6378137.0;
		b = 6356752.3142;
	}
	else 
	{
		a = 6378245.0;
		b = 6356863.0187330473;
//		_e2 = 0.006738525;
//		 c = 6399698.901;
	}
	e2 = (a*a-b*b)/a/a;
	_e2 = (a*a-b*b)/b/b;
	c = a*a/b;

	// Calculate XX
	bata0 = 1 - 3./4*_e2 + 45./64*_e2*_e2 - 175./256*pow(_e2,3) + 11025./16384*pow(_e2,4);
	C0 = bata0 * c;
	//
	m0 = a*(1-e2);
	m2 = 3.0/2*e2*m0;
	m4 = 5.0/4*e2*m2;
	m6 = 7.0/6*e2*m4;
	m8 = 9.0/8*e2*m6;

	a2 = m2/2 + m4/2 + 15.0/32*m6 + 7.0/16*m8;
	a4 = m4/8 + 3.0/16*m6 + 7.0/32*m8;
	a6 = m6/32 + m8/16;
	a8 = m8/128;

	x = X;
	Daihao = int(Y/1000000);
	y = Y - Daihao*1000000 - 500000;
	Bf=X/C0;//111134.8611;//C0;//
	cnt=0;
	do
	{
		Bf1=Bf;
		FBf1=-a2/2*sin(2*Bf1)+a4/4*sin(4*Bf1)-a6/6*sin(6*Bf1);//+a8/8*sin(8*Bf1);
//		FBf1=-16036.4803*sin(2*Bf1)+16.8281*sin(4*Bf1)-0.022*sin(6*Bf1);
		Bf=(X-FBf1)/C0;//111134.8611;
		cnt++;
	}while( fabs(Bf-Bf1)>=0.0000001 && cnt<Maxcnt );
	if( cnt==Maxcnt && fabs(Bf-Bf1)>=0.0000001 )
	{
		MessageBox("Error in counting Bf!");
		system("pause");
		exit(0);
	}else{
		//
		Vf=sqrt(1+_e2*cos(Bf)*cos(Bf));
		Nf=c/Vf;
		Nf2=Nf*Nf;
		Mf=c/(Vf*Vf*Vf);
		tf2=tan(Bf)*tan(Bf);
		etarf2=_e2*cos(Bf)*cos(Bf);

		y2=y*y;
		tb=tan(Bf)/(Mf*Nf)*y2;
		*B=Bf-tb/2+tb/(24*Nf2)*y2*(5+3*tf2+etarf2-9*etarf2*tf2)-tb/(720*Nf2*Nf2)*y2*y2*(61+90*tf2+45*tf2*tf2);
		*(B+1)=int( *B * rad_dgr );
		*(B+2)=int( *B * rad_min - *(B+1)*60 );
		*(B+3)= *B * rad_sec - *(B+1)*3600 - *(B+2)*60 ;

		tl=y/Nf/cos(Bf);
		*L=tl-tl/(6*Nf2)*y2*(1+2*tf2+etarf2)+tl/(120*Nf2*Nf2)*y2*y2*(5+28*tf2+24*tf2*tf2+6*etarf2+8*etarf2*tf2);
//		*L=20*6-3+l;
		double  t;
		t = ( (6*Daihao-3) / rad_dgr + *L ) * rad_dgr ;
		*(L+1)= int( t );
		t = (t- *(L+1)) * 60;//rad_min ;
		*(L+2)= int( t );
		t = (t- *(L+2)) * 60;//rad_sec ;
		*(L+3)= t ;
	}
}