直角坐标系转化为大地坐标系的C++代码

大地坐标系和直角坐标系之间的转换需要涉及到椭球体的参数，例如椭球体的长半轴、短半轴、扁率等等。这里提供一个可以用于WGS84椭球体（GPS使用的椭球体）的直角坐标系转化为大地坐标系的C++代码。需要注意的是，这里的大地坐标系采用的是经度、纬度和高度。

#include <cmath>

const double a = 6378137.0;  // 椭球体长半轴
const double b = 6356752.314245;  // 椭球体短半轴
const double f = (a - b) / a;  // 扁率
const double e_sq = f * (2 - f);  // 第一偏心率的平方

struct ECEFPoint {
    double x;
    double y;
    double z;
};

struct GeodeticPoint {
    double latitude;
    double longitude;
    double altitude;
};

ECEFPoint toECEF(double latitude, double longitude, double altitude) {
    double sin_lat = std::sin(latitude);
    double cos_lat = std::cos(latitude);
    double sin_lon = std::sin(longitude);
    double cos_lon = std::cos(longitude);
    double N = a / std::sqrt(1 - e_sq * sin_lat * sin_lat);

    ECEFPoint point;
    point.x = (N + altitude) * cos_lat * cos_lon;
    point.y = (N + altitude) * cos_lat * sin_lon;
    point.z = ((1 - e_sq) * N + altitude) * sin_lat;

    return point;
}

GeodeticPoint toGeodetic(double x, double y, double z) {
    double p = std::sqrt(x * x + y * y);
    double theta = std::atan2(z * a, p * b);

    double sin_theta = std::sin(theta);
    double cos_theta = std::cos(theta);

    double latitude = std::atan2(z + e_sq * b * sin_theta * sin_theta * sin_theta,
                                 p - e_sq * a * cos_theta * cos_theta * cos_theta);
    double longitude = std::atan2(y, x);
    double sin_lat = std::sin(latitude);
    double N = a / std::sqrt(1 - e_sq * sin_lat * sin_lat);
    double altitude = p / std::cos(latitude) - N;

    GeodeticPoint point;
    point.latitude = latitude;
    point.longitude = longitude;
    point.altitude = altitude;

    return point;
}

使用时可以按照以下方式进行调用：

double latitude = 39.9087;  // 纬度
double longitude = 116.3975;  // 经度
double altitude = 50;  // 高度
ECEFPoint ecef_point = toECEF(latitude, longitude, altitude);
GeodeticPoint geodetic_point = toGeodetic(ecef_point.x, ecef_point.y, ecef_point.z);