gistfile1.txt

static const double EARTH_RADIUS = 6371400;

//! WGS84 ellipsoid params (c) JCoord
static const double WGS84_MAJ = 6378137.000;
static const double WGS84_MIN = 6356752.3141;
static const double WGS84_ECC = ((WGS84_MAJ * WGS84_MAJ) - (WGS84_MIN * WGS84_MIN)) / (WGS84_MAJ * WGS84_MAJ);



inline void toUTM(double lat, double lon, double &east, double &north)
{
    double UTM_F0 = 0.9996;
    double a = WGS84_MAJ;
    double eSquared = WGS84_ECC;

    double latRad = lat * (M_PI / 180.0);
    double lonRad = lon * (M_PI / 180.0);

    int longitudeZone = (int) floor((lon + 180.0) / 6.0) + 1;

    // Special zone for Norway
    if (lat >= 56.0 && lat < 64.0 && lon >= 3.0 && lon < 12.0)
        longitudeZone = 32;


    // Special zones for Svalbard
    if (lat >= 72.0 && lat < 84.0)
    {
        if (lon >= 0.0 && lon < 9.0)
            longitudeZone = 31;
        else if (lon >= 9.0 && lon < 21.0)
            longitudeZone = 33;
        else if (lon >= 21.0 && lon < 33.0)
            longitudeZone = 35;
        else if (lon >= 33.0 && lon < 42.0)
            longitudeZone = 37;
    }

    double longitudeOrigin = (longitudeZone - 1) * 6 - 180 + 3;
    double longitudeOriginRad = longitudeOrigin * (M_PI / 180.0);

    double ePrimeSquared = (eSquared) / (1 - eSquared);

    double n = a / sqrt(1 - eSquared * sin(latRad) * sin(latRad));
    double t = tan(latRad) * tan(latRad);
    double c = ePrimeSquared * cos(latRad) * cos(latRad);
    double A = cos(latRad) * (lonRad - longitudeOriginRad);

    double M =
        a
            * ((1 - eSquared / 4 - 3 * eSquared * eSquared / 64 - 5 * eSquared
                * eSquared * eSquared / 256)
                * latRad
                - (3 * eSquared / 8 + 3 * eSquared * eSquared / 32 + 45
                    * eSquared * eSquared * eSquared / 1024)
                * sin(2 * latRad)
                + (15 * eSquared * eSquared / 256 + 45 * eSquared * eSquared
                    * eSquared / 1024) * sin(4 * latRad) - (35
                * eSquared * eSquared * eSquared / 3072)
                * sin(6 * latRad));

    east =
        (UTM_F0
            * n
            * (A + (1 - t + c) * pow(A, 3.0) / 6 + (5 - 18 * t + t * t
                + 72 * c - 58 * ePrimeSquared)
                * pow(A, 5.0) / 120) + 500000.0);

    north =
        (UTM_F0 * (M + n
            * tan(latRad)
            * (A * A / 2 + (5 - t + (9 * c) + (4 * c * c)) * pow(A, 4.0)
                / 24 + (61 - (58 * t) + (t * t) + (600 * c) - (330 * ePrimeSquared))
                * pow(A, 6.0) / 720)));

    // Adjust for the southern hemisphere
    if (lat < 0)
        north += 10000000.0;

}