docs: 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'))