大地测量_坐标转换_c++代码

#include<iostream> #include<string> #include<sstream> #define pi 3.14159265358979323846 using namespace std; struct dms { int d; int m; double s; }; bool DmstoD(double &L, dms &lon);; bool XYZtoGeo(double &x, double &y, double &z, dms &lon, dms &lat,double &H) { cout << "请输入X(单位为m): " ; cin >> x; cout << "请输入Y(单位为m): " ; cin >> y; cout << "请输入Z(单位为m): " ; cin >> z; cout << "X是 " << x << "°" << endl; cout << "Y是 " << y << "°" << endl; cout << "Z是 " << z << "°" << endl; //WGS-84长半轴a，短半轴b double a = 6378137; double b = 6356752.314200000; double dtr = pi / 180.0; double ee = (a*a - b * b) / (a*a); //计算经度L double tp = sqrt(x*x + y * y); double L = acos(x / tp); if (y < 0) L = 2 * pi - L; //计算纬度B double B0 = atan(z / tp); int i = 1; while (1) //对B进行迭代 { cout << "第 " << i << " 次迭代： " << endl; i += 1; double W = sqrt(1 - ee * sin(B0)*sin(B0)); double N = a / W; //卯酉圈曲率半径 double tp1 = z + N * ee*sin(B0); cout << fixed << "B0 = " << B0 << endl; cout << "W = " << B0 << endl; cout << "N = " << B0 << endl; double B = atan(tp1 / tp); cout << "B = " << B << endl; if (abs(B0 - B) < 1.0e-15) break; cout << "------------------------------" << endl; B0 = B; H = tp / cos(B) - N; } DmstoD(L, lon); DmstoD(B0, lat); cout << "经度是: " << lon.d << "°" << lon.m << "'" << lon.s << "''" << endl; cout << "纬度是: " << lat.d << "°" << lat.m << "'" << lat.s << "''" << endl; cout << "高程是:" << H << "m" << endl; return 1; }