利用C++做一个大地坐标与空间直角坐标（包括克拉索夫斯基椭球、IUGG1975椭球、CGCS2000椭球）互转换算例的程序

以下是一个基于C++的大地坐标与空间直角坐标互转换程序，包括了对克拉索夫斯基椭球、IUGG1975椭球、CGCS2000椭球的支持：

#include <iostream>
#include <cmath>

const double PI = 3.141592653589793238462643383279502884;

// 定义椭球体参数结构体
struct Ellipsoid {
    double a;   // 长半轴
    double b;   // 短半轴
    double f;   // 扁率
    double e2;  // 第一偏心率平方
    double e2p; // 第二偏心率平方
};

// 克拉索夫斯基椭球体参数
const Ellipsoid KRASOVSKY = { 6378245.0, 6356863.0188, 1 / 298.3, 0.00669342162297, 0.00673852541468 };

// IUGG1975椭球体参数
const Ellipsoid IUGG1975 = { 6378140.0, 6356755.2882, 1 / 298.257, 0.006694543, 0.006739501 };

// CGCS2000椭球体参数
const Ellipsoid CGCS2000 = { 6378137.0, 6356752.3141, 1 / 298.257222101, 0.00669438002290, 0.00673949677548 };

// 弧度转角度
double rad2deg(double rad) {
    return rad * 180 / PI;
}

// 角度转弧度
double deg2rad(double deg) {
    return deg * PI / 180;
}

// 计算子午线弧长
double getMeridianArc(double a, double b, double f, double lat) {
    double e2 = 2 * f - f * f;
    double n = (a - b) / (a + b);
    double n2 = n * n;
    double n3 = n2 * n;
    double n4 = n3 * n;
    double n5 = n4 * n;
    double A = a * (1 - n + 5.0 / 4 * (n2 - n3) + 81.0 / 64 * (n4 - n5));
    double B = 3.0 / 2 * a * (n - n2 + 7.0 / 8 * (n3 - n4) + 55.0 / 64 * n5);
    double C = 15.0 / 16 * a * (n2 - n3 + 3.0 / 4 * (n4 - n5));
    double D = 35.0 / 48 * a * (n3 - n4 + 11.0 / 16 * n5);
    double E = 315.0 / 512 * a * (n4 - n5);
    double S = A * lat - B * sin(2 * lat) + C * sin(4 * lat) - D * sin(6 * lat) + E * sin(8 * lat);
    return S;
}

// 大地坐标转空间直角坐标
void geodetic2cartesian(double lat, double lon, double h, Ellipsoid& ellipsoid, double& x, double& y, double& z) {
    double radLat = deg2rad(lat);
    double radLon = deg2rad(lon);
    double sinLat = sin(radLat);
    double cosLat = cos(radLat);
    double sinLon = sin(radLon);
    double cosLon = cos(radLon);
    double N = ellipsoid.a / sqrt(1 - ellipsoid.e2 * sinLat * sinLat);
    x = (N + h) * cosLat * cosLon;
    y = (N + h) * cosLat * sinLon;
    z = (N * (1 - ellipsoid.e2) + h) * sinLat;
}

// 空间直角坐标转大地坐标
void cartesian2geodetic(double x, double y, double z, Ellipsoid& ellipsoid, double& lat, double& lon, double& h) {
    double a = ellipsoid.a;
    double b = ellipsoid.b;
    double e2 = ellipsoid.e2;
    double ep2 = ellipsoid.e2p;
    double p = sqrt(x * x + y * y);
    double theta = atan(z * a / (p * b));
    double sinTheta = sin(theta);
    double cosTheta = cos(theta);
    lat = atan((z + ep2 * b * sinTheta * sinTheta * sinTheta) / (p - e2 * a * cosTheta * cosTheta * cosTheta));
    lon = atan(y / x);
    double sinLat = sin(lat);
    double N = a / sqrt(1 - e2 * sinLat * sinLat);
    h = p / cos(lat) - N;
    double deltaLat;
    do {
        deltaLat = lat;
        double M = getMeridianArc(a, b, ellipsoid.f, lat);
        lat = atan((z + ep2 * N * sinLat) / p / (1 - e2 * N / (N + h)));
    } while (fabs(lat - deltaLat) > 1e-10);
}

