code: def haversine(a0,a1):
      
  
      a = [np.asarray(a0,dtype=float),np.asarray(a1,dtype=float)]
  
      for i in range(2):
          if len(a[i].shape) == 1:
              a[i] = np.reshape(a[i],(1,a[i].shape[0]))
  
          if a[i].shape[1] != 2:
              logger.error("Inputs to haversine have the wrong shape!")
              return
  
          a[i] = np.deg2rad(a[i])
  
      distance_km = np.zeros((a[0].shape[0],a[1].shape[0]))
  
      for i in range(a[1].shape[0]):
          b = np.sin((a[0][:,1] - a[1][i,1])/2.)**2 + np.cos(a[0][:,1]) * np.cos(a[1][i,1]) * np.sin((a[0][:,0] - a[1][i,0])/2.)**2
          distance_km[:,i] = 6371.000 * 2 * np.arctan2( np.sqrt(b), np.sqrt(1-b) )
  
      return distance_km
docs: Compute the distance in km between two sets of points in long/lat.
  
      One dimension of a* should be 2; longitude then latitude. Uses haversine
      formula.
  
      Parameters
      ----------
      a0: array/list of at most 2 dimensions
          One dimension must be 2
      a1: array/list of at most 2 dimensions
          One dimension must be 2
  
      Returns
      -------
      distance_km:  np.array
          2-dimensional array of distances in km between points in a0, a1
  
  
      Examples
      --------
  
      haversine([10.1,52.6],[[10.8,52.1],[-34,56.]])
  
      returns array([[   73.15416698,  2836.6707696 ]])