用c语言写一个大地坐标与空间直角坐标相互转换的代码

以下是一个C语言编写的大地坐标与空间直角坐标相互转换的代码示例：

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

#define PI 3.141592653589793

/* 地球椭球体参数 */
#define a 6378137.0        /* 长半轴 */
#define b 6356752.314245   /* 短半轴 */
#define f (a - b) / a      /* 扁率 */

/* 弧度转度数 */
#define Rad2Deg (180.0 / PI)

/* 度数转弧度 */
#define Deg2Rad (PI / 180.0)

/* 由经纬度（单位为度）计算大地坐标 */
void GeodeticToRectangular(double BLH[3], double XYZ[3])
{
    double L = BLH[0] * Deg2Rad;   /* 经度 */
    double B = BLH[1] * Deg2Rad;   /* 纬度 */
    double H = BLH[2];             /* 高程 */

    double e2 = f * (2 - f);       /* 第一偏心率的平方 */
    double N = a / sqrt(1 - e2 * sin(B) * sin(B)); /* 卯酉圈曲率半径 */
    double X = (N + H) * cos(B) * cos(L);
    double Y = (N + H) * cos(B) * sin(L);
    double Z = (N * (1 - e2) + H) * sin(B);

    XYZ[0] = X;
    XYZ[1] = Y;
    XYZ[2] = Z;
}

/* 由大地坐标计算经纬度（单位为度） */
void RectangularToGeodetic(double XYZ[3], double BLH[3])
{
    double X = XYZ[0];
    double Y = XYZ[1];
    double Z = XYZ[2];

    double e2 = f * (2 - f);           /* 第一偏心率的平方 */
    double r2 = X * X + Y * Y;         /* 横坐标与纵坐标的平方和 */
    double v = a * a / sqrt(r2 * (a * a / b * b) + Z * Z);  /* 卯酉圈曲率半径 */
    double alpha = atan(Z / sqrt(r2)); /* 点的地心角 */
    double sinalpha = sin(alpha);
    double cossalpha = cos(alpha);
    double B = atan((Z + v * e2 * sinalpha * sinalpha * sinalpha) / (r2 - a * a / b * b * v * cossalpha * cossalpha * cossalpha)); /* 纬度 */
    double L = atan(Y / X);           /* 经度 */
    double H = sqrt(r2 / cos(B) - a * a / sqrt(1 - e2 * sin(B) * sin(B))); /* 高程 */

    BLH[0] = L * Rad2Deg;
    BLH[1] = B * Rad2Deg;
    BLH[2] = H;
}

int main()
{
    /* 经纬度大地坐标转换为空间直角坐标 */
    double BLH[3] = {120.189, 30.183, 50.0};  /* 经度（度）、纬度（度）、高程（米）*/
    double XYZ[3] = {0.0};  /* 横坐标、纵坐标、高程 */

    GeodeticToRectangular(BLH, XYZ);

    printf("BLH = %lf, %lf, %lf\n", BLH[0], BLH[1], BLH[2]);
    printf("XYZ = %lf, %lf, %lf\n", XYZ[0], XYZ[1], XYZ[2]);

    /* 空间直角坐标转换为经纬度大地坐标 */
    double XYZ2[3] = {1000000.0, 300000.0, 500000.0}; /* 横坐标、纵坐标、高程 */
    double BLH2[3] = {0.0};  /* 经度（度）、纬度（度）、高程（米）*/

    RectangularToGeodetic(XYZ2, BLH2);

    printf("XYZ = %lf, %lf, %lf\n", XYZ2[0], XYZ2[1], XYZ2[2]);
    printf("BLH = %lf, %lf, %lf\n", BLH2[0], BLH2[1], BLH2[2]);

    return 0;
}