def _evaluate(self,*args,**kwargs):
        """
        NAME:
           __call__ (_evaluate)
        PURPOSE:
           evaluate the actions (jr,lz,jz)
        INPUT:
           Either:
              a) R,vR,vT,z,vz[,phi]:
                 1) floats: phase-space value for single object (phi is optional) (each can be a Quantity)
                 2) numpy.ndarray: [N] phase-space values for N objects (each can be a Quantity)
              b) Orbit instance: initial condition used if that's it, orbit(t) if there is a time given as well as the second argument 
           delta= (object-wide default) can be used to override the object-wide focal length; can also be an array with length N to allow different delta for different phase-space points
           u0= (None) if object-wide option useu0 is set, u0 to use (if useu0 and useu0 is None, a good value will be computed)
           c= (object-wide default, bool) True/False to override the object-wide setting for whether or not to use the C implementation
           order= (object-wide default, int) number of points to use in the Gauss-Legendre numerical integration of the relevant action integrals  
           When not using C:
              fixed_quad= (False) if True, use Gaussian quadrature (scipy.integrate.fixed_quad instead of scipy.integrate.quad)
              scipy.integrate.fixed_quad or .quad keywords
        OUTPUT:
           (jr,lz,jz)
        HISTORY:
           2012-11-27 - Written - Bovy (IAS)
           2017-12-27 - Allowed individual delta for each point - Bovy (UofT)
        """
        delta= kwargs.pop('delta',self._delta)
        order= kwargs.get('order',self._order)
        if ((self._c and not ('c' in kwargs and not kwargs['c']))\
                or (ext_loaded and (('c' in kwargs and kwargs['c'])))) \
                and _check_c(self._pot):
            if len(args) == 5: #R,vR.vT, z, vz
                R,vR,vT, z, vz= args
            elif len(args) == 6: #R,vR.vT, z, vz, phi
                R,vR,vT, z, vz, phi= args
            else:
                self._parse_eval_args(*args)
                R= self._eval_R
                vR= self._eval_vR
                vT= self._eval_vT
                z= self._eval_z
                vz= self._eval_vz
            if isinstance(R,float):
                R= nu.array([R])
                vR= nu.array([vR])
                vT= nu.array([vT])
                z= nu.array([z])
                vz= nu.array([vz])
            Lz= R*vT
            if self._useu0:
                #First calculate u0
                if 'u0' in kwargs:
                    u0= nu.asarray(kwargs['u0'])
                else:
                    E= nu.array([_evaluatePotentials(self._pot,R[ii],z[ii])
                                 +vR[ii]**2./2.+vz[ii]**2./2.+vT[ii]**2./2. for ii in range(len(R))])
                    u0= actionAngleStaeckel_c.actionAngleStaeckel_calcu0(\
                        E,Lz,self._pot,delta)[0]
                kwargs.pop('u0',None)
            else:
                u0= None
            jr, jz, err= actionAngleStaeckel_c.actionAngleStaeckel_c(\
                self._pot,delta,R,vR,vT,z,vz,u0=u0,order=order)
            if err == 0:
                return (jr,Lz,jz)
            else: #pragma: no cover
                raise RuntimeError("C-code for calculation actions failed; try with c=False")
        else:
            if 'c' in kwargs and kwargs['c'] and not self._c: #pragma: no cover
                warnings.warn("C module not used because potential does not have a C implementation",galpyWarning)
            kwargs.pop('c',None)
            if (len(args) == 5 or len(args) == 6) \
                    and isinstance(args[0],nu.ndarray):
                ojr= nu.zeros((len(args[0])))
                olz= nu.zeros((len(args[0])))
                ojz= nu.zeros((len(args[0])))
                for ii in range(len(args[0])):
                    if len(args) == 5:
                        targs= (args[0][ii],args[1][ii],args[2][ii],
                                args[3][ii],args[4][ii])
                    elif len(args) == 6:
                        targs= (args[0][ii],args[1][ii],args[2][ii],
                                args[3][ii],args[4][ii],args[5][ii])
                    tkwargs= copy.copy(kwargs)
                    try:
                        tkwargs['delta']= delta[ii]
                    except TypeError:
                        tkwargs['delta']= delta
                    tjr,tlz,tjz= self(*targs,**tkwargs)
                    ojr[ii]= tjr
                    ojz[ii]= tjz
                    olz[ii]= tlz
                return (ojr,olz,ojz)
            else:
                #Set up the actionAngleStaeckelSingle object
                aASingle= actionAngleStaeckelSingle(*args,pot=self._pot,
                                                     delta=delta)
                return (aASingle.JR(**copy.copy(kwargs)),
                        aASingle._R*aASingle._vT,
                        aASingle.Jz(**copy.copy(kwargs)))
NAME:
           __call__ (_evaluate)
        PURPOSE:
           evaluate the actions (jr,lz,jz)
        INPUT:
           Either:
              a) R,vR,vT,z,vz[,phi]:
                 1) floats: phase-space value for single object (phi is optional) (each can be a Quantity)
                 2) numpy.ndarray: [N] phase-space values for N objects (each can be a Quantity)
              b) Orbit instance: initial condition used if that's it, orbit(t) if there is a time given as well as the second argument 
           delta= (object-wide default) can be used to override the object-wide focal length; can also be an array with length N to allow different delta for different phase-space points
           u0= (None) if object-wide option useu0 is set, u0 to use (if useu0 and useu0 is None, a good value will be computed)
           c= (object-wide default, bool) True/False to override the object-wide setting for whether or not to use the C implementation
           order= (object-wide default, int) number of points to use in the Gauss-Legendre numerical integration of the relevant action integrals  
           When not using C:
              fixed_quad= (False) if True, use Gaussian quadrature (scipy.integrate.fixed_quad instead of scipy.integrate.quad)
              scipy.integrate.fixed_quad or .quad keywords
        OUTPUT:
           (jr,lz,jz)
        HISTORY:
           2012-11-27 - Written - Bovy (IAS)
           2017-12-27 - Allowed individual delta for each point - Bovy (UofT)