code: def skewvT(self,R,romberg=False,nsigma=None,phi=0.):
          
          surfmass= self.surfacemass(R,romberg=romberg,nsigma=nsigma,
                                     use_physical=False)
          vt= self._vmomentsurfacemass(R,0,1,romberg=romberg,nsigma=nsigma)\
              /surfmass
          vt2= self._vmomentsurfacemass(R,0,2,romberg=romberg,nsigma=nsigma)\
              /surfmass
          vt3= self._vmomentsurfacemass(R,0,3,romberg=romberg,nsigma=nsigma)\
              /surfmass
          s2= vt2-vt**2.
          return (vt3-3.*vt*vt2+2.*vt**3.)*s2**(-1.5)
docs: NAME:
  
             skewvT
  
          PURPOSE:
  
             calculate skew in vT 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:
  
             skewvT
  
          HISTORY:
  
             2011-12-07 - Written - Bovy (NYU)