code: def project(lon0, lat0, lon1, lat1):
      
  
      costh = np.cos(np.pi / 2. - lat0)
      cosphi = np.cos(lon0)
  
      sinth = np.sin(np.pi / 2. - lat0)
      sinphi = np.sin(lon0)
  
      xyz = lonlat_to_xyz(lon1, lat1)
      x1 = xyz[0]
      y1 = xyz[1]
      z1 = xyz[2]
  
      x1p = x1 * costh * cosphi + y1 * costh * sinphi - z1 * sinth
      y1p = -x1 * sinphi + y1 * cosphi
      z1p = x1 * sinth * cosphi + y1 * sinth * sinphi + z1 * costh
  
      r = np.arctan2(np.sqrt(x1p ** 2 + y1p ** 2), z1p)
      phi = np.arctan2(y1p, x1p)
  
      return r * np.cos(phi), r * np.sin(phi)
docs: This function performs a stereographic projection on the unit
      vector (lon1,lat1) with the pole defined at the reference unit
      vector (lon0,lat0).