gpsdlg.cpp

本程序是用于WGS-84向北京1954坐标系转换的。要想进行转换需要至少3个控制点。这些控制点要已知他们在WGS-84和北京1954坐标系空间直角三维坐标。程序通过解算

// GPSDlg.cpp : implementation file
//

#include "stdafx.h"
#include "GPS.h"
#include "GPSDlg.h"
#include "math.h"
#define PI 3.14159265358979323846264338327950
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// CAboutDlg dialog used for App About

class CAboutDlg : public CDialog
{
public:
	CAboutDlg();

// Dialog Data
	//{{AFX_DATA(CAboutDlg)
	enum { IDD = IDD_ABOUTBOX };
	//}}AFX_DATA

	// ClassWizard generated virtual function overrides
	//{{AFX_VIRTUAL(CAboutDlg)
	protected:
	virtual void DoDataExchange(CDataExchange* pDX);    // DDX/DDV support
	//}}AFX_VIRTUAL

// Implementation
protected:
	//{{AFX_MSG(CAboutDlg)
	//}}AFX_MSG
	DECLARE_MESSAGE_MAP()
};

CAboutDlg::CAboutDlg() : CDialog(CAboutDlg::IDD)
{
	//{{AFX_DATA_INIT(CAboutDlg)
	//}}AFX_DATA_INIT
}

void CAboutDlg::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CAboutDlg)
	//}}AFX_DATA_MAP
}

BEGIN_MESSAGE_MAP(CAboutDlg, CDialog)
	//{{AFX_MSG_MAP(CAboutDlg)
		// No message handlers
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CGPSDlg dialog

CGPSDlg::CGPSDlg(CWnd* pParent /*=NULL*/)
	: CDialog(CGPSDlg::IDD, pParent)
{
	//{{AFX_DATA_INIT(CGPSDlg)
	m_B = 30.302622;
	m_BB = 0.0;
	m_H = 47.360;
	m_HH = 0.0;
	m_L = 114.212196;
	m_LL = 0.0;
	m_X = 0.0;
	m_XX = 0.0;
	m_Y = 0.0;
	m_YY = 0.0;
	m_Z = 0.0;
	m_L0 = 0.0;
	//}}AFX_DATA_INIT
	// Note that LoadIcon does not require a subsequent DestroyIcon in Win32
	m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);
}

void CGPSDlg::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CGPSDlg)
	DDX_Text(pDX, IDC_B, m_B);
	DDX_Text(pDX, IDC_BB, m_BB);
	DDX_Text(pDX, IDC_H, m_H);
	DDX_Text(pDX, IDC_HH, m_HH);
	DDX_Text(pDX, IDC_L, m_L);
	DDX_Text(pDX, IDC_LL, m_LL);
	DDX_Text(pDX, IDC_X, m_X);
	DDX_Text(pDX, IDC_xx, m_XX);
	DDX_Text(pDX, IDC_Y, m_Y);
	DDX_Text(pDX, IDC_yy, m_YY);
	DDX_Text(pDX, IDC_Z, m_Z);
	DDX_Text(pDX, IDC_L0, m_L0);
	//}}AFX_DATA_MAP
}

BEGIN_MESSAGE_MAP(CGPSDlg, CDialog)
	//{{AFX_MSG_MAP(CGPSDlg)
	ON_WM_SYSCOMMAND()
	ON_WM_PAINT()
	ON_WM_QUERYDRAGICON()
	ON_BN_CLICKED(IDOK, OnCalculate)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CGPSDlg message handlers

