def capacity_nzs_vm4_2011(sl, fd, h_l=0, h_b=0, vertical_load=1, slope=0, verbose=0, **kwargs):
    """
    calculates the capacity according to
     Appendix B verification method 4 of the NZ building code

    :param sl: Soil object
    :param fd: Foundation object
    :param h_l: Horizontal load parallel to length
    :param h_b: Horizontal load parallel to width
    :param vertical_load: Vertical load
    :param slope: ground slope
    :param verbose: verbosity
    :return: ultimate bearing stress
    """
    # Need to make adjustments if sand  has DR<40% or
    # clay has liquidity indices greater than 0.7

    if not kwargs.get("disable_requires", False):
        models.check_required(sl, ["phi_r", "cohesion", "unit_dry_weight"])
        models.check_required(fd, ["length", "width", "depth"])

    horizontal_load = np.sqrt(h_l ** 2 + h_b ** 2)

    h_eff_b = kwargs.get("h_eff_b", 0)
    h_eff_l = kwargs.get("h_eff_l", 0)
    loc_v_l = kwargs.get("loc_v_l", fd.length / 2)
    loc_v_b = kwargs.get("loc_v_b", fd.width / 2)

    ecc_b = h_b * h_eff_b / vertical_load
    ecc_l = h_l * h_eff_l / vertical_load

    width_eff = min(fd.width, 2 *
                    (loc_v_b + ecc_b), 2 * (fd.width - loc_v_b - ecc_b))
    length_eff = min(fd.length, 2 *
                     (loc_v_l + ecc_l), 2 * (fd.length - loc_v_l - ecc_l))
    area_foundation = length_eff * width_eff

    # check para 3.4.1
    if width_eff / 2 < fd.width / 6:
        raise DesignError("failed on eccentricity")

    # LOAD FACTORS:
    fd.nq_factor = ((np.tan(np.pi / 4 + sl.phi_r / 2)) ** 2 * np.exp(np.pi * np.tan(sl.phi_r)))
    if sl.phi_r == 0:
        fd.nc_factor = 5.14
    else:
        fd.nc_factor = (fd.nq_factor - 1) / np.tan(sl.phi_r)
    fd.ng_factor = 2.0 * (fd.nq_factor - 1) * np.tan(sl.phi_r)

    # shape factors:
    s_c = 1.0 + fd.nq_factor / fd.nc_factor * width_eff / length_eff
    s_q = 1 + width_eff / length_eff * np.tan(sl.phi_r)
    s_g = max(1.0 - 0.4 * width_eff / length_eff, 0.6)  # add limit of 0.6 based on Vesics

    # depth factors:
    if fd.depth / width_eff > 1:
        k = np.arctan(fd.depth / width_eff)
    else:
        k = fd.depth / width_eff
    if sl.phi_r == 0:
        d_c = 1 + 0.4 * k
        d_q = 1.0
    else:
        d_q = (1 + 2 * np.tan(sl.phi_r) *
               (1 - np.sin(sl.phi_r)) ** 2 * k)
        d_c = d_q - (1 - d_q) / (fd.nq_factor * np.tan(sl.phi_r))
    d_g = 1.0

    # load inclination factors:
    if sl.phi_r == 0:
        i_c = 0.5 * (1 + np.sqrt(1 - horizontal_load / (area_foundation * sl.cohesion)))
        i_q = 1.0
        i_g = 1.0
    else:
        if h_b == 0:
            i_q = 1 - horizontal_load / (vertical_load + area_foundation * sl.cohesion /
                                         np.tan(sl.phi_r))
            i_g = i_q
        elif h_b > 0 and h_l == 0:
            i_q = ((1 - 0.7 * horizontal_load / (vertical_load + area_foundation * sl.cohesion /
                                                 np.tan(sl.phi_r))) ** 3)
            i_g = ((1 - horizontal_load / (vertical_load + area_foundation * sl.cohesion /
                                           np.tan(sl.phi_r))) ** 3)
        else:
            raise DesignError("not setup for bi-directional loading")
        i_c = (i_q * fd.nq_factor - 1) / (fd.nq_factor - 1)

    # ground slope factors:
    g_c = 1 - slope * (1.0 - fd.depth / (2 * width_eff)) / 150
    g_q = (1 - np.tan(slope * (1 - fd.depth / (2 * width_eff)))) ** 2
    g_g = g_q

    # stress at footing base:
    q_d = sl.unit_dry_weight * fd.depth

    if verbose:
        log("Nc: ", fd.nc_factor)
        log("Nq: ", fd.nq_factor)
        log("Ng: ", fd.ng_factor)
        log("H: ", horizontal_load)
        log("s_c: ", s_c)
        log("s_q: ", s_q)
        log("s_g: ", s_g)
        log("d_c: ", d_c)
        log("d_q: ", d_q)
        log("d_g: ", d_g)
        log("i_c: ", i_c)
        log("i_q: ", i_q)
        log("i_g: ", i_g)
        log("g_c: ", g_c)
        log("g_q: ", g_q)
        log("g_g: ", g_g)

    # Capacity
    fd.q_ult = (sl.cohesion * fd.nc_factor * s_c * d_c * i_c * g_c +
                q_d * fd.nq_factor * s_q * d_q * i_q * g_q +
                0.5 * width_eff * sl.unit_dry_weight *
                fd.ng_factor * s_g * d_g * i_g * g_g)
    if verbose:
        log("q_ult: ", fd.q_ult)
    return fd.q_ult
calculates the capacity according to
     Appendix B verification method 4 of the NZ building code

    :param sl: Soil object
    :param fd: Foundation object
    :param h_l: Horizontal load parallel to length
    :param h_b: Horizontal load parallel to width
    :param vertical_load: Vertical load
    :param slope: ground slope
    :param verbose: verbosity
    :return: ultimate bearing stress