code: def pole_from_endpoints(coord1, coord2):
      
      c1 = coord1.cartesian / coord1.cartesian.norm()
  
      coord2 = coord2.transform_to(coord1.frame)
      c2 = coord2.cartesian / coord2.cartesian.norm()
  
      pole = c1.cross(c2)
      pole = pole / pole.norm()
      return coord1.frame.realize_frame(pole)
docs: Compute the pole from a great circle that connects the two specified
      coordinates.
  
      This assumes a right-handed rule from coord1 to coord2: the pole is the
      north pole under that assumption.
  
      Parameters
      ----------
      coord1 : `~astropy.coordinates.SkyCoord`
          Coordinate of one point on a great circle.
      coord2 : `~astropy.coordinates.SkyCoord`
          Coordinate of the other point on a great circle.
  
      Returns
      -------
      pole : `~astropy.coordinates.SkyCoord`
          The coordinates of the pole.