帮我写个WGS84坐标转UTM坐标的c语言代码并加上注释

以下是WGS84坐标转UTM坐标的C语言代码及注释：

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

#define PI 3.14159265358979323846
#define a 6378137.0 // 长半轴
#define f 0.0033528106647474805 // 扁率
#define k0 0.9996 // 比例因子

double deg2rad(double deg) {
    return deg * PI / 180.0;
}

double rad2deg(double rad) {
    return rad * 180.0 / PI;
}

void WGS84toUTM(double lat, double lon, double* east, double* north, int* zone) {
    double b = a * (1 - f); // 短半轴
    double e = sqrt(1 - b * b / a / a); // 第一偏心率
    double e2 = e * e; // 第二偏心率

    // 计算所在的UTM带号
    *zone = (int)floor((lon + 180.0) / 6.0) + 1;
    if (lat >= 56.0 && lat < 64.0 && lon >= 3.0 && lon < 12.0) {
        *zone = 32;
    }
    if (lat >= 72.0 && lat < 84.0) {
        if (lon >= 0.0 && lon < 9.0) {
            *zone = 31;
        } else if (lon >= 9.0 && lon < 21.0) {
            *zone = 33;
        } else if (lon >= 21.0 && lon < 33.0) {
            *zone = 35;
        } else if (lon >= 33.0 && lon < 42.0) {
            *zone = 37;
        }
    }

    double lon0 = (*zone - 1) * 6.0 - 180.0 + 3.0; // 该带号中央经线
    double lonr = deg2rad(lon - lon0); // 经度差弧度值

    // 计算归化纬度
    double tanlat = tan(deg2rad(lat));
    double coslat = cos(deg2rad(lat));
    double sinlat = sin(deg2rad(lat));
    double N = a / sqrt(1 - e2 * sinlat * sinlat); // 卯酉圈曲率半径
    double T = tanlat * tanlat;
    double C = e2 * coslat * coslat / (1 - e2);
    double A = cos(deg2rad(lonr));
    double M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * lat
                   - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * sin(2 * deg2rad(lat))
                   + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * sin(4 * deg2rad(lat))
                   - (35 * e2 * e2 * e2 / 3072) * sin(6 * deg2rad(lat)));

    // 计算UTM坐标
    *east = k0 * N * (A * M + (1 - T + C) * M * M * M / 6.0
                      + (5 - 18 * T + T * T + 72 * C - 58 * e2) * M * M * M * M * M / 120.0) + 500000.0;
    *north = k0 * (N * tanlat * (lonr * lonr / 2.0 + (5 - T + 9 * C + 4 * C * C) * lonr * lonr * lonr * lonr / 24.0
                                + (61 - 58 * T + T * T + 600 * C - 330 * e2) * lonr * lonr * lonr * lonr * lonr * lonr / 720.0));

    if (lat < 0) {
        *north += 10000000.0; // 南半球加上一百万米
    }
}

int main() {
    double lat = 31.2304; // 纬度
    double lon = 121.4737; // 经度
    double east, north; // UTM坐标
    int zone; // UTM带号
    WGS84toUTM(lat, lon, &east, &north, &zone);
    printf("Easting: %.3f\nNorthing: %.3f\nZone: %d\n", east, north, zone);
    return 0;
}

该代码实现了将WGS84坐标系下的经纬度转换为UTM坐标系下的东北坐标，并且输出所在UTM带号。其中涉及到了UTM坐标系的定义和计算公式，代码中有详细的注释阐明每一步的计算过程。