code: def _zforce(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._zforce(rr,zz) for rr,zz in zip(R,z)])
                  return out
              if R > 16.*self._hr or R > 6.: return self._kp.zforce(R,z)
              if R < 1.: R4max= 1.
              else: R4max= R
              kmax= self._kmaxFac*self._beta
              maxj0zeroIndx= nu.argmin((self._j0zeros-kmax*R4max)**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= ks*special.jn(0,ks*R)*(self._alpha**2.+ks**2.)**-1.5*(nu.exp(-ks*nu.fabs(z))-nu.exp(-self._beta*nu.fabs(z)))/(self._beta**2.-ks**2.)
              if z > 0.:
                  return -2.*nu.pi*self._alpha*self._beta*nu.sum(weights*evalInt)
              else:
                  return 2.*nu.pi*self._alpha*self._beta*nu.sum(weights*evalInt)
docs: NAME:
             zforce
          PURPOSE:
             evaluate vertical force K_z  (R,z)
          INPUT:
             R - Cylindrical Galactocentric radius
             z - vertical height
             phi - azimuth
             t - time
          OUTPUT:
             K_z (R,z)
          HISTORY:
             2010-04-16 - Written - Bovy (NYU)
          DOCTEST: