/*    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 $
northMostLat: 84 $
EastWestSize: 180 $
NorthSouthSize: 90 $
bandHeightMGRS: 8 $
zoneWidthUTM: 6 $


toRadians(deg) := ev(deg * %pi / 180, numer) $


precisionRound(number, precision) := block([factor],
    factor: 10 ^ precision,
    return(ev(round(number * factor) / factor, numer))
)$


getZoneBand(latit, longitude) := block([ZoneNumber, MGRSband: ""],
    ZoneNumber: floor((longitude + EastWestSize) / zoneWidthUTM) + 1,
    if longitude >= 8 and longitude <= 13 and latitude > 54 and latitude < 58 then
        ZoneNumber: 32,
    if longitude >= 3 and longitude < 12 and latitude >= 56 and latitude < 64 then
        ZoneNumber: 32
    else
        /*  zones above Norway disappear  */
        if latitude >= 72 and latitude < 84 then
            if longitude >= 0 and longitude < 9 then
                ZoneNumber: 31
            else
                if longitude >= 9 and longitude < 21 then
                    ZoneNumber: 33
                else
                    if longitude >= 21 and longitude < 33 then
                        ZoneNumber: 35
                    else
                        if longitude >= 33 and longitude < 42 then
                            ZoneNumber: 37,
    if latit < southMostLat then
        if longitude < 0 then MGRSband: "A" else MGRSband: "B",
    if latit > northMostLat then
        if longitude < 0 then MGRSband: "Y" else MGRSband: "Z",
    if slength(MGRSband) = 0 then block([x],
        x: floor((latit - southMostLat + bandHeightMGRS) / bandHeightMGRS),
        MGRSband: charat("CDEFGHJKLMNPQRSTUVWXX", x)
    ),
    [ZoneNumber, MGRSband]
)$


convertLatLngToUtm(latitude, longitude, precision) := block(
    [small_a, f, K, false_E, false_N, ZoneNumber, LatRad, LongRad,
    eccSquared, A, UTMEasting, UTMNorthing, MGRSband,
    LongOrigin, LongOriginRad, eccPrimeSquared, N, T, C, M
    ],
/*  error checking  */
    if not(numberp(latitude) or numberp(longitude)) then
        return ("Not a number"),
    if latitude < southMostLat then
        return ("Latitude is too small"),
    if latitude > northMostLat then
        return ("Latitude is too big"),
    if longitude < -EastWestSize then
        return ("East longitude is too small"),
    if longitude >= EastWestSize then
        return ("West longitude is too big"),

    small_a : 6378137,
    f :  1/298.257223563,  /*  WGS84  */
    K : 0.9996,
    false_E : 500000,
    false_N : 10000000,

    LatRad: toRadians(latitude),
    LongRad: toRadians(longitude),
    eccSquared: f * (2 - f),

    [ZoneNumber, MGRSband]: getZoneBand(latitude, longitude),
    LongOrigin: (ZoneNumber - 1) * zoneWidthUTM - EastWestSize + zoneWidthUTM / 2,
    LongOriginRad: toRadians(LongOrigin),
    eccPrimeSquared: (eccSquared) / (1 - eccSquared),
    N: ev(small_a / sqrt(1 - eccSquared * sin(LatRad) * sin(LatRad)), numer),
    T: ev(tan(LatRad) * tan(LatRad), numer),
    C: ev(eccPrimeSquared * cos(LatRad) * cos(LatRad), numer),
    A: ev(cos(LatRad) * (LongRad - LongOriginRad), numer),

    M: ev(small_a * ((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)), numer),

    UTMEasting: K * N * (A + (1 - T + C) * A * A * A / 6
        + (5 - 18 * T + T * T + 72 * C - 58 * eccPrimeSquared) * A * A * A * A * A / 120)
        + false_E,
    UTMNorthing: ev(K * (M + N * 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)), numer),
    if latitude < 0 then
        UTMNorthing: UTMNorthing + false_N,
    [ZoneNumber, MGRSband, precisionRound(UTMEasting, precision), precisionRound(UTMNorthing, precision)]
)$


convertLatLngToMgrs(latDecDegrees, longDecDegrees) := block(
    [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: charlist(concat(kMSInorthTemp, substring(kMSInorthTemp, 1, kMSIrow2start + 1))),
    [Zone, MGRSband, Easting, Northing]: convertLatLngToUtm(latDecDegrees, longDecDegrees, 2),
    kMSIeast: kMSIeastLetters[floor(Easting / kMSIsize) +
        (Zone, 3)* kMSIeastRowSize],
    kMSInorth: kMSInorthLetters[mod(floor(Northing / kMSIsize), 20) +
        mod(Zone, 2) * kMSIrow2start + 1],
    temp: sconcat(mod(floor(Easting), kMSIsize)),
    EastStr: concat(substring("0000", slength(temp)), temp),
    temp: sconcat(mod(floor(Northing), kMSIsize)),
    NorthStr: concat(substring("0000", slength(temp)), temp),
    [Zone, MGRSband, sconcat(kMSIeast, kMSInorth), EastStr, NorthStr]
)$