BOOL CGPSDlg::OnInitDialog()
{
	CDialog::OnInitDialog();

	// Add "About..." menu item to system menu.

	// IDM_ABOUTBOX must be in the system command range.
	ASSERT((IDM_ABOUTBOX & 0xFFF0) == IDM_ABOUTBOX);
	ASSERT(IDM_ABOUTBOX < 0xF000);

	CMenu* pSysMenu = GetSystemMenu(FALSE);
	if (pSysMenu != NULL)
	{
		CString strAboutMenu;
		strAboutMenu.LoadString(IDS_ABOUTBOX);
		if (!strAboutMenu.IsEmpty())
		{
			pSysMenu->AppendMenu(MF_SEPARATOR);
			pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);
		}
	}

	// Set the icon for this dialog.  The framework does this automatically
	//  when the application's main window is not a dialog
	SetIcon(m_hIcon, TRUE);			// Set big icon
	SetIcon(m_hIcon, FALSE);		// Set small icon
	
	// TODO: Add extra initialization here
	
	return TRUE;  // return TRUE  unless you set the focus to a control
}

void CGPSDlg::OnSysCommand(UINT nID, LPARAM lParam)
{
	if ((nID & 0xFFF0) == IDM_ABOUTBOX)
	{
		CAboutDlg dlgAbout;
		dlgAbout.DoModal();
	}
	else
	{
		CDialog::OnSysCommand(nID, lParam);
	}
}

// If you add a minimize button to your dialog, you will need the code below
//  to draw the icon.  For MFC applications using the document/view model,
//  this is automatically done for you by the framework.

void CGPSDlg::OnPaint() 
{
	if (IsIconic())
	{
		CPaintDC dc(this); // device context for painting

		SendMessage(WM_ICONERASEBKGND, (WPARAM) dc.GetSafeHdc(), 0);

		// Center icon in client rectangle
		int cxIcon = GetSystemMetrics(SM_CXICON);
		int cyIcon = GetSystemMetrics(SM_CYICON);
		CRect rect;
		GetClientRect(&rect);
		int x = (rect.Width() - cxIcon + 1) / 2;
		int y = (rect.Height() - cyIcon + 1) / 2;

		// Draw the icon
		dc.DrawIcon(x, y, m_hIcon);
	}
	else
	{
		CDialog::OnPaint();
	}
}

// The system calls this to obtain the cursor to display while the user drags
//  the minimized window.
HCURSOR CGPSDlg::OnQueryDragIcon()
{
	return (HCURSOR) m_hIcon;
}

double curve2(double a,double e,double B) //计算卯酉圈曲率半径N
{
	double N;
	double n0=a;
	double n1=e*e*n0/2;
	double n2=3*e*e*n1/4;
	double n3=5*e*e*n2/6;
	double n4=7*e*e*n3/8;
    double SINB;
	SINB=sin(B);
	N=n0+n1*pow(SINB,2)+n2*pow(SINB,4)+n3*pow(SINB,6)+n4*pow(SINB,8); 
	return N;
}

double Radian1(double a,double e,double B)  //计算X(自赤道量起的子午线弧长)
{
	double X;
	double m0=a*(1-e*e);
	double m1=3*e*e*m0/2;
	double m2=5*e*e*m1/4;
	double m3=7*e*e*m2/6;
	double m4=9*e*e*m3/8;
	double a0=m0+m1/2+3*m2/8+5*m3/16+35*m4/128;
	double a1=m1/2+m2/2+15*m3/32+7*m4/16;
	double a2=m2/8+3*m3/16+7*m4/32;
	double a3=m3/32+m4/16;
	double a4=m4/128;
	X=a0*B-a1*sin(2*B)/2+a2*sin(4*B)/4-a3*sin(6*B)/6+a4*sin(8*B);
	return X;
}


long double dmstorad(long double dms) //度分秒转换成弧度
{ 
	long double Degree, Minute; 
	long double Second; 
	int Sign; 
	long double rad; 
	if(dms >= 0) 
		Sign = 1; 
	else 
		Sign = -1; 
	dms = fabs(dms); 
	Degree = floor(dms); 
	Minute = floor(fmod(dms * 100.0, 100.0)); 
	Second = fmod(dms * 10000.0, 100.0); 
	rad = Sign * (Degree + Minute / 60.0 + Second / 3600.0) * PI / 180.0; 
	return rad; 
}


