def orbit_posvel(Ms,eccs,semimajors,mreds,obspos=None):
Es = Efn(Ms,eccs) <EOL>rs = semimajors*(<NUM_LIT:1>-eccs*np.cos(Es))<EOL>nus = <NUM_LIT:2> * np.arctan2(np.sqrt(<NUM_LIT:1>+eccs)*np.sin(Es/<NUM_LIT:2>),np.sqrt(<NUM_LIT:1>-eccs)*np.cos(Es/<NUM_LIT:2>))<EOL>xs = semimajors*(np.cos(Es) - eccs)         <EOL>ys = semimajors*np.sqrt(<NUM_LIT:1>-eccs**<NUM_LIT:2>)*np.sin(Es)  <EOL>Edots = np.sqrt(G*mreds*MSUN/(semimajors*AU)**<NUM_LIT:3>)/(<NUM_LIT:1>-eccs*np.cos(Es))<EOL>xdots = -semimajors*AU*np.sin(Es)*Edots/<NUM_LIT>  <EOL>ydots = semimajors*AU*np.sqrt(<NUM_LIT:1>-eccs**<NUM_LIT:2>)*np.cos(Es)*Edots/<NUM_LIT> <EOL>n = np.size(xs)<EOL>orbpos = SkyCoord(xs,ys,<NUM_LIT:0>*u.AU,representation='<STR_LIT>',unit='<STR_LIT>')<EOL>orbvel = SkyCoord(xdots,ydots,<NUM_LIT:0>*u.km/u.s,representation='<STR_LIT>',unit='<STR_LIT>')<EOL>if obspos is None:<EOL><INDENT>obspos = random_spherepos(n) <EOL><DEDENT>if type(obspos) == type((<NUM_LIT:1>,<NUM_LIT:2>,<NUM_LIT:3>)):<EOL><INDENT>obspos = SkyCoord(obspos[<NUM_LIT:0>],obspos[<NUM_LIT:1>],obspos[<NUM_LIT:2>],<EOL>representation='<STR_LIT>').represent_as('<STR_LIT>')<EOL><DEDENT>if not hasattr(obspos,'<STR_LIT>'): <EOL><INDENT>obspos = obspos.represent_as('<STR_LIT>')<EOL><DEDENT>psi = rand.random(n)*<NUM_LIT:2>*np.pi  <EOL>x,y,z = orbitproject(orbpos.x,orbpos.y,obspos.theta,obspos.phi,psi)<EOL>vx,vy,vz = orbitproject(orbvel.x,orbvel.y,obspos.theta,obspos.phi,psi)<EOL>return (SkyCoord(x,y,z,representation='<STR_LIT>'),<EOL>SkyCoord(vx,vy,vz,representation='<STR_LIT>'))<EOL>
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.