code: def reference_to_greatcircle(reference_frame, greatcircle_frame):
      
  
      # Define rotation matrices along the position angle vector, and
      # relative to the origin.
      pole = greatcircle_frame.pole.transform_to(coord.ICRS)
      ra0 = greatcircle_frame.ra0
      center = greatcircle_frame.center
      R_rot = rotation_matrix(greatcircle_frame.rotation, 'z')
  
      if not np.isnan(ra0):
          xaxis = np.array([np.cos(ra0), np.sin(ra0), 0.])
          zaxis = pole.cartesian.xyz.value
          if np.abs(zaxis[2]) >= 1e-15:
              xaxis[2] = -(zaxis[0]*xaxis[0] + zaxis[1]*xaxis[1]) / zaxis[2] # what?
          else:
              xaxis[2] = 0.
          xaxis = xaxis / np.sqrt(np.sum(xaxis**2))
          yaxis = np.cross(zaxis, xaxis)
          R = np.stack((xaxis, yaxis, zaxis))
  
      elif center is not None:
          R1 = rotation_matrix(pole.ra, 'z')
          R2 = rotation_matrix(90*u.deg - pole.dec, 'y')
          Rtmp = matrix_product(R2, R1)
  
          rot = center.cartesian.transform(Rtmp)
          rot_lon = rot.represent_as(coord.UnitSphericalRepresentation).lon
          R3 = rotation_matrix(rot_lon, 'z')
          R = matrix_product(R3, R2, R1)
  
      else:
          R1 = rotation_matrix(pole.ra, 'z')
          R2 = rotation_matrix(pole.dec, 'y')
          R = matrix_product(R2, R1)
  
      return matrix_product(R_rot, R)
docs: Convert a reference coordinate to a great circle frame.