code: def _Rforce(self,R,z,phi=0.,t=0.):
          
          if True:
              if isinstance(R,nu.ndarray):
                  if not isinstance(z,nu.ndarray): z= nu.ones_like(R)*z
                  out= nu.array([self._Rforce(rr,zz) for rr,zz in zip(R,z)])
                  return out
              if (R > 16.*self._hr or R > 6.) and hasattr(self,'_kp'): return self._kp.Rforce(R,z)
              if R < 1.: R4max= 1.
              else: R4max= R
              kmax= self._kmaxFac*self._beta
              kmax= 2.*self._kmaxFac*self._beta
              maxj1zeroIndx= nu.argmin((self._j1zeros-kmax*R4max)**2.) #close enough
              ks= nu.array([0.5*(self._glx+1.)*self._dj1zeros[ii+1] + self._j1zeros[ii] for ii in range(maxj1zeroIndx)]).flatten()
              weights= nu.array([self._glw*self._dj1zeros[ii+1] for ii in range(maxj1zeroIndx)]).flatten()
              evalInt= ks*special.jn(1,ks*R)*(self._alpha**2.+ks**2.)**-1.5*(self._beta*nu.exp(-ks*nu.fabs(z))-ks*nu.exp(-self._beta*nu.fabs(z)))/(self._beta**2.-ks**2.)
              return -2.*nu.pi*self._alpha*nu.sum(weights*evalInt)
docs: NAME:
             Rforce
          PURPOSE:
             evaluate radial force K_R  (R,z)
          INPUT:
             R - Cylindrical Galactocentric radius
             z - vertical height
             phi - azimuth
             t - time
          OUTPUT:
             K_R (R,z)
          HISTORY:
             2010-04-16 - Written - Bovy (NYU)
          DOCTEST: