code: def _zforce(self,R,z,phi=0,t=0):
          
          r= numpy.sqrt(R**2.+z**2.)
          out= self._scf.zforce(R,z,phi=phi,use_physical=False)
          for a,s,ds,H,dH in zip(self._Sigma_amp,self._Sigma,self._dSigmadR,
                               self._Hz,self._dHzdz):
              out-= 4.*numpy.pi*a*(ds(r)*H(z)*z/r+s(r)*dH(z))
          return out
docs: NAME:
             _zforce
          PURPOSE:
             evaluate the vertical force at (R,z, phi)
          INPUT:
             R - Cylindrical Galactocentric radius
             z - vertical height
             phi - azimuth
             t - time
          OUTPUT:
             vertical force at (R,z, phi)
          HISTORY:
             2016-12-26 - Written - Bovy (UofT/CCA)