大地坐标转大地空间直角坐标，经纬度计算距离，源代码如下：
源代码下载地址：
https://download.csdn.net/download/c_bright/11203522
// 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); }
