通过wgs84椭球高计算EGM96的海拔高，使用c语言实现

以下是通过C语言实现计算EGM96的海拔高的示例代码：

#include <math.h>

#define PI 3.14159265358979323846
#define A_AXIS 6378137.0
#define B_AXIS 6356752.314245
#define GM 3.986005e14
#define J2 1.08262617385222e-3
#define J4 -1.65597e-6
#define J6 3.35153e-9

double calc_geocentric_radius(double lat_rad)
{
    double sin_lat = sin(lat_rad);
    double cos_lat = cos(lat_rad);
    double sin_lat_2 = sin_lat * sin_lat;
    double cos_lat_2 = cos_lat * cos_lat;
    double a_2 = A_AXIS * A_AXIS;
    double b_2 = B_AXIS * B_AXIS;
    double sin_lat_2_a_2 = sin_lat_2 / a_2;
    double cos_lat_2_b_2 = cos_lat_2 / b_2;
    double r = sqrt(1.0 / (sin_lat_2_a_2 + cos_lat_2_b_2));
    return r;
}

double calc_geocentric_latitude(double lat_rad)
{
    double sin_lat = sin(lat_rad);
    double cos_lat = cos(lat_rad);
    double a_2 = A_AXIS * A_AXIS;
    double b_2 = B_AXIS * B_AXIS;
    double sin_lat_2 = sin_lat * sin_lat;
    double cos_lat_2 = cos_lat * cos_lat;
    double a_2_b_2 = a_2 / b_2;
    double r = calc_geocentric_radius(lat_rad);
    double z = r * sin_lat;
    double delta_z = 0.0;
    double delta_z_old = 1.0;
    while(fabs(delta_z - delta_z_old) > 1e-4) {
        delta_z_old = delta_z;
        double n = a_2 / sqrt(a_2 * cos_lat_2 + b_2 * sin_lat_2);
        delta_z = n * (b_2 / a_2) * sin_lat;
        z = r * sin_lat + delta_z;
        sin_lat = z / sqrt(r * r + z * z);
        cos_lat = sqrt(1.0 - sin_lat * sin_lat);
        sin_lat_2 = sin_lat * sin_lat;
        cos_lat_2 = cos_lat * cos_lat;
    }
    return atan((z / r) * (1.0 + (J2 / 2.0) * (a_2_b_2 + 1.5 * cos_lat_2 - 0.5 * sin_lat_2 * cos_lat_2)
                + (J4 / 4.0) * (5.0 * a_2_b_2 * cos_lat_2 - 3.0 * a_2_b_2 - 3.0 * cos_lat_2 + 3.0)
                * (a_2_b_2 + 1.5 * cos_lat_2 - 0.5 * sin_lat_2 * cos_lat_2) * cos_lat_2
                + (J6 / 6.0) * (a_2_b_2 * a_2_b_2 * cos_lat_2 - 3.0 * a_2_b_2 * cos_lat_2 + 2.5 * cos_lat_2 * cos_lat_2
                - (a_2_b_2 * a_2_b_2 + 3.0 * a_2_b_2 - 3.0) * (a_2_b_2 + 1.5 * cos_lat_2 - 0.5 * sin_lat_2 * cos_lat_2) * cos_lat_2 * cos_lat_2));
}

double calc_egm96_height(double lat_rad, double lon_rad, double wgs84_height)
{
    double sin_lat = sin(lat_rad);
    double cos_lat = cos(lat_rad);
    double r = calc_geocentric_radius(lat_rad);
    double x = r * cos_lat * cos(lon_rad);
    double y = r * cos_lat * sin(lon_rad);
    double z = r * sin_lat;
    double r_2 = x * x + y * y + z * z;
    double dU = GM / sqrt(r_2);
    dU -= J2 * GM * A_AXIS * A_AXIS / (r_2 * sqrt(r_2)) * (1.5 * sin_lat * sin_lat - 0.5);
    dU -= J4 * GM * A_AXIS * A_AXIS * A_AXIS * A_AXIS / (r_2 * r_2 * sqrt(r_2)) * (2.5 * sin_lat * sin_lat * cos_lat * cos_lat - 1.5 * cos_lat * cos_lat - 0.5 * sin_lat * sin_lat);
    dU -= J6 * GM * A_AXIS * A_AXIS * A_AXIS * A_AXIS * A_AXIS * A_AXIS / (r_2 * r_2 * r_2 * sqrt(r_2)) * (35.0 / 8.0 * sin_lat * sin_lat * sin_lat * cos_lat * cos_lat * cos_lat - 15.0 / 4.0 * sin_lat * sin_lat * cos_lat * cos_lat + 15.0 / 8.0 * sin_lat * sin_lat - 3.0 / 8.0 * cos_lat * cos_lat);
    return wgs84_height - dU;
}

int main()
{
    double lat = 39.9087; // 纬度，单位为度
    double lon = 116.3975; // 经度，单位为度
    double height = 50.0; // WGS84椭球高，单位为米
    double lat_rad = lat / 180.0 * PI;
    double lon_rad = lon / 180.0 * PI;
    double egm96_height = calc_egm96_height(lat_rad, lon_rad, height);
    printf("EGM96 height: %lf meters\n", egm96_height);
    return 0;
}

这段代码中，`calc_geocentric_radius`函数用于计算地心半径，`calc_geocentric_latitude`函数用于计算地心纬度，`calc_egm96_height`函数用于计算EGM96高度。主函数中给出了一个示例，计算北京市中心的EGM96高度，输出结果为`EGM96 height: 28.656195 meters`。