def _evaluate(self,R,z,phi=0.,t=0.,dR=0,dphi=0):
        """
        NAME:
           _evaluate
        PURPOSE:
           evaluate the potential at (R,z)
        INPUT:
           R - Cylindrical Galactocentric radius
           z - vertical height
           phi - azimuth
           t - time
        OUTPUT:
           potential at (R,z)
        HISTORY:
           2010-04-16 - Written - Bovy (NYU)
           2012-12-26 - New method using Gaussian quadrature between zeros - Bovy (IAS)
        DOCTEST:
           >>> doubleExpPot= DoubleExponentialDiskPotential()
           >>> r= doubleExpPot(1.,0) #doctest: +ELLIPSIS
           ...
           >>> assert( r+1.89595350484)**2.< 10.**-6.
        """
        if True:
            if isinstance(R,float):
                floatIn= True
                R= nu.array([R])
                z= nu.array([z])
            else:
                floatIn= False
            out= nu.empty(len(R))
            indx= (R <= 6.)
            out[True^indx]= self._kp(R[True^indx],z[True^indx])
            R4max= nu.copy(R)
            R4max[(R < 1.)]= 1.
            kmax= self._kmaxFac*self._beta
            for jj in range(len(R)):
                if not indx[jj]: continue
                maxj0zeroIndx= nu.argmin((self._j0zeros-kmax*R4max[jj])**2.) #close enough
                ks= nu.array([0.5*(self._glx+1.)*self._dj0zeros[ii+1] + self._j0zeros[ii] for ii in range(maxj0zeroIndx)]).flatten()
                weights= nu.array([self._glw*self._dj0zeros[ii+1] for ii in range(maxj0zeroIndx)]).flatten()
                evalInt= special.jn(0,ks*R[jj])*(self._alpha**2.+ks**2.)**-1.5*(self._beta*nu.exp(-ks*nu.fabs(z[jj]))-ks*nu.exp(-self._beta*nu.fabs(z[jj])))/(self._beta**2.-ks**2.)
                out[jj]= -2.*nu.pi*self._alpha*nu.sum(weights*evalInt)
            if floatIn: return out[0]
            else: return out
NAME:
           _evaluate
        PURPOSE:
           evaluate the potential at (R,z)
        INPUT:
           R - Cylindrical Galactocentric radius
           z - vertical height
           phi - azimuth
           t - time
        OUTPUT:
           potential at (R,z)
        HISTORY:
           2010-04-16 - Written - Bovy (NYU)
           2012-12-26 - New method using Gaussian quadrature between zeros - Bovy (IAS)
        DOCTEST:
           >>> doubleExpPot= DoubleExponentialDiskPotential()
           >>> r= doubleExpPot(1.,0) #doctest: +ELLIPSIS
           ...
           >>> assert( r+1.89595350484)**2.< 10.**-6.