大地坐标转大地空间直角坐标，经纬度计算距离，源代码如下：
源代码下载地址：
// CoordinateTransformationDlg.cpp : 实现文件
//#include "stdafx.h"
#include "CoordinateTransformation.h"
#include "CoordinateTransformationDlg.h"
#include "afxdialogex.h"
#include <math.h>
#include <iostream>
using namespace std;#ifdef _DEBUG
#define new DEBUG_NEW
#endif // 数学符号pi
#ifndef PI
#define PI 3.1415926535897932384626433832795
#endif
// 数学符号pi
#ifndef EARTH_RADIUS
#define EARTH_RADIUS 6378137.0
#endif
// 用于应用程序“关于”菜单项的 CAboutDlg 对话框class CAboutDlg : public CDialogEx
{
public:
CAboutDlg();// 对话框数据
enum { IDD = IDD_ABOUTBOX }; protected:
virtual void DoDataExchange(CDataExchange* pDX); // DDX/DDV 支持// 实现
protected:
DECLARE_MESSAGE_MAP()
};CAboutDlg::CAboutDlg() : CDialogEx(CAboutDlg::IDD)
{
}void CAboutDlg::DoDataExchange(CDataExchange* pDX)
{
CDialogEx::DoDataExchange(pDX);
}BEGIN_MESSAGE_MAP(CAboutDlg, CDialogEx)
END_MESSAGE_MAP() // CCoordinateTransformationDlg 对话框
#define PI 3.1415926535897323
double a,b,c,e2,ep2;
CCoordinateTransformationDlg::CCoordinateTransformationDlg(CWnd* pParent /*=NULL*/)
: CDialogEx(CCoordinateTransformationDlg::IDD, pParent)
, m_dflon(0)
, m_dfLat(0)
, m_dfMaxAxis(6378137)
, m_dfMinAxis(6356752.3142)
, m_dfCanxinX(0)
, m_dfCanxinY(0)
, m_dfCanxinZ(0)
, m_dfPingmianRoteAngle(105.6)
, m_dfCanxinPMRoteX(0)
, m_dfCanxinPMRoteZ(0)
, m_dfCanxinPMRoteY(0)
, m_dfAngle(0)
, m_dfMinite(0)
, m_dfSeconde(0)
, m_dfBL(0)
, m_dfTuoqiuGao(0)
, m_dfP1Longitude(0)
, m_dfP2Longitude(0)
, m_dfP1Latitude(0)
, m_dfP2Latitude(0)
, m_dfP1P2Distance(0)
{
m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);
}void CCoordinateTransformationDlg::DoDataExchange(CDataExchange* pDX)
{
CDialogEx::DoDataExchange(pDX);
DDX_Text(pDX, IDC_EDIT1_LONTITUDE, m_dflon);
DDX_Text(pDX, IDC_EDIT2_LATITUDE, m_dfLat);
DDX_Text(pDX, IDC_EDIT3_MAJORSEMIAXIS, m_dfMaxAxis);
DDX_Text(pDX, IDC_EDIT4_SEMI_MINOR_AXIS, m_dfMinAxis);
DDX_Text(pDX, IDC_EDIT5_CANXIN_X, m_dfCanxinX);
DDX_Text(pDX, IDC_EDIT6_CANXIN_Y, m_dfCanxinY);
DDX_Text(pDX, IDC_EDIT7_CANXIN_Z, m_dfCanxinZ);
DDX_Text(pDX, IDC_EDIT8_PINGMIAN_MOVE_ANGLE, m_dfPingmianRoteAngle);
DDX_Text(pDX, IDC_EDIT9_PINGMIANMOVE_CANXIN_X, m_dfCanxinPMRoteX);
DDX_Text(pDX, IDC_EDIT10_PINGMIAN_MOVE_CANXINZ, m_dfCanxinPMRoteZ);
DDX_Text(pDX, IDC_EDIT11_PINGMIAN_MOVE_CANXINY, m_dfCanxinPMRoteY);
DDX_Text(pDX, IDC_EDIT12_DU, m_dfAngle);
DDX_Text(pDX, IDC_EDIT13_FEN, m_dfMinite);
DDX_Text(pDX, IDC_EDIT14_MIAO, m_dfSeconde);
DDX_Text(pDX, IDC_EDIT15_NEW_DU, m_dfBL);
DDX_Text(pDX, IDC_EDIT16_TUOQIU_GAO, m_dfTuoqiuGao);
DDX_Control(pDX, IDC_COMBO1_TUOQIU_NAME, m_ComboTuoQiuName);
DDX_Text(pDX, IDC_EDIT1_POINT1_LONG, m_dfP1Longitude);
DDX_Text(pDX, IDC_EDIT3_POINT2_LONG, m_dfP2Longitude);
DDX_Text(pDX, IDC_EDIT2_POINT1_LAT, m_dfP1Latitude);
DDX_Text(pDX, IDC_EDIT4_POINT2_LAT, m_dfP2Latitude);
DDX_Text(pDX, IDC_EDIT5_LATLONG_DISTANCE, m_dfP1P2Distance);
}BEGIN_MESSAGE_MAP(CCoordinateTransformationDlg, CDialogEx)
ON_WM_SYSCOMMAND()
ON_WM_PAINT()
ON_WM_QUERYDRAGICON()
ON_BN_CLICKED(IDC_BUTTON1_BL_TO_CANXIN_XYZ, &CCoordinateTransformationDlg::OnBnClickedButton1BlToCanxinXyz)
ON_BN_CLICKED(IDC_BUTTON2_CANXIN_XYZ_TO_CANXIN_PINGMIANMOVE_XYZ, &CCoordinateTransformationDlg::OnBnClickedButton2CanxinXyzToCanxinPingmianmoveXyz)
ON_BN_CLICKED(IDC_BUTTON3_DUFENMIAO_TO_DU, &CCoordinateTransformationDlg::OnBnClickedButton3DufenmiaoToDu)
ON_BN_CLICKED(IDC_BUTTON4_DU_TO_DUFENMIAO, &CCoordinateTransformationDlg::OnBnClickedButton4DuToDufenmiao)
ON_CBN_SELCHANGE(IDC_COMBO1_TUOQIU_NAME, &CCoordinateTransformationDlg::OnCbnSelchangeCombo1TuoqiuName)
ON_BN_CLICKED(IDC_BUTTON1_COMPUTE_LATLONG_DISTAND, &CCoordinateTransformationDlg::OnBnClickedButton1ComputeLatlongDistand)
ON_BN_CLICKED(IDC_BTN_XYZ_TO_BL, &CCoordinateTransformationDlg::OnBnClickedBtnXyzToBl)
END_MESSAGE_MAP() // CCoordinateTransformationDlg 消息处理程序
BOOL CCoordinateTransformationDlg::OnInitDialog()
{
CDialogEx::OnInitDialog(); // 将“关于...”菜单项添加到系统菜单中。
// IDM_ABOUTBOX 必须在系统命令范围内。
ASSERT((IDM_ABOUTBOX & 0xFFF0) == IDM_ABOUTBOX);
ASSERT(IDM_ABOUTBOX < 0xF000); CMenu* pSysMenu = GetSystemMenu(FALSE);
if (pSysMenu != NULL)
{
BOOL bNameValid;
CString strAboutMenu;
bNameValid = strAboutMenu.LoadString(IDS_ABOUTBOX);
ASSERT(bNameValid);
if (!strAboutMenu.IsEmpty())
{
pSysMenu->AppendMenu(MF_SEPARATOR);
pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);
}
} // 设置此对话框的图标。当应用程序主窗口不是对话框时，框架将自动
// 执行此操作
SetIcon(m_hIcon, TRUE); // 设置大图标
SetIcon(m_hIcon, FALSE); // 设置小图标 // TODO: 在此添加额外的初始化代码
// printf("请选择椭球参数(输入椭球序号):\n");
// printf("1.克拉索夫斯基椭球参数\n");
// printf("2.IUGG_1975椭球参数\n");
// printf("3.CGCS_2000椭球参数\n");
// printf("0.其他椭球参数(自行输入)\n");
// scanf("%d",&n); m_ComboTuoQiuName.AddString(_T("WGS84 椭球"));
m_ComboTuoQiuName.AddString(_T("克拉索夫斯基椭球"));
m_ComboTuoQiuName.AddString(_T("IUGG_1975椭球"));
m_ComboTuoQiuName.AddString(_T("CGCS_2000椭球"));
m_ComboTuoQiuName.AddString(_T("其他椭球参数(自行输入)"));
m_ComboTuoQiuName.SetCurSel(0); CWnd* wnd = GetDlgItem(IDC_EDIT3_MAJORSEMIAXIS);
wnd->EnableWindow(FALSE);
wnd = GetDlgItem(IDC_EDIT4_SEMI_MINOR_AXIS);
wnd->EnableWindow(FALSE);
int n = m_ComboTuoQiuName.GetCurSel();
switch(n)
{
case 0:
{
a=6378137;b=6356752.3142;
c=a*a/b;
ep2=(a*a-b*b)/(b*b);
e2=(a*a-b*b)/(a*a);
break;
}
case 1:
{
a=6378245.0;b=6356863.0188;c=6399698.9018;e2=0.00669342162297;ep2=0.00673852541468;break;
}
case 2:a=6378140.0;b=6356755.2882;c=6399596.6520;e2=0.00669438499959;ep2=0.00673950181947;break;
case 3:a=6378137.0;b=6356752.3141;c=6399593.6259;e2=0.00669438002290;ep2=0.00673949677547;break;
default:
break;
}
return TRUE; // 除非将焦点设置到控件，否则返回 TRUE
}void CCoordinateTransformationDlg::OnSysCommand(UINT nID, LPARAM lParam)
{
if ((nID & 0xFFF0) == IDM_ABOUTBOX)
{
CAboutDlg dlgAbout;
dlgAbout.DoModal();
}
else
{
CDialogEx::OnSysCommand(nID, lParam);
}
}// 如果向对话框添加最小化按钮，则需要下面的代码
// 来绘制该图标。对于使用文档/视图模型的 MFC 应用程序，
// 这将由框架自动完成。void CCoordinateTransformationDlg::OnPaint()
{
if (IsIconic())
{
CPaintDC dc(this); // 用于绘制的设备上下文 SendMessage(WM_ICONERASEBKGND, reinterpret_cast<WPARAM>(dc.GetSafeHdc()), 0);
// 使图标在工作区矩形中居中
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; // 绘制图标
dc.DrawIcon(x, y, m_hIcon);
}
else
{
CDialogEx::OnPaint();
}
}//当用户拖动最小化窗口时系统调用此函数取得光标
//显示。
HCURSOR CCoordinateTransformationDlg::OnQueryDragIcon()
{
return static_cast<HCURSOR>(m_hIcon);
}double RAD(double d,double f,double m)
{
double e;
double sign=(d<0.0)?-1.0:1.0;
if(d==0)
{
sign=(f<0.0)?-1.0:1.0;
if(f==0)
{
sign=(m<0.0)?-1.0:1.0;
}
}
if(d<0)
d=d*(-1.0);
if(f<0)
f=f*(-1.0);
if(m<0)
m=m*(-1.0); e=sign*(d*3600+f*60+m)*PI/(3600*180);
return e;
}
double RAD(double du)
{ double d = int(du);
double dfTmp = du - d;
double dfMiniteTmp = (dfTmp * 60);
double f = (int) dfMiniteTmp;
double m = (dfMiniteTmp - f)*60; double e;
double sign=(d<0.0)?-1.0:1.0;
if(d==0)
{
sign=(f<0.0)?-1.0:1.0;
if(f==0)
{
sign=(m<0.0)?-1.0:1.0;
}
}
if(d<0)
d=d*(-1.0);
if(f<0)
f=f*(-1.0);
if(m<0)
m=m*(-1.0); e=sign*(d*3600+f*60+m)*PI/(3600*180);
return e;
}
void RBD(double hd)
{
int t;
int d,f;
double m;
double sign=(hd<0.0)?-1.0:1.0;
if(hd<0)
hd=fabs(hd);
hd=hd*3600*180/PI;
t=int(hd/3600);
d=sign*t;
hd=hd-t*3600;
f=int(hd/60);
m=hd-f*60;
//printf("%d'%d'%lf'\n",d,f,m);
}
void BLH_XYZ()
{
double B,L,H,N,W;
double d,f,m;
double X,Y,Z;
printf(" 请输入大地坐标（输入格式为角度(例如：30'40'50')）：\n");
printf(" 大地经度 L=");
scanf("%lf'%lf'%lf'",&d,&f,&m);
L=RAD(d,f,m);
printf(" 大地纬度 B=");
scanf("%lf'%lf'%lf'",&d,&f,&m);
B=RAD(d,f,m);
printf(" 大地高 H=");
scanf("%lf",&H); W=sqrt(1-e2*sin(B)*sin(B));
N=a/W;
X=(N+H)*cos(B)*cos(L);
Y=(N+H)*cos(B)*sin(L);
Z=(N*(1-e2)+H)*sin(B); printf("\n\n 转换后得到大地空间直角坐标为：\n\n");
printf("X=%lf\nY=%lf\nZ=%lf\n\n",X,Y,Z);
}
void XYZ_BLH()
{
double B,L,H,N,W;
double X,Y,Z;
double tgB0,tgB1; printf(" 请输入大地空间直角坐标：\n");
printf(" X=");
scanf("%lf",&X);
printf(" Y=");
scanf("%lf",&Y);
printf(" Z=");
scanf("%lf",&Z); printf("\n\n 转换后得到大地坐标为：\n\n");
L=atan(Y/X);
printf(" 大地经度为： L=");
RBD(L);
printf("\n"); tgB0=Z/sqrt(X*X+Y*Y);
tgB1=(1/sqrt(X*X+Y*Y))*(Z+a*e2*tgB0/sqrt(1+tgB0*tgB0-e2*tgB0*tgB0)); while(fabs(tgB0-tgB1)>5*pow(10.0,-10))
{
tgB0=tgB1;
tgB1=(1/sqrt(X*X+Y*Y))*(Z+a*e2*tgB0/sqrt(1+tgB0*tgB0-e2*tgB0*tgB0));
}
B=atan(tgB1);
printf(" 大地纬度为：B=");
RBD(B);
printf("\n");
W=sqrt(1-e2*sin(B)*sin(B));
N=a/W;
H=sqrt(X*X+Y*Y)/cos(B)-N;
printf(" 大地高为：H=%lf\n\n",H);
}
void B_ZS()
{
double L1,B1,A1,s,d,f,mi;
double u1,u2,m,M,k2,alfa,bt,r,kp2,alfap,btp,rp;
double sgm0,sgm1,lmd,lmd1,lmd2,A2,B2,l,L2; printf("请输入已知点的大地坐标（输入格式为角度(例如：30'40'50')，下同）：\nL1=");
scanf("%lf'%lf'%lf'",&d,&f,&mi);
L1=RAD(d,f,mi);
printf("\nB1=");
scanf("%lf'%lf'%lf'",&d,&f,&mi);
B1=RAD(d,f,mi);
printf("请输入大地方位角：\nA1=");
scanf("%lf'%lf'%lf'",&d,&f,&mi);
A1=RAD(d,f,mi);
printf("请输入该点至另一点的大地线长：\ns=");
scanf("%lf",&s); u1=atan(sqrt(1-e2)*tan(B1));
m=asin(cos(u1)*sin(A1));
M=atan(tan(u1)/cos(A1));
m=(m>0)?m:m+2*PI;
M=(M>0)?M:M+PI;
k2=ep2*cos(m)*cos(m);
alfa=(1-k2/4+7*k2*k2/64-15*k2*k2*k2/256)/b;
bt=k2/4-k2*k2/8+37*k2*k2*k2/512;
r=k2*k2/128-k2*k2*k2/128; sgm0=alfa*s;
sgm1=alfa*s+bt*sin(sgm0)*cos(2*M+sgm0)+r*sin(2*sgm0)*cos(4*M+2*sgm0);
while(fabs(sgm0-sgm1)>2.8*PI/180*pow(10.0,-7))
{
sgm0=sgm1;
sgm1=alfa*s+bt*sin(sgm0)*cos(2*M+sgm0)+r*sin(2*sgm0)*cos(4*M+2*sgm0);
}
sgm0=sgm1; A2=atan(tan(m)/cos(M+sgm0));
A2=(A2>0)?A2:A2+PI;
A2=(A1>PI)?A2:A2+PI;
u2=atan(-cos(A2)*tan(M+sgm0));
lmd1=atan(sin(u1)*tan(A1));
lmd1=(lmd1>0)?lmd1:lmd1+PI;
lmd1=(m<PI)?lmd1:lmd1+PI;
lmd2=atan(sin(u2)*tan(A2));
lmd2=(lmd2>0)?lmd2:lmd2+PI;
lmd2=(m<PI)?(((M+sgm0)<PI)?lmd2:lmd2+PI):(((M+sgm0)>PI)?lmd2:lmd2+PI);
lmd=lmd2-lmd1;
B2=atan(sqrt(1+ep2)*tan(u2)); kp2=e2*cos(m)*cos(m);
alfap=(e2/2+e2*e2/8+e2*e2*e2/16)-e2/16*(1+e2)*kp2+3*e2*kp2*kp2/128;
btp=e2*(1+e2)*kp2/16-e2*kp2*kp2/32;
rp=e2*kp2*kp2/256; l=lmd-sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos(2*M+sgm0)+rp*sin(2*sgm0)*cos(4*M+2*sgm0));
L2=L1+l; printf("\n\n得到另一点的大地坐标和大地线在该点的大地方位角为：\n\n");
printf("L2=");
RBD(L2);printf("\n");
printf("B2=");
RBD(B2);printf("\n");
printf("A2=");
RBD(A2);
printf("\n");
}
void B_FS()
{
double L1,B1,L2,B2,s,A1,A2,du,f,mi,m0,m,M;
double l,u1,u2,alfa,bt,r,lmd0,dit_lmd,lmd,sgm,dit_sgm,sgm0,sgm1,alfap,btp,rp,k2,kp2; printf("请输入第一个点大地坐标（输入格式为角度(例如：30'40'50')，下同）：\n大地经度 L1=");
scanf("%lf'%lf'%lf'",&du,&f,&mi);
L1=RAD(du,f,mi);
printf("大地纬度 B1=");
scanf("%lf'%lf'%lf'",&du,&f,&mi);
B1=RAD(du,f,mi);
printf("\n请输入第二个点大地坐标：\n大地经度：L2=");
scanf("%lf'%lf'%lf'",&du,&f,&mi);
L2=RAD(du,f,mi);
printf("大地纬度：B2=");
scanf("%lf'%lf'%lf'",&du,&f,&mi);
B2=RAD(du,f,mi); l=L2-L1;
u1=atan(sqrt(1-e2)*tan(B1));
u2=atan(sqrt(1-e2)*tan(B2));
sgm0=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l));
m0=asin(cos(u1)*cos(u2)*sin(l)/sin(sgm0));
dit_lmd=0.003351831*sgm0*sin(m0);
lmd0=l+dit_lmd;
dit_sgm=sin(m0)*dit_lmd;
sgm1=sgm0+dit_sgm;
m=asin(cos(u1)*cos(u2)*sin(lmd0)/sin(sgm1));
A1=atan(sin(lmd0)/(cos(u1)*tan(u2)-sin(u1)*cos(lmd0)));
A1=(A1>0)?A1:A1+PI;
A1=(m>0)?A1:A1+PI; M=atan(sin(u1)*tan(A1)/sin(m));
M=(M>0)?M:M+PI;
k2=ep2*cos(m)*cos(m);
alfa=(1-k2/4+7*k2*k2/64-15*k2*k2*k2/256)/b;
bt=k2/4-k2*k2/8+37*k2*k2*k2/512;
r=k2*k2/128-k2*k2*k2/128; kp2=e2*cos(m)*cos(m);
alfap=(e2/2+e2*e2/8+e2*e2*e2/16)-e2/16*(1+e2)*kp2+3*e2*kp2*kp2/128;
btp=e2*(1+e2)*kp2/16-e2*kp2*kp2/32;
rp=e2*kp2*kp2/256; sgm0=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l));
sgm1=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l+sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos(2*M+sgm0))));
while(fabs(sgm0-sgm1)>1*PI/180*pow(10.0,-8))
{
sgm0=sgm1;
sgm1=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l+sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos(2*M+sgm0))));
}
sgm=sgm1;
lmd=l+sin(m)*(alfap*sgm+btp*sin(sgm)*cos(2*M+sgm)); s=(sgm-bt*sin(sgm)*cos(2*M+sgm)-r*sin(2*sgm)*cos(4*M+2*sgm))/alfa;
A1=atan(sin(lmd)/(cos(u1)*tan(u2)-sin(u1)*cos(lmd)));
A1=(A1>0)?A1:A1+PI;
A1=(m>0)?A1:A1+PI;
A2=atan(sin(lmd)/(sin(u2)*cos(lmd)-tan(u1)*cos(u2)));
A2=(A2>0)?A2:A2+PI;
A2=(m<0)?A2:A2+PI; printf("\n\n得到两点间大地线长S和大地正反方位角A1、A2如下：\n\n");
printf("s=%lf\n",s);
printf("A1=");
RBD(A1);printf("\n");
printf("A2=");
RBD(A2);printf("\n");
}
void GUS_ZS()
{
double B,L,x3,x6,y3,y6,Y3,Y6,du,f,mi,X,N,n,t;
double At,Bt,Ct,Dt,m3,m6,l3,l6,W,L03,L06;
int DH3,DH6; printf("请输入大地坐标（输入格式为角度(例如：30'40'50')）：\n大地经度 L=");
scanf("%lf'%lf'%lf'",&du,&f,&mi);
L=RAD(du,f,mi);
printf("\n大地纬度 B=");
scanf("%lf'%lf'%lf'",&du,&f,&mi);
B=RAD(du,f,mi); At=1+3*e2/4+45*e2*e2/64+175*e2*e2*e2/256;
Bt=3*e2/4+15*e2*e2/16+525*e2*e2*e2/512;
Ct=15*e2*e2/64+105*e2*e2*e2/256;
Dt=35*e2*e2*e2/512;
X=a*(1-e2)*(At*B-Bt*sin(2*B)/2+Ct*sin(4*B)/4-Dt*sin(6*B)/6); W=sqrt(1-e2*sin(B)*sin(B));
N=a/W;
n=sqrt(ep2)*cos(B);
t=tan(B); DH3=(L-(1.5*PI/180))/(3*PI/180)+1;
DH6=L/(6*PI/180)+1;
L03=DH3*(3*PI/180);
L06=DH6*(6*PI/180)-(3*PI/180);
l3=L-L03;
l6=L-L06;
m3=cos(B)*l3;
m6=cos(B)*l6; x3=X+N*t*(m3*m3/2+(5-t*t+9*n*n+4*n*n*n*n)*m3*m3*m3*m3/24+(61-58*t*t+t*t*t*t)*m3*m3*m3*m3*m3*m3/720);
x6=X+N*t*(m6*m6/2+(5-t*t+9*n*n+4*n*n*n*n)*m6*m6*m6*m6/24+(61-58*t*t+t*t*t*t)*m6*m6*m6*m6*m6*m6/720);
y3=N*(m3+(1-t*t+n*n)*m3*m3*m3/6+(5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t)*m3*m3*m3*m3*m3/120);
y6=N*(m6+(1-t*t+n*n)*m6*m6*m6/6+(5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t)*m6*m6*m6*m6*m6/120);
Y3=DH3*1000000+500000+y3;
Y6=DH6*1000000+500000+y6; printf("\n\n 得到的高斯平面坐标为:\n\n");
printf(" 对于3度带:\n 纵坐标x=%.3lf\n 横坐标y=%.3lf(通用坐标Y=%.3lf)\n\n",x3,y3,Y3);
printf(" 对于6度带:\n 纵坐标x=%.3lf\n 横坐标y=%.3lf(通用坐标Y=%.3lf)\n\n",x6,y6,Y6);
}
void GUS_FS()
{
double x,y,Y,B,B0,B1,Bf,Vf,tf,Nf,nf,L,At,Bt,Ct,Dt,L3,L6;
long DH; printf(" 请输入高斯平面坐标：\n\n");
printf(" 纵坐标 X=");
scanf("%lf",&x);printf("\n");
printf(" 自然坐标 y=");
scanf("%lf",&y);printf("\n");
printf(" 通用坐标 Y=");
scanf("%lf",&Y);printf("\n"); At=1+3*e2/4+45*e2*e2/64+175*e2*e2*e2/256;
Bt=3*e2/4+15*e2*e2/16+525*e2*e2*e2/512;
Ct=15*e2*e2/64+105*e2*e2*e2/256;
Dt=35*e2*e2*e2/512;
B0=x/(a*(1-e2)*At);
B1=(x-a*(1-e2)*(-Bt*sin(2*B0)/2+Ct*sin(4*B0)/4-Dt*sin(6*B0)/6))/(a*(1-e2)*At);
while(fabs(B1-B0)>1*pow(10.0,-8))
{
B0=B1;
B1=(x-a*(1-e2)*(-Bt*sin(2*B0)/2+Ct*sin(4*B0)/4-Dt*sin(6*B0)/6))/(a*(1-e2)*At);
}
Bf=B1; nf=sqrt(ep2)*cos(Bf);
tf=tan(Bf);
Vf=sqrt(1+ep2*cos(Bf)*cos(Bf));
Nf=c/Vf;
B=Bf-Vf*Vf*tf/2*((y/Nf)*(y/Nf)-(5+3*tf*tf+nf*nf-9*nf*nf*tf*tf)*pow((y/Nf),4)/12+(61+90*tf*tf+45*tf*tf)*pow((y/Nf),6)/360);
L=((y/Nf)-(1+2*tf*tf+nf*nf)*(y/Nf)*(y/Nf)*(y/Nf)/6+(5+28*tf*tf+24*pow(tf,4)+6*nf*nf+8*nf*nf*tf*tf)*pow((y/Nf),5)/120)/cos(Bf); DH=Y/1000000;
L3=3*PI/180*double(DH)+L;
L6=6*PI/180*double(DH)-3*PI/180+L; printf("\n\n 得到的大地坐标为:\n\n");
printf(" 大地纬度 B=");
RBD(B);printf("\n");
printf(" 若为6度带，大地经度 L=");
RBD(L6);printf("\n");
printf(" 若为3度带，大地经度 L=");
RBD(L3);printf("\n");
}void CCoordinateTransformationDlg::OnBnClickedButton1BlToCanxinXyz()
{
//经纬度转参心空间直角坐标
UpdateData(TRUE);
// double dfe = sqrt(m_dfMaxAxis*m_dfMaxAxis - m_dfMinAxis*m_dfMinAxis)/m_dfMaxAxis;
// double dfW = sqrt(1-dfe*dfe*sin(m_dfLat)*sin(m_dfLat));
// double dfN = m_dfMaxAxis / dfW;
//
// m_dfCanxinX = (dfN + m_dfTuoqiuGao) * cos(m_dfLat) * cos(m_dflon);
// m_dfCanxinY = (dfN + m_dfTuoqiuGao) * cos(m_dfLat) * sin(m_dflon);
// m_dfCanxinZ = (dfN*(1-dfe*dfe)+ m_dfTuoqiuGao) * sin(m_dfLat); //
if (m_dfMaxAxis < 100 || m_dfMinAxis < 100)
{
AfxMessageBox(_T("请输入规范的椭球参数！"));
return;
} a=m_dfMaxAxis;b=m_dfMinAxis;
c=a*a/b;
ep2=(a*a-b*b)/(b*b);
e2=(a*a-b*b)/(a*a); double B,L,H,N,W;
double d,f,m;
double X,Y,Z;
L=RAD(m_dflon); B=RAD(m_dfLat);
//printf(" 大地高 H=");
H = m_dfTuoqiuGao; W=sqrt(1-e2*sin(B)*sin(B));
N=a/W;
X=(N+H)*cos(B)*cos(L);
Y=(N+H)*cos(B)*sin(L);
Z=(N*(1-e2)+H)*sin(B); m_dfCanxinX = X;
m_dfCanxinY = Y;
m_dfCanxinZ = Z;
UpdateData(FALSE);
} void CCoordinateTransformationDlg::OnBnClickedButton2CanxinXyzToCanxinPingmianmoveXyz()
{
UpdateData(TRUE);
if (m_dfMaxAxis < 100 || m_dfMinAxis < 100)
{
AfxMessageBox(_T("请输入规范的椭球参数！"));
return;
}
if (m_dfPingmianRoteAngle == 0)
{
m_dfCanxinPMRoteZ = m_dfCanxinX;
m_dfCanxinPMRoteX = m_dfCanxinY;
m_dfCanxinPMRoteY = m_dfCanxinZ;
}
else
{
double dfX0 = 0, dfY0 = 0;
// double dfHudu = m_dfPingmianRoteAngle;
double dfHudu = RAD(-m_dfPingmianRoteAngle);
// m_dfCanxinPMRoteX = (m_dfCanxinX - dfX0) * cos(dfHudu) - (m_dfCanxinY - dfY0) * sin(dfHudu) + dfX0;
// m_dfCanxinPMRoteZ = (m_dfCanxinX - dfX0) * sin(dfHudu) + (m_dfCanxinY - dfY0) * cos(dfHudu) + dfY0;
// m_dfCanxinPMRoteY = m_dfCanxinZ; m_dfCanxinPMRoteZ = (m_dfCanxinY * cos(dfHudu) - m_dfCanxinX * sin(dfHudu))/(cos(dfHudu)*cos(dfHudu) - sin(dfHudu)*sin(dfHudu));
m_dfCanxinPMRoteX = (m_dfCanxinY - m_dfCanxinPMRoteZ*cos(dfHudu))/sin(dfHudu);
m_dfCanxinPMRoteY = m_dfCanxinZ;
}
UpdateData(FALSE);
} void CCoordinateTransformationDlg::OnBnClickedButton3DufenmiaoToDu()
{
UpdateData(TRUE);
m_dfBL = m_dfAngle + m_dfMinite/60 + m_dfSeconde/3600;
UpdateData(FALSE);
} void CCoordinateTransformationDlg::OnBnClickedButton4DuToDufenmiao()
{
UpdateData(TRUE);
m_dfAngle = int(m_dfBL);
double dfTmp = m_dfBL - m_dfAngle;
double dfMiniteTmp = (dfTmp * 60);
m_dfMinite = (int) dfMiniteTmp;
m_dfSeconde = (dfMiniteTmp - m_dfMinite)*60; UpdateData(FALSE);
} void CCoordinateTransformationDlg::OnCbnSelchangeCombo1TuoqiuName()
{
UpdateData(TRUE);
CWnd* wnd =NULL;
int n = m_ComboTuoQiuName.GetCurSel();
switch(n)
{
case 0:
{
a=6378137;b=6356752.3142;
c=a*a/b;
ep2=(a*a-b*b)/(b*b);
e2=(a*a-b*b)/(a*a);
m_dfMaxAxis = a;
m_dfMinAxis = b;
wnd = GetDlgItem(IDC_EDIT3_MAJORSEMIAXIS);
wnd->EnableWindow(FALSE);
wnd = GetDlgItem(IDC_EDIT4_SEMI_MINOR_AXIS);
wnd->EnableWindow(FALSE);
break;
}
case 1:
{
a=6378245.0;b=6356863.0188;c=6399698.9018;e2=0.00669342162297;ep2=0.00673852541468;
m_dfMaxAxis = a;
m_dfMinAxis = b;
wnd = GetDlgItem(IDC_EDIT3_MAJORSEMIAXIS);
wnd->EnableWindow(FALSE);
wnd = GetDlgItem(IDC_EDIT4_SEMI_MINOR_AXIS);
wnd->EnableWindow(FALSE);
break;
}
case 2:
{
a=6378140.0;b=6356755.2882;c=6399596.6520;e2=0.00669438499959;ep2=0.00673950181947;
m_dfMaxAxis = a;
m_dfMinAxis = b;
wnd = GetDlgItem(IDC_EDIT3_MAJORSEMIAXIS);
wnd->EnableWindow(FALSE);
wnd = GetDlgItem(IDC_EDIT4_SEMI_MINOR_AXIS);
wnd->EnableWindow(FALSE);
break;
}
case 3:
{
a=6378137.0;b=6356752.3141;c=6399593.6259;e2=0.00669438002290;ep2=0.00673949677547;
m_dfMaxAxis = a;
m_dfMinAxis = b;
wnd = GetDlgItem(IDC_EDIT3_MAJORSEMIAXIS);
wnd->EnableWindow(FALSE);
wnd = GetDlgItem(IDC_EDIT4_SEMI_MINOR_AXIS);
wnd->EnableWindow(FALSE);
break;
}
default:
{
wnd = GetDlgItem(IDC_EDIT3_MAJORSEMIAXIS);
wnd->EnableWindow(TRUE);
wnd = GetDlgItem(IDC_EDIT4_SEMI_MINOR_AXIS);
wnd->EnableWindow(TRUE);
m_dfMaxAxis = 0;m_dfMinAxis = 0;
break;
}
}
UpdateData(FALSE);
} static inline double rad( double degree )
{
return PI * degree / 180.0;
}
static inline double rad2degree( double rad )
{
return rad * 180.0 / PI;
}
static inline double haverSin(double x)
{
double v = sin(x / 2.0);
return v * v;
}double GetRadioDistance(double lon1, double lat1, double lon2, double lat2)
{
if (lon1 < -180 || lon1 > 180)
{
return 0;
}
if (lon2 < -180 || lon2 > 180)
{
return 0;
}
if (lat1 < -90 || lat1 > 90)
{
return 0;
}
if (lat2 < -90 || lat2 > 90)
{
return 0;
}
double radlon1 = rad(lon1);
double radlat1 = rad(lat1);
double radlon2 = rad(lon2);
double radlat2 = rad(lat2); double b = fabs(radlat1 - radlat2);
double a = fabs(radlon1 - radlon2); double h = haverSin(b) + cos(radlat1)*cos(radlat2)*haverSin(a);
h = sqrt(h);
h = asin(h);
double distance = 2 * EARTH_RADIUS * h;
return distance;
}double GetLength(double x1, double y1, double x2, double y2)
{
double dfLength = 0.0;
double cx = x1 - x2;
double cy = y1 - y2;
dfLength = sqrt(cx*cx + cy*cy); return dfLength;
}
void CCoordinateTransformationDlg::OnBnClickedButton1ComputeLatlongDistand()
{
UpdateData(TRUE);
bool bLatlong = true;
if (m_dfP1Longitude < -180 || m_dfP1Longitude > 180)
{
bLatlong = 0;
}
if (m_dfP2Longitude < -180 || m_dfP2Longitude > 180)
{
bLatlong = 0;
}
if (m_dfP1Latitude < -90 || m_dfP1Latitude > 90)
{
bLatlong = 0;
}
if (m_dfP2Latitude < -90 || m_dfP2Latitude > 90)
{
bLatlong = 0;
}
if (!bLatlong)
{
AfxMessageBox(_T("兄台，这哪里是经度啊？分明是想忽悠我嘛"));
return;
} m_dfP1P2Distance = GetRadioDistance(m_dfP1Longitude, m_dfP1Latitude, m_dfP2Longitude, m_dfP2Latitude);
UpdateData(FALSE);
} void CCoordinateTransformationDlg::OnBnClickedBtnXyzToBl()
{ //经纬度转参心空间直角坐标
UpdateData(TRUE); if (m_dfMaxAxis < 100 || m_dfMinAxis < 100)
{
AfxMessageBox(_T("请输入规范的椭球参数！"));
return;
}
double B,L,H,N,W;
double X = m_dfCanxinX,Y=m_dfCanxinY,Z=m_dfCanxinZ;
double tgB0,tgB1;
L=atan(Y/X); m_dflon = rad2degree(L);
printf(" 大地经度为： L=");
RBD(L);
printf("\n"); tgB0=Z/sqrt(X*X+Y*Y);
tgB1=(1/sqrt(X*X+Y*Y))*(Z+a*e2*tgB0/sqrt(1+tgB0*tgB0-e2*tgB0*tgB0)); while(fabs(tgB0-tgB1)>5*pow(10.0,-10))
{
tgB0=tgB1;
tgB1=(1/sqrt(X*X+Y*Y))*(Z+a*e2*tgB0/sqrt(1+tgB0*tgB0-e2*tgB0*tgB0));
}
B=atan(tgB1);
m_dfLat = rad2degree(B);
//printf(" 大地纬度为：B=");
RBD(B);
//printf("\n");
W=sqrt(1-e2*sin(B)*sin(B));
N=a/W;
H=sqrt(X*X+Y*Y)/cos(B)-N;
m_dfTuoqiuGao = H;
//printf(" 大地高为：H=%lf\n\n",H);
UpdateData(FALSE);
}
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
