code: def _surfdens(self,R,z,phi=0.,t=0.):
          
          r= numpy.sqrt(R**2.+z**2.)
          Rma= numpy.sqrt(R**2.-self.a**2.+0j)
          if Rma == 0.:
              return (-12.*self.a**3-5.*self.a*z**2
                        +numpy.sqrt(1.+z**2/self.a**2)\
                           *(12.*self.a**3-self.a*z**2+2/self.a*z**4))\
                            /30./numpy.pi*z**-5.
          else:
              return self.a*((2.*self.a**2.+R**2.)*Rma**-5\
                                 *(numpy.arctan(z/Rma)-numpy.arctan(self.a*z/r/Rma))
                             +z*(5.*self.a**3.*r-4.*self.a**4
                                 +self.a**2*(2.*r**2.+R**2)
                                 -self.a*r*(5.*R**2.+3.*z**2.)+R**2.*r**2.)
                             /(self.a**2.-R**2.)**2.
                             /(r**2-self.a**2.)**2.).real/4./numpy.pi
docs: 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)