def surfacemass(self,R,romberg=False,nsigma=None,relative=False):
        """
        NAME:

           surfacemass

        PURPOSE:

           calculate the surface-mass at R by marginalizing over velocity

        INPUT:

           R - radius at which to calculate the surfacemass density (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:

           surface mass at R

        HISTORY:

           2010-03-XX - Written - Bovy (NYU)

        """
        if nsigma == None:
            nsigma= _NSIGMA
        logSigmaR= self.targetSurfacemass(R,log=True,use_physical=False)
        sigmaR2= self.targetSigma2(R,use_physical=False)
        sigmaR1= sc.sqrt(sigmaR2)
        logsigmaR2= sc.log(sigmaR2)
        if relative:
            norm= 1.
        else:
            norm= sc.exp(logSigmaR)
        #Use the asymmetric drift equation to estimate va
        va= sigmaR2/2./R**self._beta*(1./self._gamma**2.-1.
                                      -R*self._surfaceSigmaProfile.surfacemassDerivative(R,log=True)
                                      -R*self._surfaceSigmaProfile.sigma2Derivative(R,log=True))
        if math.fabs(va) > sigmaR1: va = 0.#To avoid craziness near the center
        if romberg:
            return sc.real(bovy_dblquad(_surfaceIntegrand,
                                        self._gamma*(R**self._beta-va)/sigmaR1-nsigma,
                                        self._gamma*(R**self._beta-va)/sigmaR1+nsigma,
                                        lambda x: 0., lambda x: nsigma,
                                        [R,self,logSigmaR,logsigmaR2,sigmaR1,
                                         self._gamma],
                                        tol=10.**-8)/sc.pi*norm)
        else:
            return integrate.dblquad(_surfaceIntegrand,
                                     self._gamma*(R**self._beta-va)/sigmaR1-nsigma,
                                     self._gamma*(R**self._beta-va)/sigmaR1+nsigma,
                                     lambda x: 0., lambda x: nsigma,
                                     (R,self,logSigmaR,logsigmaR2,sigmaR1,
                                      self._gamma),
                                     epsrel=_EPSREL)[0]/sc.pi*norm
NAME:

           surfacemass

        PURPOSE:

           calculate the surface-mass at R by marginalizing over velocity

        INPUT:

           R - radius at which to calculate the surfacemass density (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:

           surface mass at R

        HISTORY:

           2010-03-XX - Written - Bovy (NYU)