请用C语言计算子午线弧长

在C语言中，计算子午线弧长通常涉及到地理学中的球面三角法，特别是在经纬度坐标系下。对于这个问题，我们可以使用Haversine公式，它是一个近似的算法，用于估算地球上两点之间的最短距离。以下是Haversine公式的简化版C语言函数：

#include <math.h>

#define EARTH_RADIUS 6371 // 地球平均半径，单位公里

double haversin(double lat1, double lon1, double lat2, double lon2) {
    double dLat = (lat2 - lat1) * M_PI / 180;
    double dLon = (lon2 - lon1) * M_PI / 180;

    double a = pow(sin(dLat / 2), 2) + cos(lat1 * M_PI / 180) * cos(lat2 * M_PI / 180) * pow(sin(dLon / 2), 2);
    double c = 2 * atan2(sqrt(a), sqrt(1 - a));

    return EARTH_RADIUS * c; // 返回弧长，单位公里
}

double subparallel_arc_length(double lat1, double lon1, double lat2, double lon2) {
    if (lat1 == lat2) { // 如果经度相同，即在同一纬线上，返回0
        return 0;
    } else {
        return haversin(lat1, lon1, lat2, lon2); // 考虑到纬度变化带来的弧长
    }
}

这个函数`subparallel_arc_length`接收两个起点经纬度`(lat1, lon1)`和终点经纬度`(lat2, lon2)`，并返回两点之间的子午线弧长。