空间直角坐标与大地坐标转换源程序

以下是一个基于C语言的空间直角坐标与大地坐标转换源程序：

#include <stdio.h>
#include <math.h>

#define pi 3.14159265358979323846

// 定义椭球体参数结构体
typedef struct {
    double a;   // 长半轴
    double b;   // 短半轴
    double f;   // 扁率
} Ellipsoid;

// 定义大地坐标结构体
typedef struct {
    double B;   // 纬度
    double L;   // 经度
    double H;   // 高程
} GeodeticCoord;

// 定义空间直角坐标结构体
typedef struct {
    double X;
    double Y;
    double Z;
} CartesianCoord;

// 计算椭球体的第一偏心率
double Eccentricity(Ellipsoid ellipsoid) {
    double e = sqrt(1 - pow(ellipsoid.b / ellipsoid.a, 2));
    return e;
}

// 计算椭球体的第二偏心率
double SecondEccentricity(Ellipsoid ellipsoid) {
    double e = Eccentricity(ellipsoid);
    double ep = sqrt(pow(ellipsoid.a / ellipsoid.b, 2) - 1);
    return ep;
}

// 计算子午圈半径
double MeridianRadius(double lat, Ellipsoid ellipsoid) {
    double a = ellipsoid.a;
    double b = ellipsoid.b;
    double e2 = pow(Eccentricity(ellipsoid), 2);
    double M = a * (1 - e2);
    double N = a / sqrt(1 - e2 * pow(sin(lat), 2));
    double Rm = M / pow(1 - e2 * pow(sin(lat), 2), 1.5);
    return Rm;
}

// 计算卯酉圈半径
double PrimeVerticalRadius(double lat, Ellipsoid ellipsoid) {
    double a = ellipsoid.a;
    double b = ellipsoid.b;
    double e2 = pow(Eccentricity(ellipsoid), 2);
    double M = a * (1 - e2);
    double N = a / sqrt(1 - e2 * pow(sin(lat), 2));
    double Rn = N / sqrt(1 + e2 * pow(cos(lat), 2));
    return Rn;
}

// 大地坐标转空间直角坐标
CartesianCoord GeodeticToCartesian(GeodeticCoord geodeticCoord, Ellipsoid ellipsoid) {
    double B = geodeticCoord.B * pi / 180;
    double L = geodeticCoord.L * pi / 180;
    double H = geodeticCoord.H;
    double e2 = pow(Eccentricity(ellipsoid), 2);
    double N = ellipsoid.a / sqrt(1 - e2 * pow(sin(B), 2));
    double X = (N + H) * cos(B) * cos(L);
    double Y = (N + H) * cos(B) * sin(L);
    double Z = (N * (1 - e2) + H) * sin(B);
    CartesianCoord cartesianCoord = {X, Y, Z};
    return cartesianCoord;
}

// 空间直角坐标转大地坐标
GeodeticCoord CartesianToGeodetic(CartesianCoord cartesianCoord, Ellipsoid ellipsoid) {
    double X = cartesianCoord.X;
    double Y = cartesianCoord.Y;
    double Z = cartesianCoord.Z;
    double p = sqrt(pow(X, 2) + pow(Y, 2));
    double theta = atan2(Z * ellipsoid.a, p * ellipsoid.b);
    double B = atan2(Z + pow(SecondEccentricity(ellipsoid), 2) * ellipsoid.b * pow(sin(theta), 3), p - pow(Eccentricity(ellipsoid), 2) * ellipsoid.a * pow(cos(theta), 3));
    double L = atan2(Y, X);
    double H = p / cos(B) - ellipsoid.a / sqrt(1 - pow(Eccentricity(ellipsoid), 2) * pow(sin(B), 2));
    GeodeticCoord geodeticCoord = {B * 180 / pi, L * 180 / pi, H};
    return geodeticCoord;
}

int main() {
    // 定义WGS84椭球体参数
    Ellipsoid ellipsoid = {6378137.0, 6356752.3142, 1.0 / 298.257223563};

    // 定义大地坐标
    GeodeticCoord geodeticCoord1 = {39.9075, 116.39723, 0};
    printf("大地坐标: B=%.6f, L=%.6f, H=%.6f\n", geodeticCoord1.B, geodeticCoord1.L, geodeticCoord1.H);

    // 大地坐标转空间直角坐标
    CartesianCoord cartesianCoord1 = GeodeticToCartesian(geodeticCoord1, ellipsoid);
    printf("空间直角坐标: X=%.6f, Y=%.6f, Z=%.6f\n", cartesianCoord1.X, cartesianCoord1.Y, cartesianCoord1.Z);

    // 空间直角坐标转大地坐标
    GeodeticCoord geodeticCoord2 = CartesianToGeodetic(cartesianCoord1, ellipsoid);
    printf("大地坐标: B=%.6f, L=%.6f, H=%.6f\n", geodeticCoord2.B, geodeticCoord2.L, geodeticCoord2.H);

    return 0;
}

其中，WGS84椭球体参数为：

- 长半轴：6378137.0米
- 短半轴：6356752.3142米
- 扁率：1/298.257223563

你可以根据实际需要修改椭球体参数。