code: def skewvR(self,R,romberg=False,nsigma=None,phi=0.):
          
          surfmass= self.surfacemass(R,romberg=romberg,nsigma=nsigma,
                                     use_physical=False)
          vr= self._vmomentsurfacemass(R,1,0,romberg=romberg,nsigma=nsigma)\
              /surfmass
          vr2= self._vmomentsurfacemass(R,2,0,romberg=romberg,nsigma=nsigma)\
              /surfmass
          vr3= self._vmomentsurfacemass(R,3,0,romberg=romberg,nsigma=nsigma)\
              /surfmass
          s2= vr2-vr**2.
          return (vr3-3.*vr*vr2+2.*vr**3.)*s2**(-1.5)
docs: NAME:
  
             skewvR
  
          PURPOSE:
  
             calculate skew in vR at R by marginalizing over velocity
  
          INPUT:
  
             R - radius at which to calculate <vR> (can be Quantity)
  
          OPTIONAL INPUT:
  
             nsigma - number of sigma to integrate the velocities over
  
          KEYWORDS:
  
             romberg - if True, use a romberg integrator (default: False)
  
          OUTPUT:
  
             skewvR
  
          HISTORY:
  
             2011-12-07 - Written - Bovy (NYU)