def _surfdens(self,R,z,phi=0.,t=0.):
        """
        NAME:
           _surfdens
        PURPOSE:
           evaluate the surface density for this potential
        INPUT:
           R - Galactocentric cylindrical radius
           z - vertical height
           phi - azimuth
           t - time
        OUTPUT:
           the surface density
        HISTORY:
           2018-08-19 - Written - Bovy (UofT)
        """
        r= numpy.sqrt(R**2.+z**2.)
        Rma= numpy.sqrt(R**2.-self.a**2.+0j)
        if Rma == 0.:
            za2= (z/self.a)**2
            return self.a*(2.+numpy.sqrt(za2+1.)*(za2-2.))/6./numpy.pi/z**3
        else:
            return (self.a*Rma**-3\
                        *(numpy.arctan(self.a*z/r/Rma)-numpy.arctan(z/Rma))
                    +z/(r+self.a)/(R**2.-self.a**2.)).real/2./numpy.pi
NAME:
           _surfdens
        PURPOSE:
           evaluate the surface density for this potential
        INPUT:
           R - Galactocentric cylindrical radius
           z - vertical height
           phi - azimuth
           t - time
        OUTPUT:
           the surface density
        HISTORY:
           2018-08-19 - Written - Bovy (UofT)