docs: def _gccalc(lon, lat, azimuth, maxdist=None):
      """
      Original javascript on http://williams.best.vwh.net/gccalc.htm
      Translated into python by Thomas Lecocq
      This function is a black box, because trigonometry is difficult
      
      """
      glat1 = lat * np.pi / 180.
      glon1 = lon * np.pi / 180.
      s = maxdist / 1.852243
      faz = azimuth * np.pi / 180.
   
      EPS = 0.00000000005
      if ((np.abs(np.cos(glat1)) < EPS) and not (np.abs(np.sin(faz)) < EPS)):
          raise CourseException("Only North-South courses are meaningful")
  
      a = 6378.137 / 1.852243
      f = 1 / 298.257223563
      r = 1 - f
      tu = r * np.tan(glat1)
      sf = np.sin(faz)
      cf = np.cos(faz)
      if (cf == 0):
          b = 0.
      else:
          b = 2. * np.arctan2 (tu, cf)
  
      cu = 1. / np.sqrt(1 + tu * tu)
      su = tu * cu
      sa = cu * sf
      c2a = 1 - sa * sa
      x = 1. + np.sqrt(1. + c2a * (1. / (r * r) - 1.))
      x = (x - 2.) / x
      c = 1. - x
      c = (x * x / 4. + 1.) / c
      d = (0.375 * x * x - 1.) * x
      tu = s / (r * a * c)
      y = tu
      c = y + 1
      while (np.abs (y - c) > EPS):
          sy = np.sin(y)
          cy = np.cos(y)
          cz = np.cos(b + y)
          e = 2. * cz * cz - 1.
          c = y
          x = e * cy
          y = e + e - 1.
          y = (((sy * sy * 4. - 3.) * y * cz * d / 6. + x) *
              d / 4. - cz) * sy * d + tu
  
      b = cu * cy * cf - su * sy
      c = r * np.sqrt(sa * sa + b * b)
      d = su * cy + cu * sy * cf
      glat2 = (np.arctan2(d, c) + np.pi) % (2*np.pi) - np.pi
      c = cu * cy - su * sy * cf
      x = np.arctan2(sy * sf, c)
      c = ((-3. * c2a + 4.) * f + 4.) * c2a * f / 16.
      d = ((e * cy * c + cz) * sy * c + y) * sa
      glon2 = ((glon1 + x - (1. - c) * d * f + np.pi) % (2*np.pi)) - np.pi
  
      baz = (np.arctan2(sa, b) + np.pi) % (2 * np.pi)
  
      glon2 *= 180./np.pi
      glat2 *= 180./np.pi
      baz *= 180./np.pi
      return (glon2, glat2, baz)