LiveCode function for converting a UTM location to latitude/longitude

utmtolatlong.livecodescript

script "UTM To LatLong Library"
/**
Summary: Converts UTM coordindates to lat/long

Parameters:
pZoneNumber: A number between 1 and 60.
pZoneLetter: A-Z, Zones A-M are soutern, N-Z are northern/
pEasting: Number
pNorthing: Number

Description:

put utmToLatlong(30, "U", 655874, 5469610) into tLatLong
-- 49.359238,-0.853344

Returns: Point
*/
function utmToLatLong pZoneNumber, pZoneLetter, pEasting, pNorthing
  local tNorthernHemisphere
  local tLatitude, tLongitude
  local t_a, t_e, t_e1sq, t_k0
  local tArc, t_mu, t_ei, t_ca, t_cb, t_cc, t_cd, t_phi1
  local t_n0, t_r0, t_fact1, t_a1, t_dd0, t_fact2
  local t_t0, t_Q0, t_fact3, t_fact4, t_lof1, t_lof2, t_lof3, t_a2, t_a3
  local t_zoneCM

  put pZoneLetter is among the chars of "NOPQRSTUVWXYZ" into tNorthernHemisphere
  if not tNorthernHemisphere then
    put 10000000 - pNorthing into pNorthing
  end if

  put 6378137 into t_a
  put 0.081819191 into t_e
  put 0.006739497 into t_e1sq
  put 0.9996 into t_k0

  put pNorthing / t_k0 into tArc
  put tArc / (t_a * (1-t_e^2 / 4.0 - 3 * t_e^4 / 64.0 - 5 * t_e^6 / 256.0)) into t_mu

  put (1- (1 - t_e * t_e)^(1 / 2.0)) / (1 + (1 - t_e * t_e)^(1 / 2.0)) into t_ei

  put 3 * t_ei / 2 - 27 * t_ei^3 / 32.0 into t_ca

  put 21 * t_ei^2 / 16 - 55 * t_ei^4 / 32 into t_cb
  put 151 * t_ei^3 / 96 into t_cc
  put 1097 * t_ei^4 / 512 into t_cd
  put t_mu + t_ca * sin(2 * t_mu) + t_cb * sin(4 * t_mu) + t_cc * sin(6 * t_mu) + t_cd * sin(8 * t_mu) into t_phi1

  put t_a / (1 - (t_e * sin(t_phi1))^2)^(1 / 2.0) into t_n0

  put t_a * (1 - t_e * t_e) / (1 - (t_e * sin(t_phi1))^2)^(3 / 2.0) into t_r0
  put t_n0 * tan(t_phi1) / t_r0 into t_fact1

  put 500000 - pEasting into t_a1
  put t_a1 / (t_n0 * t_k0) into t_dd0
  put t_dd0 * t_dd0 / 2 into t_fact2

  put tan(t_phi1)^2 into t_t0
  put t_e1sq * cos(t_phi1)^2 into t_Q0
  put (5 + 3 * t_t0 + 10 * t_Q0 - 4 * t_Q0 * t_Q0 - 9 * t_e1sq) * t_dd0^4 / 24 into t_fact3

  put (61 + 90 * t_t0 + 298 * t_Q0 + 45 * t_t0 * t_t0 - 252 * t_e1sq - 3 * t_Q0 * t_Q0) * t_dd0^6 / 720 into t_fact4

  put t_a1 / (t_n0 * t_k0) into t_lof1
  put (1 + 2 * t_t0 + t_Q0) * t_dd0^3 / 6.0 into t_lof2
  put (5 - 2 * t_Q0 + 28 * t_t0 - 3 * t_Q0^2 + 8 * t_e1sq + 24 * t_t0^2) * t_dd0^5 / 120 into t_lof3
  put (t_lof1 - t_lof2 + t_lof3) / cos(t_phi1) into t_a2
  put t_a2 * 180 / pi into t_a3

  put 180 * (t_phi1 - t_fact1 * (t_fact2 + t_fact3 + t_fact4)) / pi into tLatitude

  if not tNorthernHemisphere then
    put -tLatitude into tLatitude
  end if

  if pZoneNumber > 0 then
    put 6 * pZoneNumber - 183 into t_zoneCM
  else
    put 3.0 into t_zoneCM
  end if

  put t_zoneCM - t_a3 into tLongitude

  return tLatitude & "," & tLongitude
end utmToLatLong