def _zforce(self,R,z,phi=0,t=0):
        """
        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)
        """
        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