code: def angular_separation(lonp1, latp1, lonp2, latp2):
      
  
      xp1, yp1, zp1 = lonlat2xyz(lonp1, latp1)
      xp2, yp2, zp2 = lonlat2xyz(lonp2, latp2)
  
      ## dot products to obtain angles
  
      angles = np.arccos((xp1 * xp2 + yp1 * yp2 + zp1 * zp2))
  
      ## As this is a unit sphere, angle = length
  
      return angles
docs: Compute the angles between lon / lat points p1 and p2 given in radians.
      On the unit sphere, this also corresponds to the great circle distance.
      p1 and p2 can be numpy arrays of the same length.