(*    Snyder, John P., (1987),        *)
(*        Map Projections - A Working Manual,        *)
(*        U.S. Geological Survey Professional Paper 1395,        *)
(*        Government Printing Office, Washington.        *)
(*    pp. 61, 63, 265. URL https://pubs.usgs.gov/pp/1395/report.pdf        *)

(*  Calculate UTM & MGRS co-ordinates.  Source-code based on  *)
(*  https://github.com/shahid28/utm-latlng/blob/master/UTMLatLng.js  *)


southMostLat = -80;
zoneWidthUTM = 6 ;
northMostLat = 84 ;
EastWestSize = 180 ;
NorthSouthSize = 90 ;
bandHeightMGRS = 8 ;
zoneWidthUTM = 6 ;


precisionRound(number_, precision_) := Module(
    {}, 
    N(number, precision)
);


getZoneBand(latit_, longitude_) := Module(
    {x, ZoneNumber, MGRSband = ""},
    ZoneNumber = Floor((longitude + EastWestSize) / zoneWidthUTM) + 1;
    If( And(longitude >= 8, longitude <= 13, latit > 54, latit < 58) ,
        ZoneNumber = 32,
    );
    If( And(longitude >= 3, longitude < 12, latit >= 56, latit < 64) ,
        ZoneNumber = 32 ,
        (*  zones above Norway disappear  *)
        If( And(latit >= 72, latit < 84) ,
            If( And(longitude >= 0, longitude < 9) ,
                ZoneNumber = 31 ,
                If( And(longitude >= 9, longitude < 21) ,
                    ZoneNumber = 33 ,
                    If( And(longitude >= 21, longitude < 33) ,
                        ZoneNumber = 35 ,
                        If( And(longitude >= 33, longitude < 42) ,
                            ZoneNumber = 37, );
                    )
                )
            )
        )
    );
    If( latit < southMostLat ,
        If( longitude < 0, MGRSband = "A", MGRSband = "B"),
    );
    If( latit > northMostLat ,
        If( longitude < 0, MGRSband = "Y", MGRSband = "Z"),
    );
    If( Length(MGRSband) == 0,
        x = Floor((latit - southMostLat + bandHeightMGRS) / bandHeightMGRS);
        MGRSband = StringTake("CDEFGHJKLMNPQRSTUVWXX", {x}),
    );
    {ZoneNumber, MGRSband}
);


convertLatLngToUtm(latitude_, longitude_, precision_Integer) := Module(
    {smallA, f, K, falseE, falseN, LatRad, LongRad, eccSquared,
        ZoneNumber, MGRSband, LongOrigin, LongOriginRad, eccPrimeSquared,
        N, T, C, A, M, UTMEasting, UTMNorthing
    },

(*  error checking  *)
    If( Not(Or(NumberQ(latitude), NumberQ(longitude) ) ) ,
        Return("Not a number"), );
    If( latitude < southMostLat ,
        Return("Latitude is too small"), );
    If( latitude > northMostLat ,
        Return("Latitude is too big"), );
    If( longitude < -EastWestSize ,
        Return("East longitude is too small"), );
    If( longitude >= EastWestSize ,
        Return("West longitude is too big"), );

    smallA = 6378137;
    f =  1/298.257223563;  (*  WGS84  *)
    K = 0.9996;
    falseE = 500000;
    falseN = 10000000;
    LatRad = latitude Degree;
    LongRad = longitude Degree;
    eccSquared = f * (2 - f);

    {ZoneNumber, MGRSband} = getZoneBand(latitude, longitude);
    LongOrigin = (ZoneNumber - 1) * zoneWidthUTM - EastWestSize + zoneWidthUTM / 2;
    LongOriginRad = LongOrigin Degree;
    eccPrimeSquared = (eccSquared) / (1 - eccSquared);
    bN = smallA / Sqrt(1 - eccSquared * Sin(LatRad) * Sin(LatRad));
    T = Tan(LatRad) * Tan(LatRad);
    C = eccPrimeSquared * Cos(LatRad) * Cos(LatRad);
    A = Cos(LatRad) * (LongRad - LongOriginRad);

    M = smallA * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 -
          5 * eccSquared * eccSquared * eccSquared / 256) * LatRad
        - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 +
          45 * eccSquared * eccSquared * eccSquared / 1024) * Sin[2 * LatRad]
        + (15 * eccSquared * eccSquared / 256 +
          45 * eccSquared * eccSquared * eccSquared / 1024) * Sin[4 * LatRad]
        - (35 * eccSquared * eccSquared * eccSquared / 3072) * Sin[6 * LatRad]);

    UTMEasting = K * bN * (A + (1 - T + C) * A * A * A / 6
        + (5 - 18 * T + T * T + 72 * C - 58 * eccPrimeSquared) * A * A * A * A * A / 120)
        + falseE;
    UTMNorthing = K * (M + bN * Tan(LatRad) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24
        + (61 - 58 * T + T * T + 600 * C - 330 * eccPrimeSquared) * A * A * A * A * A * A / 720));
    If( latitude < 0,
        UTMNorthing = UTMNorthing + falseN, );
    {ZoneNumber, MGRSband, precisionRound(UTMEasting, precision), precisionRound(UTMNorthing, precision)}
);


convertLatLngToMgrs(latDecDegrees_, longDecDegrees_) := Module(
    {temp, Zone, MGRSband, Easting, Northing, kMSInorthLetters, kMSIeast, kMSInorth, 
        kMSIsize = 100000, kMSIeastRowSize = 8, kMSIrow2start = 15, EastStr, NorthStr,
        kMSIeastLetters = {"S", "T", "U", "V", "W", "X", "Y", "Z", "A", "B",
            "C", "D", "E", "F", "G", "H", "J", "K", "L", "M", "N", "P", "Q", "R"},
        kMSInorthTemp = "FGHJKLMNPQRSTUVABCDE"
    },
    kMSInorthLetters = StringJoin(kMSInorthTemp, StringTake(kMSInorthTemp,
        {1, kMSIrow2start}));
    {Zone, MGRSband, Easting, Northing} = convertLatLngToUtm(latDecDegrees, longDecDegrees, 2);
    kMSIeast = kMSIeastLetters[[Floor(Easting / kMSIsize) +
        Mod(Zone, 3)* kMSIeastRowSize]];
    temp = Mod(Floor(Northing / kMSIsize), 20) + Mod(Zone, 2) * kMSIrow2start + 1;
    kMSInorth = StringTake(kMSInorthLetters, {temp});
    temp = ToString(Mod(Floor(Easting), kMSIsize));
    EastStr = StringTake(StringJoin("00000", temp), {-5, -1});
    temp = ToString(Mod(Floor(Northing), kMSIsize));
    NorthStr = StringTake(StringJoin("00000", temp), {-5, -1});
    {Zone, MGRSband, StringJoin(kMSIeast, kMSInorth), EastStr, NorthStr}
);