returns positions in projected AU and velocities in km/s for given mean anomalies

    Returns positions and velocities as SkyCoord objects.  Uses
    ``orbitutils.kepler.Efn`` to calculate eccentric anomalies using interpolation.    

    Parameters
    ----------
    Ms, eccs, semimajors, mreds : float or array-like
        Mean anomalies, eccentricities, semimajor axes (AU), reduced masses (Msun).

    obspos : ``None``, (x,y,z) tuple or ``SkyCoord`` object
        Locations of observers for which to return coordinates.
        If ``None`` then populate randomly on sphere.  If (x,y,z) or
        ``SkyCoord`` object provided, then use those.

    Returns
    -------
    pos,vel : ``SkyCoord``
        Objects representing the positions and velocities, the coordinates
        of which are ``Quantity`` objects that have units.  Positions are in
        projected AU and velocities in km/s.

def orbit_posvel(Ms,eccs,semimajors,mreds,obspos=None):
    """returns positions in projected AU and velocities in km/s for given mean anomalies

    Returns positions and velocities as SkyCoord objects.  Uses
    ``orbitutils.kepler.Efn`` to calculate eccentric anomalies using interpolation.    

    Parameters
    ----------
    Ms, eccs, semimajors, mreds : float or array-like
        Mean anomalies, eccentricities, semimajor axes (AU), reduced masses (Msun).

    obspos : ``None``, (x,y,z) tuple or ``SkyCoord`` object
        Locations of observers for which to return coordinates.
        If ``None`` then populate randomly on sphere.  If (x,y,z) or
        ``SkyCoord`` object provided, then use those.

    Returns
    -------
    pos,vel : ``SkyCoord``
        Objects representing the positions and velocities, the coordinates
        of which are ``Quantity`` objects that have units.  Positions are in
        projected AU and velocities in km/s. 
    """

    Es = Efn(Ms,eccs) #eccentric anomalies by interpolation

    rs = semimajors*(1-eccs*np.cos(Es))
    nus = 2 * np.arctan2(np.sqrt(1+eccs)*np.sin(Es/2),np.sqrt(1-eccs)*np.cos(Es/2))

    xs = semimajors*(np.cos(Es) - eccs)         #AU
    ys = semimajors*np.sqrt(1-eccs**2)*np.sin(Es)  #AU
    
    Edots = np.sqrt(G*mreds*MSUN/(semimajors*AU)**3)/(1-eccs*np.cos(Es))
        
    xdots = -semimajors*AU*np.sin(Es)*Edots/1e5  #km/s
    ydots = semimajors*AU*np.sqrt(1-eccs**2)*np.cos(Es)*Edots/1e5 # km/s
        
    n = np.size(xs)

    orbpos = SkyCoord(xs,ys,0*u.AU,representation='cartesian',unit='AU')
    orbvel = SkyCoord(xdots,ydots,0*u.km/u.s,representation='cartesian',unit='km/s')
    if obspos is None:
        obspos = random_spherepos(n) #observer position
    if type(obspos) == type((1,2,3)):
        obspos = SkyCoord(obspos[0],obspos[1],obspos[2],
                          representation='cartesian').represent_as('physicsspherical')

    if not hasattr(obspos,'theta'): #if obspos not physics spherical, make it 
        obspos = obspos.represent_as('physicsspherical')
        
    #random orientation of the sky 'x-y' coordinates
    psi = rand.random(n)*2*np.pi  

    #transform positions and velocities into observer coordinates
    x,y,z = orbitproject(orbpos.x,orbpos.y,obspos.theta,obspos.phi,psi)
    vx,vy,vz = orbitproject(orbvel.x,orbvel.y,obspos.theta,obspos.phi,psi)

    return (SkyCoord(x,y,z,representation='cartesian'),
            SkyCoord(vx,vy,vz,representation='cartesian'))