docs: def great_circle_distance(lon1, lat1, lon2, lat2):
      """Calculate the great circle distance between one or multiple pairs of
      points given in spherical coordinates. Spherical coordinates are expected
      in degrees. Angle definition follows standard longitude/latitude definition.
      This uses the arctan version of the great-circle distance function
      (en.wikipedia.org/wiki/Great-circle_distance) for increased
      numerical stability.
  
      Parameters
      ----------
      lon1: float scalar or numpy array
          Longitude coordinate(s) of the first element(s) of the point
          pair(s), given in degrees.
      lat1: float scalar or numpy array
          Latitude coordinate(s) of the first element(s) of the point
          pair(s), given in degrees.
      lon2: float scalar or numpy array
          Longitude coordinate(s) of the second element(s) of the point
          pair(s), given in degrees.
      lat2: float scalar or numpy array
          Latitude coordinate(s) of the second element(s) of the point
          pair(s), given in degrees.
  
      Calculation of distances follows numpy elementwise semantics, so if
      an array of length N is passed, all input parameters need to be
      arrays of length N or scalars.
  
  
      Returns
      -------
      distance: float scalar or numpy array
          The great circle distance(s) (in degrees) between the
          given pair(s) of points.
  
      """
      # Convert to radians:
      lat1 = np.array(lat1)*np.pi/180.0
      lat2 = np.array(lat2)*np.pi/180.0
      dlon = (lon1-lon2)*np.pi/180.0
  
      # Evaluate trigonometric functions that need to be evaluated more
      # than once:
      c1 = np.cos(lat1)
      s1 = np.sin(lat1)
      c2 = np.cos(lat2)
      s2 = np.sin(lat2)
      cd = np.cos(dlon)
  
      # This uses the arctan version of the great-circle distance function
      # from en.wikipedia.org/wiki/Great-circle_distance for increased
      # numerical stability.
      # Formula can be obtained from [2] combining eqns. (14)-(16)
      # for spherical geometry (f=0).
  
      return 180.0 / np.pi * np.arctan2(np.sqrt((c2*np.sin(dlon))**2 + (c1*s2-s1*c2*cd)**2), s1*s2+c1*c2*cd)