-- gridcalc.exu - A series of routines for calculating Maidenhead grid
-- squares, distances, and bearings.

-- Written by Michael J. Sabal
-- adapted on JavaScript by Chris Veness
    -- http://www.movable-type.co.uk/scripts/LatLong.html
-- based on the Haversine algorithms
-- 27 February 2005

include misc.e

-------------------------------------------------------------------------------

global function even(atom a)

  sequence s
  s = sprintf("%d",a)
  if find(s[length(s)],{'0','2','4','6','8'}) > 0 then
    return 1
  else
    return 0    
  end if

end function

-------------------------------------------------------------------------------

global function odd(atom a)
  return not even(a)
end function
  
-------------------------------------------------------------------------------

global function deg2rad(atom deg)
  return deg * PI / 180
end function

-------------------------------------------------------------------------------

global function rad2deg(atom rad)
  return rad * 180 / PI
end function

-------------------------------------------------------------------------------

global function miles2km(atom miles)
  return miles * 1.6093
end function

-------------------------------------------------------------------------------

global function km2miles(atom km)
  return km * 0.6214
end function

-------------------------------------------------------------------------------

global function grid2deg(sequence grid)

  atom lat, long, val
  lat = 0
  long = 0
  
  if odd(length(grid)) or length(grid) < 4 then
    return {0,0}
  end if

  if not integer(grid[1]) or not integer(grid[2]) or not integer(grid[3]) or
     not integer(grid[4]) then
    return {0,0}
  end if

  if grid[1] >= 'a' and grid[1] <= 'z' then
    grid[1] = grid[1] - 32
  end if
  if grid[2] >= 'a' and grid[2] <= 'z' then
    grid[2] = grid[2] - 32
  end if
  
  long = ((grid[1] - 'A') * 20) - 180
  lat = ((grid[2] - 'A') * 10) - 90
  
  val = ((grid[3] - '0') * 2)
  long = long + val
  val = (grid[4] - '0')
  lat = lat + val
  
  if length(grid) = 4 then
    return {lat+.5,long+1} -- return centered in the grid
  end if
  
  if grid[5] >= 'A' and grid[5] <= 'Z' then
    grid[5] = grid[5] + 32
  end if
  if grid[6] >= 'A' and grid[6] <= 'Z' then
    grid[6] = grid[6] + 32
  end if

  val = ((grid[5] - 'a') * (1 / 12)) + (1 / 24)
  long = long + val

  val = ((grid[6] - 'a') * (1/24)) + (1/48)
  lat = lat + val
  
  return {lat,long}

end function

-------------------------------------------------------------------------------

global function deg2grid(atom lat, atom long)

  atom tlat, tlong
  sequence grid  
  grid = ""
  
  tlong = floor(((long + 180)) / 20) + 'A'
  tlat = floor((lat + 90) / 10) + 'A'
   
  grid = tlong & tlat
  
  tlong = floor(remainder((long+180),20)/2) + '0'
  tlat = floor(remainder(lat+90,10)) + '0'
  
  grid = grid & tlong & tlat
  
  tlong = floor((((long+180)) - floor((long+180))) * 60 / 5) + 'a'
  tlat = floor(((lat+90) - floor(lat+90)) * 60 / 2.5) + 'a'
  
  if odd(floor((long+180))) then
    tlong = tlong + 12
  end if
  
  grid = grid & tlong & tlat
  
  return grid
  
end function

-------------------------------------------------------------------------------

function sf(atom val, integer fig)
  -- The Haversine algorithms work best when there are only 4 significant
  -- digits
  
  atom scale, mult
  scale = log(val) * 0.43429448190325181667
  if not integer(scale) then
    scale = floor(scale+1) -- ceiling()
  end if
  mult = power(10,fig-scale)
  return floor(((val*mult)/mult)+.5) -- round()
  
end function  

-------------------------------------------------------------------------------

global function dist(atom lat1, atom long1, atom lat2, atom long2)
    -- each value should be in decimal degrees.

-- http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1

  atom R, dLat, dLong, a, c, d
  R = 6371 -- Earth's radius in km
  dLat = lat2 - lat1
  dLong = long2 - long1

  if dLat = 0 and dLong = 0 then
    return 0
  end if
  
  a = (sin(deg2rad(dLat/2)) * sin (deg2rad(dLat/2))) + 
       (cos(deg2rad(lat1)) * cos(deg2rad(lat2)) * 
       sin(deg2rad(dLong/2)) * sin(deg2rad(dLong/2)))
  c = 2 * arctan(sqrt(a) / sqrt(1-a))
  d = R * c
  
  return d
  
end function

-------------------------------------------------------------------------------

global function bearing(atom lat1, atom long1, atom lat2, atom long2)
  -- Algorithm based on Ed Williams' Aviation Formulary 
  -- http://williams.best.vwh.net/avform.htm#Crs
  
  atom y,x,b
  
  if lat1 = lat2 and long1 = long2 then
    return 0
  end if
  
  y = sin(deg2rad(long1-long2)) * cos(deg2rad(lat2))
  x = (cos(deg2rad(lat1)) * sin(deg2rad(lat2))) - 
       (sin(deg2rad(lat1)) * cos(deg2rad(lat2)) *
       cos(deg2rad(long1 - long2)))
  b = arctan(-y / x)    -- use -y b/c Williams treats West as +ve
  
  if lat2 < lat1 then
    b = b + PI -- remember that tan repeats every 180deg.
  end if
  
  return rad2deg(remainder(b+(2*PI),2*PI))

end function

-------------------------------------------------------------------------------