code: def oortC(self,R,romberg=False,nsigma=None,phi=0.):
          
          #2C= -meanvR/R-dmeanvphi/R/dphi+dmeanvR/dR
          meanvr= self.meanvR(R,romberg=romberg,nsigma=nsigma,phi=phi,
                              use_physical=False)
          dmeanvphiRdphi= 0. #We know this, since the DF does not depend on phi
          surfmass= self._vmomentsurfacemass(R,0,0,phi=phi,romberg=romberg,nsigma=nsigma)
          dmeanvRdR= self._vmomentsurfacemass(R,1,0,deriv='R',phi=phi,romberg=romberg,nsigma=nsigma)/\
              surfmass #other terms is zero because f is even in vR
          return 0.5*(-meanvr/R-dmeanvphiRdphi/R+dmeanvRdR)
docs: NAME:
  
             oortC
  
          PURPOSE:
  
             calculate the Oort function C
  
          INPUT:
  
             R - radius at which to calculate C (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:
  
             Oort C at R
  
          HISTORY:
  
             2011-04-19 - Written - Bovy (NYU)
  
          BUGS:
  
             could be made more efficient, e.g., surfacemass is calculated multiple times
             we know this is zero, but it is calculated anyway (bug or feature?)