code: def surfacemass_z(self,R,nz=7,zmax=1.,fixed_quad=True,fixed_order=8,
                        **kwargs):
          
          if fixed_quad:
              return 2.*integrate.fixed_quad(lambda x: self.density(R*numpy.ones(fixed_order),x,use_physical=False),
                                             0.,.5,n=fixed_order)[0]
          zs= numpy.linspace(0.,zmax,nz)
          sf= numpy.array([self.density(R,z,use_physical=False,
                                        **kwargs) for z in zs])
          lsf= numpy.log(sf)
          #Interpolate
          lsfInterp= interpolate.UnivariateSpline(zs,
                                                  lsf,
                                                  k=3)
          #Integrate
          return 2.*integrate.quad((lambda x: numpy.exp(lsfInterp(x))),
                                   0.,1.)[0]
docs: NAME:
  
             surfacemass_z
  
          PURPOSE:
  
             calculate the vertically-integrated surface density
  
          INPUT:
  
             R - Galactocentric radius (can be Quantity)
  
             fixed_quad= if True (default), use Gauss-Legendre integration
  
             fixed_order= (20), order of GL integration to use
  
             nz= number of zs to use to estimate
  
             zmax= maximum z to use (can be Quantity)
  
             density kwargs
  
          OUTPUT:
  
             \Sigma(R)
  
          HISTORY:
  
             2012-08-30 - Written - Bovy (IAS)