long double radtodms(double rad) //弧度转换成度分秒
{ 
	double Degree, Minute; 
	double Second; 
	int Sign; 
	double dms; 
	if(rad >= 0) 
		Sign = 1; 
	else 
		Sign = -1; 
	rad = fabs(rad * 180.0 / PI); 
	Degree = floor(rad); 
	Minute = floor(fmod(rad * 60.0, 60.0)); 
	Second = fmod(rad * 3600.0, 60.0); 
	dms = Sign * (Degree + Minute / 100.0 + Second / 10000.0); 
	return dms; 
} 

void CGPSDlg::OnCalculate() 
{
	// TODO: Add your control notification handler code here
	
	//1.BLH转化XYZ
	UpdateData(TRUE);
	long double a,e;
	a=6378137.0;
	e=sqrtl(0.0066943799013);
    long double m_B1=dmstorad(m_B);
	long double m_L1=dmstorad(m_L);
	long double W = sqrtl(1 - e*e*sinl(m_B1)*sinl(m_B1));
	long double N = a / W;
	m_X=(N + m_H) * cosl(m_B1) * cosl(m_L1);
	m_Y=(N + m_H) * cosl(m_B1) * sinl(m_L1);
	m_Z=(N * (1- e*e) + m_H) * sinl(m_B1);
	
	
	//2.XYZ转化BLH
	//a=6378245.0;
	//e=sqrtl(0.006693421622966);
	long double t0 , p , k , e1_2 , c;	
	e1_2 = e*e / (1 - e*e);
	c = a * sqrtl(1 + e1_2);
	k  = 1 + e1_2;
	
	long double R_2 = m_X*m_X + m_Y*m_Y;
	
	t0 = m_Z / sqrtl(R_2);	
	p  = c * e*e / sqrtl(R_2);	
	const long double precision = PI / 180 / 3600 / 100000000;
	long double tanBmin = tanl(precision);
	
	long double t1,t2;
	long double fx , dfx;
	
	t1 = t0;
	do{		
		fx = t1-t0-p*t1/sqrtl(k+t1*t1);
		dfx = 1 - p*k / (k+t1*t1) / sqrtl(k+t1*t1);
		t2 = t1;
		t1 = t1 - fx / dfx;
	}
	while(fabs(t1 - t2) > tanBmin);
	m_BB = atanl(t1);
	W = sqrtl(1 - e*e*sinl(m_BB)*sinl(m_BB));
	N = a / W;
	m_HH = sqrtl(R_2) / cosl(m_BB) - N;
   	m_LL = acosl(m_X/sqrtl(R_2));
	m_BB =radtodms(m_BB);
	m_LL =radtodms(m_LL);


	//3.高斯投影
	double X,x,y,L0,L=m_LL;
	L0=114;
	double l=(dmstorad(L)-dmstorad(L0))*180*3600/PI;
	double B=dmstorad(m_BB);
	double t=tan(B);
	double n=e*cos(B)/sqrt(1-e*e);
	N=curve2(a,e,B);
	X=Radian1(a,e,B);
	double r=180*3600/PI;
	x=X+N*sin(B)*cos(B)*l*l/(2*r*r)+N*sin(B)*pow(cos(B),3)*(5-t*t+9*n*n+
		4*pow(n,4))/(24*pow(r,4))+N*sin(B)*pow(cos(B),5)*(61-58*t*t+pow(t,4))*pow(l,6)/(720*pow(r,6));
    y=500000+N*cos(B)*l/r+N*pow(cos(B),3)*(1-t*t+n*n)*pow(l,3)/(6*pow(r,3))+
		N*pow(cos(B),5)*(5-18*t*t+pow(t,4)+14*n*n-58*n*n*t*t)*pow(l,5)/(120*pow(r,5));
	m_XX=x;
	m_YY=y;
    m_L0=int(m_LL/3);
	UpdateData(FALSE);

	
}