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.

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)