c++空间直角坐标转换东北天坐标代码
时间: 2023-07-09 12:08:05 浏览: 75
以下是C++实现的空间直角坐标转换为东北天坐标的代码:
```cpp
#include <iostream>
#include <cmath>
using namespace std;
struct Point3D {
double x, y, z;
};
struct PointENU {
double e, n, u;
};
const double PI = 3.14159265358979323846;
const double a = 6378137.0; // 长半轴
const double b = 6356752.3142; // 短半轴
const double e = sqrt(pow(a, 2) - pow(b, 2)) / a; // 第一偏心率e
// 弧度制转换为角度制
double rad2deg(double rad) {
return rad * 180.0 / PI;
}
// 角度制转换为弧度制
double deg2rad(double deg) {
return deg * PI / 180.0;
}
// 计算子午圈曲率半径
double radiusOfCurvature(double B) {
double sinB = sin(B);
return a * (1 - pow(e, 2)) / pow(1 - pow(e * sinB, 2), 1.5);
}
// 空间直角坐标转换为大地坐标
Point3D XYZ2BLH(Point3D p) {
Point3D result;
double B, L, H;
double rho, N, sinB;
double cosB, sinL, cosL;
rho = sqrt(pow(p.x, 2) + pow(p.y, 2));
N = radiusOfCurvature(B);
sinB = p.z / sqrt(pow(p.x, 2) + pow(p.y, 2) + pow(p.z, 2));
B = atan(p.z / sqrt(pow(p.x, 2) + pow(p.y, 2)) / (1 - pow(e, 2)));
do {
B = atan((p.z + N * e * e * sinB) / rho);
N = radiusOfCurvature(B);
sinB = (p.z + N * e * e * sinB) / sqrt(pow(p.x, 2) + pow(p.y, 2));
} while (abs(B - atan(p.z / rho / (1 - pow(e, 2))))) > 1E-10;
L = atan(p.y / p.x);
cosB = cos(B);
sinB = sin(B);
sinL = sin(L);
cosL = cos(L);
H = rho / cosB - N;
result.x = rad2deg(B);
result.y = rad2deg(L);
result.z = H;
return result;
}
// 大地坐标转换为东北天坐标
PointENU BLH2ENU(Point3D p, Point3D origin) {
PointENU result;
double B, L, H;
double x0, y0, z0;
double cosB, sinB, cosL, sinL;
B = deg2rad(p.x);
L = deg2rad(p.y);
H = p.z;
x0 = origin.x;
y0 = origin.y;
z0 = origin.z;
cosB = cos(B);
sinB = sin(B);
cosL = cos(L);
sinL = sin(L);
double dx = -(cosB * cosL * x0 + cosB * sinL * y0 + sinB * z0);
double dy = -(-sinL * x0 + cosL * y0);
double dz = -(sinB * cosL * x0 + sinB * sinL * y0 - cosB * z0);
result.e = -sinL * dx + cosL * dy;
result.n = -sinB * cosL * dx - sinB * sinL * dy + cosB * dz;
result.u = cosB * cosL * dx + cosB * sinL * dy + sinB * dz;
return result;
}
int main() {
Point3D point{ 2000000, 3000000, 4000000 }; // 空间直角坐标
Point3D origin{ 31.2304, 121.4737, 10 }; // 原点大地坐标
Point3D p = XYZ2BLH(point); // 空间直角坐标转换为大地坐标
PointENU enu = BLH2ENU(p, origin); // 大地坐标转换为东北天坐标
cout << "E: " << enu.e << endl;
cout << "N: " << enu.n << endl;
cout << "U: " << enu.u << endl;
return 0;
}
```
需要注意的是,由于地球不是完全的椭球体,所以在计算过程中可能存在误差。此外,该代码中假设输入的空间直角坐标是以米为单位的。
相关推荐
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)