code: def _zforce(self, R, z, phi=0., t=0.):
          
          if not self.isNonAxi and phi is None:
              phi= 0.
          r, theta, phi = bovy_coords.cyl_to_spher(R,z,phi)
          #x = z
          dr_dz = nu.divide(z,r); dtheta_dz = nu.divide(-R,r**2); dphi_dz = 0
          return self._computeforceArray(dr_dz, dtheta_dz, dphi_dz, R,z,phi)
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-06-06 - Written - Aladdin