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

这个算法涉及到了比较复杂的数学公式和计算方法，需要一定的数学基础。以下给出基于克拉索夫斯基椭球的算法实现。

首先，我们需要定义一些常量和变量：

const double PI = 3.14159265358979323846; // 圆周率
const double a = 6378245.0; // 长半轴
const double b = 6356863.0187730473; // 短半轴
const double f = (a - b) / a; // 扁率
const double e2 = 1 - (b * b) / (a * a); // 第一偏心率的平方
const double e12 = (a * a) / (b * b) - 1; // 第二偏心率的平方

接下来是大地坐标转空间直角坐标的实现：

void GeodeticToCartesian(double B, double L, double H, double &X, double &Y, double &Z)
{
    double sinB = sin(B * PI / 180);
    double cosB = cos(B * PI / 180);
    double sinL = sin(L * PI / 180);
    double cosL = cos(L * PI / 180);
    double N = a / sqrt(1 - e2 * sinB * sinB);
    X = (N + H) * cosB * cosL;
    Y = (N + H) * cosB * sinL;
    Z = (N * (1 - e2) + H) * sinB;
}

其中，B、L、H分别为纬度、经度和高程，X、Y、Z为空间直角坐标系下的坐标。

接下来是空间直角坐标转大地坐标的实现：

void CartesianToGeodetic(double X, double Y, double Z, double &B, double &L, double &H)
{
    double p = sqrt(X * X + Y * Y);
    double theta = atan((Z * a) / (p * b));
    double sinTheta = sin(theta);
    double cosTheta = cos(theta);
    double tanB = (Z + e12 * b * sinTheta * sinTheta * sinTheta) / (p - e2 * a * cosTheta * cosTheta * cosTheta);
    B = atan(tanB);
    double sinB = sin(B);
    double cosB = cos(B);
    L = atan2(Y, X);
    double N = a / sqrt(1 - e2 * sinB * sinB);
    H = p / cosB - N;
}

其中，X、Y、Z为空间直角坐标系下的坐标，B、L、H分别为纬度、经度和高程。

以上是基于克拉索夫斯基椭球的算法实现，对于其他椭球，只需要将定义常量的部分进行修改即可。同时，需要注意角度的单位，一般来说，大地坐标的角度单位为度，空间直角坐标的角度单位为弧度。