用c写utm转换wpg84

可以使用以下代码来实现：

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

#define pi 3.14159265358979323846

double deg2rad(double deg) {
return deg * pi / 180.;
}

double rad2deg(double rad) {
return rad * 180. / pi;
}

void utm2wgs84(double easting, double northing, int zone, int southhemi, double *lat, double *lon) {
double k = .9996;
double a = 6378137.;
double f = 1. / 298.257223563;
double b = a * (1. - f);
double e = sqrt(1. - (b / a) * (b / a));
double e2 = e * e;
double e4 = e2 * e2;
double e6 = e4 * e2;
double N1, T1, C1, R1, D, M;
double mu, phi1, phi1rad, x, y, eta2, nu2, p, t, c, s, A, M1;
double lat, lon, lonrad;
double x, y;
double dx = 500000.;
double dy = 10000000.;
double k = 1.;
double eps = 1e-12;
double lat1, lat1rad, lat2, lat2rad, lat3, lat3rad, lat4, lat4rad;
double lon1, lon1rad, lon2, lon2rad, lon3, lon3rad, lon4, lon4rad;
double dlat, dlon;
int i;

if (southhemi) {
lat = -80.;
} else {
lat = .;
}
lon = (zone - 1) * 6 - 180 + 3;
lonrad = deg2rad(lon);

x = dx + (zone - 1) * 100000. + 500000.;
y = dy;
x = easting - x;
y = northing - y;

if (southhemi) {
y = -y;
}

mu = y / (a * (1. - e2 / 4. - 3. * e4 / 64. - 5. * e6 / 256.));
phi1 = mu + (3. * e2 / 2. - 27. * e4 / 32. * sin(2. * mu) + 269. * e6 / 512. * sin(4. * mu) - 6607. * e4 / 245760. * sin(6. * mu) + 2205954. * e6 / 368640. * sin(8. * mu));
phi1rad = deg2rad(phi1);

N1 = a / sqrt(1. - e2 * sin(phi1rad) * sin(phi1rad));
T1 = tan(phi1rad) * tan(phi1rad);
C1 = e2 * cos(phi1rad) * cos(phi1rad);
R1 = a * (1. - e2) / pow(1. - e2 * sin(phi1rad) * sin(phi1rad), 1.5);
D = x / (N1 * k);
M = phi1rad - (N1 * tan(phi1rad) / R1) * (D * D / 2. - (5. + 3. * T1 + 10. * C1 - 4. * C1 * C1 - 9. * e2) * D * D * D * D / 24. + (61. + 90. * T1 + 298. * C1 + 45. * T1 * T1 - 252. * e2 - 3. * C1 * C1) * D * D * D * D * D * D / 720.);

lat1 = rad2deg(M);
lon1 = lon + rad2deg(D / cos(phi1rad));

lat1rad = deg2rad(lat1);
nu2 = e2 * cos(lat1rad) * cos(lat1rad);
eta2 = e2 * sin(lat1rad) * sin(lat1rad);
p = a * (1. - e2) / pow(1. - e2 * sin(lat1rad) * sin(lat1rad), 1.5);
t = tan(lat1rad) * tan(lat1rad);
c = eta2 * cos(lat1rad) * cos(lat1rad);
s = pow(cos(lat1rad), 4.);
A = (easting - x) / (k * p);
M1 = lat1rad + (1. - t + c) * pow(A, 2.) / 2. + (5. - 18. * t + t * t + 72. * c - 58. * nu2) * pow(A, 4.) / 24. - (61. - 479. * t + 179. * t * t - t * t * t) * pow(A, 6.) / 720.;

lat2 = rad2deg(M1);
lon2 = lon + A - (1. + 2. * t + c) * pow(A, 3.) / 6. + (5. - 2. * c + 28. * t - 3. * c * c + 8. * nu2 + 24. * t * t) * pow(A, 5.) / 120.;

lat2rad = deg2rad(lat2);
dlat = lat2 - lat1;
dlon = lon2 - lon1;

for (i = ; i < 5; i++) {
lat3 = lat2 - (p * sin(dlat * pi / 180.) / nu2) * (dlon * pi / 180.);
lon3 = lon2 - (p * cos(lat2rad) / (nu2 * cos(lat1rad))) * (dlon * pi / 180.);
lat3rad = deg2rad(lat3);
nu2 = e2 * cos(lat3rad) * cos(lat3rad);
eta2 = e2 * sin(lat3rad) * sin(lat3rad);
p = a * (1. - e2) / pow(1. - e2 * sin(lat3rad) * sin(lat3rad), 1.5);
dlat = lat3 - lat2;
dlon = lon3 - lon2;
lat2 = lat3;
lon2 = lon3;
}

*lat = lat3;
*lon = lon3;
}

int main() {
double easting = 500000.;
double northing = 4649776.;
int zone = 50;
int southhemi = ;
double lat, lon;

utm2wgs84(easting, northing, zone, southhemi, &lat, &lon);

printf("WGS84: %f, %f\n", lat, lon);

return ;
}