code: def anglesep_meeus(lon0: float, lat0: float,
                     lon1: float, lat1: float, deg: bool = True) -> float:
      
  
      if deg:
          lon0 = radians(lon0)
          lat0 = radians(lat0)
          lon1 = radians(lon1)
          lat1 = radians(lat1)
  
      sep_rad = 2 * arcsin(sqrt(haversine(lat0 - lat1) +
                                cos(lat0) * cos(lat1) * haversine(lon0 - lon1)))
  
      if deg:
          return degrees(sep_rad)
      else:
          return sep_rad
docs: Parameters
      ----------
  
      lon0 : float or numpy.ndarray of float
          longitude of first point
      lat0 : float or numpy.ndarray of float
          latitude of first point
      lon1 : float or numpy.ndarray of float
          longitude of second point
      lat1 : float or numpy.ndarray of float
          latitude of second point
      deg : bool, optional
            degrees input/output  (False: radians in/out)
  
      Returns
      -------
  
      sep_rad : float or numpy.ndarray of float
          angular separation
  
  
      Meeus p. 109
  
      from "Astronomical Algorithms" by Jean Meeus Ch. 16 p. 111 (16.5)
      gives angular distance in degrees between two rightAscension,Declination
      points in the sky.  Neglecting atmospheric effects, of course.
  
      Meeus haversine method is stable all the way to exactly 0 deg.
  
      either the arrays must be the same size, or one of them must be a scalar