docs: def simplify(source, kink=20):
      """
       source[] array of geojson points
       kink	in metres, kinks above this depth kept
       kink depth is the height of the triangle abc where a-b and b-c are two consecutive line segments
      """
      source_coord = map(lambda o: {"lng": o.coordinates[0], "lat": o.coordinates[1]}, source)
  
      # count, n_stack, n_dest, start, end, i, sig;
      # dev_sqr, max_dev_sqr, band_sqr;
      # x12, y12, d12, x13, y13, d13, x23, y23, d23;
      F = (math.pi / 180.0) * 0.5
      index = [] # aray of indexes of source points to include in the reduced line
      sig_start = [] # indices of start & end of working section
      sig_end = []
  
      # check for simple cases
      count = len(source_coord)
      if count < 3:
          return source_coord # one or two points
  
      # more complex case. initialize stack
  
      band_sqr = kink * 360.0 / (2.0 * math.pi * 6378137.0) # Now in degrees
      band_sqr *= band_sqr
      n_dest = 0
      sig_start[0] = 0
      sig_end[0] = count - 1
      n_stack = 1
  
      # while the stack is not empty
      while n_stack > 0:
          # ... pop the top-most entries off the stacks
          start = sig_start[n_stack - 1]
          end = sig_end[n_stack - 1]
          n_stack -= 1
  
          if (end - start) > 1: #any intermediate points ?
              # ... yes, so find most deviant intermediate point to either side of line joining start & end points
              x12 = source[end]["lng"] - source[start]["lng"]
              y12 = source[end]["lat"] - source[start]["lat"]
              if math.fabs(x12) > 180.0:
                  x12 = 360.0 - math.fabs(x12)
              x12 *= math.cos(F * (source[end]["lat"] + source[start]["lat"])) # use avg lat to reduce lng
              d12 = (x12 * x12) + (y12 * y12)
  
              i = start + 1
              sig = start
              max_dev_sqr = -1.0
              while i < end:
                  x13 = source[i]["lng"] - source[start]["lng"]
                  y13 = source[i]["lat"] - source[start]["lat"]
                  if math.fabs(x13) > 180.0:
                      x13 = 360.0 - math.fabs(x13)
                  x13 *= math.cos(F * (source[i]["lat"] + source[start]["lat"]))
                  d13 = (x13 * x13) + (y13 * y13)
                  x23 = source[i]["lng"] - source[end]["lng"]
                  y23 = source[i]["lat"] - source[end]["lat"]
                  if math.fabs(x23) > 180.0:
                      x23 = 360.0 - math.fabs(x23)
                  x23 *= math.cos(F * (source[i]["lat"] + source[end]["lat"]))
                  d23 = (x23 * x23) + (y23 * y23)
  
                  if d13 >= (d12 + d23):
                      dev_sqr = d23
                  elif d23 >= (d12 + d13):
                      dev_sqr = d13
                  else:
                      dev_sqr = (x13 * y12 - y13 * x12) * (x13 * y12 - y13 * x12) / d12 # solve triangle
                  if dev_sqr > max_dev_sqr:
                      sig = i
                      max_dev_sqr = dev_sqr
                  i += 1
  
  
              if max_dev_sqr < band_sqr: # is there a sig. intermediate point ?
              #... no, so transfer current start point
                  index[n_dest] = start
                  n_dest += 1
              else: # ... yes, so push two sub-sections on stack for further processing
                  n_stack += 1
                  sig_start[n_stack - 1] = sig
                  sig_end[n_stack - 1] = end
                  n_stack += 1
                  sig_start[n_stack - 1] = start
                  sig_end[n_stack - 1] = sig
  
          else:  # ... no intermediate points, so transfer current start point
              index[n_dest] = start
              n_dest += 1
  
      # transfer last point
      index[n_dest] = count - 1
      n_dest += 1
  
      # make return array
      r = []
      for i in range(0, n_dest):
          r.append(source_coord[index[i]])
  
      return map(lambda o:  {"type": "Point","coordinates": [o.lng, o.lat]}, r)