int main() {
    // 测试数据
    double lat = 40.0;
    double lon = 116.0;
    double h = 100.0;

    // 克拉索夫斯基椭球体参数
    Ellipsoid ellipsoid = KRASOVSKY;
    double x, y, z;
    geodetic2cartesian(lat, lon, h, ellipsoid, x, y, z);
    std::cout << "Krassovsky ellipsoid:" << std::endl;
    std::cout << "Geodetic coordinate: (" << lat << ", " << lon << ", " << h << ")" << std::endl;
    std::cout << "Cartesian coordinate: (" << x << ", " << y << ", " << z << ")" << std::endl;
    cartesian2geodetic(x, y, z, ellipsoid, lat, lon, h);
    std::cout << "Geodetic coordinate: (" << rad2deg(lat) << ", " << rad2deg(lon) << ", " << h << ")" << std::endl;

    // IUGG1975椭球体参数
    ellipsoid = IUGG1975;
    geodetic2cartesian(lat, lon, h, ellipsoid, x, y, z);
    std::cout << "IUGG1975 ellipsoid:" << std::endl;
    std::cout << "Geodetic coordinate: (" << lat << ", " << lon << ", " << h << ")" << std::endl;
    std::cout << "Cartesian coordinate: (" << x << ", " << y << ", " << z << ")" << std::endl;
    cartesian2geodetic(x, y, z, ellipsoid, lat, lon, h);
    std::cout << "Geodetic coordinate: (" << rad2deg(lat) << ", " << rad2deg(lon) << ", " << h << ")" << std::endl;

    // CGCS2000椭球体参数
    ellipsoid = CGCS2000;
    geodetic2cartesian(lat, lon, h, ellipsoid, x, y, z);
    std::cout << "CGCS2000 ellipsoid:" << std::endl;
    std::cout << "Geodetic coordinate: (" << lat << ", " << lon << ", " << h << ")" << std::endl;
    std::cout << "Cartesian coordinate: (" << x << ", " << y << ", " << z << ")" << std::endl;
    cartesian2geodetic(x, y, z, ellipsoid, lat, lon, h);
    std::cout << "Geodetic coordinate: (" << rad2deg(lat) << ", " << rad2deg(lon) << ", " << h << ")" << std::endl;

    return 0;
}

在程序中，通过定义`Ellipsoid`结构体来存储椭球体的参数，包括长半轴`a`、短半轴`b`、扁率`f`、第一偏心率平方`e2`和第二偏心率平方`e2p`。然后，定义了三个常量分别表示克拉索夫斯基椭球、IUGG1975椭球和CGCS2000椭球。

接着，定义了`rad2deg`和`deg2rad`两个函数，用于实现角度和弧度的转换。

然后，定义了`getMeridianArc`函数，用于计算子午线弧长。该函数中使用了公式计算子午线弧长，具体实现可以参考国家测绘局的技术规范。

最后，定义了`geodetic2cartesian`和`cartesian2geodetic`两个函数，分别用于实现大地坐标到空间直角坐标的转换和空间直角坐标到大地坐标的转换。其中，大地坐标到空间直角坐标的转换使用了公式：

$$
\begin{aligned}
x &= (N + h) \cos \varphi \cos \lambda \\
y &= (N + h) \cos \varphi \sin \lambda \\
z &= (N(1-e^2) + h) \sin \varphi
\end{aligned}
$$

其中，$N$表示卯酉圈曲率半径，$h$表示高程，$\varphi$表示纬度，$\lambda$表示经度，$e^2$表示第一偏心率平方。

空间直角坐标到大地坐标的转换使用了迭代法，具体实现可以参考国家测绘局的技术规范。

最后，在`main`函数中，我们分别测试了三个不同的椭球体参数下的大地坐标和空间直角坐标的转换。