private static void ConvertUtmToLatitudeLongitude(double northing, double easting, int zoneNumber, out double outLatitude,
out double outLongitude)
{
const double PI = 3.14159265;
const double Rad2Deg = 180.0/PI;
const double k0 = 0.9996;
const double a = 6378137;
const double eccSquared = 0.00669438;
const double eccPrimeSquared = (eccSquared) / (1 – eccSquared);
double e1 = (1 – Math.Sqrt(1 – eccSquared))/(1 + Math.Sqrt(1 – eccSquared));
double x = easting – 500000.0;
double y = northing;
if (zoneNumber < 0)
{
y -= 10000000.0;
}
double longOrigin = (zoneNumber – 1)*6 – 180 + 3;
double M = y/k0;
double mu = M/(a*(1 – eccSquared/4 – 3*eccSquared*eccSquared/64 – 5*eccSquared*eccSquared*eccSquared/256));
double phi1Rad = mu + (3*e1/2 – 27*e1*e1*e1/32)*Math.Sin(2*mu) +
(21*e1*e1/16 – 55*e1*e1*e1*e1/32)*Math.Sin(4*mu) +
(151*e1*e1*e1/96)*Math.Sin(6*mu);
double N1 = a/Math.Sqrt(1 – eccSquared*Math.Sin(phi1Rad)*Math.Sin(phi1Rad));
double T1 = Math.Tan(phi1Rad)*Math.Tan(phi1Rad);
double C1 = eccPrimeSquared*Math.Cos(phi1Rad)*Math.Cos(phi1Rad);
double R1 = a*(1 – eccSquared)/Math.Pow(1 – eccSquared*Math.Sin(phi1Rad)*Math.Sin(phi1Rad), 1.5);
double D = x/(N1*k0);
double Lat = phi1Rad -
(N1*Math.Tan(phi1Rad)/R1)*
(D*D/2 – (5 + 3*T1 + 10*C1 – 4*C1*C1 – 9*eccPrimeSquared)*D*D*D*D/24 +
(61 + 90*T1 + 298*C1 + 45*T1*T1 – 252*eccPrimeSquared – 3*C1*C1)*D*D*D*D*D*D/720);
Lat = Lat*Rad2Deg;
double Long = (D – (1 + 2*T1 + C1)*D*D*D/6 +
(5 – 2*C1 + 28*T1 – 3*C1*C1 + 8*eccPrimeSquared + 24*T1*T1)*D*D*D*D*D/120)/Math.Cos(phi1Rad);
Long = longOrigin + Long*Rad2Deg;
outLatitude = Lat;
outLongitude = Long;
}