docs: def orbitproject(x,y,inc,phi=0,psi=0):
      """Transform x,y planar coordinates into observer's coordinate frame.
      
      x,y are coordinates in z=0 plane (plane of the orbit)
  
      observer is at (inc, phi) on celestial sphere (angles in radians);
      psi is orientation of final x-y axes about the (inc,phi) vector.
  
      Returns x,y,z values in observer's coordinate frame, where
      x,y are now plane-of-sky coordinates and z is along the line of sight.
  
      Parameters
      ----------
      x,y : float or arrray-like
          Coordinates to transorm
  
      inc : float or array-like
          Polar angle(s) of observer (where inc=0 corresponds to north pole
          of original x-y plane).  This angle is the same as standard "inclination."
  
      phi : float or array-like, optional
          Azimuthal angle of observer around z-axis
  
      psi : float or array-like, optional
          Orientation of final observer coordinate frame (azimuthal around
          (inc,phi) vector.
  
      Returns
      -------
      x,y,z : ``ndarray``
          Coordinates in observers' frames.  x,y in "plane of sky" and z
          along line of sight.
      """
  
      x2 = x*np.cos(phi) + y*np.sin(phi)
      y2 = -x*np.sin(phi) + y*np.cos(phi)
      z2 = y2*np.sin(inc)
      y2 = y2*np.cos(inc)
  
      xf = x2*np.cos(psi) - y2*np.sin(psi)
      yf = x2*np.sin(psi) + y2*np.cos(psi)
  
      return (xf,yf,z2)