docs: def geodetic_distance(lons1, lats1, lons2, lats2, diameter=2*EARTH_RADIUS):
      """
      Calculate the geodetic distance between two points or two collections
      of points.
  
      Parameters are coordinates in decimal degrees. They could be scalar
      float numbers or numpy arrays, in which case they should "broadcast
      together".
  
      Implements http://williams.best.vwh.net/avform.htm#Dist
  
      :returns:
          Distance in km, floating point scalar or numpy array of such.
      """
      lons1, lats1, lons2, lats2 = _prepare_coords(lons1, lats1, lons2, lats2)
      distance = numpy.arcsin(numpy.sqrt(
          numpy.sin((lats1 - lats2) / 2.0) ** 2.0
          + numpy.cos(lats1) * numpy.cos(lats2)
          * numpy.sin((lons1 - lons2) / 2.0) ** 2.0
      ))
      return diameter * distance