def capacity_sp_meyerhof_and_hanna_1978(sp, fd, verbose=0):
    """
    Calculates the two-layered foundation capacity according Meyerhof and Hanna (1978)

    :param sp: Soil profile object
    :param fd: Foundation object
    :param wtl: water table level
    :param verbose: verbosity
    :return: ultimate bearing stress
    """
    assert isinstance(sp, sm.SoilProfile)

    sl_0 = sp.layer(1)
    sl_1 = sp.layer(2)
    h0 = sp.layer_depth(2)
    gwl = sp.gwl

    sl_0.nq_factor_0 = (
    (np.tan(np.pi / 4 + np.deg2rad(sl_0.phi / 2))) ** 2 * np.exp(np.pi * np.tan(np.deg2rad(sl_0.phi))))
    if sl_0.phi == 0:
        sl_0.nc_factor_0 = 5.14
    else:
        sl_0.nc_factor_0 = (sl_0.nq_factor_0 - 1) / np.tan(np.deg2rad(sl_0.phi))
    sl_0.ng_factor_0 = (sl_0.nq_factor_0 - 1) * np.tan(1.4 * np.deg2rad(sl_0.phi))

    sl_1.nq_factor_1 = (
        (np.tan(np.pi / 4 + np.deg2rad(sl_1.phi / 2))) ** 2 * np.exp(np.pi * np.tan(np.deg2rad(sl_1.phi))))
    if sl_1.phi == 0:
        sl_1.nc_factor_1 = 5.14
    else:
        sl_1.nc_factor_1 = (sl_1.nq_factor_1 - 1) / np.tan(np.deg2rad(sl_1.phi))
    sl_1.ng_factor_1 = (sl_1.nq_factor_1 - 1) * np.tan(1.4 * np.deg2rad(sl_1.phi))

    if verbose:
        log("Nc: ", sl_1.nc_factor_1)
        log("Nq: ", sl_1.nq_factor_1)
        log("Ng: ", sl_1.ng_factor_1)

    sl_0.kp_0 = (np.tan(np.pi / 4 + np.deg2rad(sl_0.phi / 2))) ** 2
    sl_1.kp_1 = (np.tan(np.pi / 4 + np.deg2rad(sl_1.phi / 2))) ** 2

    # shape factors
    if sl_0.phi >= 10:
        sl_0.s_c_0 = 1 + 0.2 * sl_0.kp_0 * (fd.width / fd.length)
        sl_0.s_q_0 = 1.0 + 0.1 * sl_0.kp_0 * (fd.width / fd.length)
    else:
        sl_0.s_c_0 = 1 + 0.2 * (fd.width / fd.length)
        sl_0.s_q_0 = 1.0
    sl_0.s_g_0 = sl_0.s_q_0

    if sl_1.phi >= 10:
        sl_1.s_c_1 = 1 + 0.2 * sl_1.kp_1 * (fd.width / fd.length)
        sl_1.s_q_1 = 1.0 + 0.1 * sl_1.kp_1 * (fd.width / fd.length)
    else:
        sl_1.s_c_1 = 1 + 0.2 * (fd.width / fd.length)
        sl_1.s_q_1 = 1.0
    sl_1.s_g_1 = sl_1.s_q_1

    # Note: this method explicitly accounts for the foundation depth, so there are no depth factors
    # TODO: inclination factors, see doi.org/10.1139/t78-060

    # Capacity
    a = 1  # assumed to  be one but can range between 1.1 and 1.27 for square footings according to Das (1999) Ch 4
    s = 1

    r = 1 + (fd.width / fd.length)

    # put the same things before that condition
    # effective weight not in the soil object

    if gwl == 0:  # case 1: GWL at surface
        q_at_interface = sl_0.unit_bouy_weight * h0
        unit_eff_weight_0_at_fd_depth = sl_0.unit_bouy_weight
        unit_eff_weight_0_at_interface = sl_0.unit_bouy_weight
        unit_eff_weight_1_below_foundation = sl_1.unit_bouy_weight

    elif 0 < gwl <= fd.depth:  # Case 2: GWL at between foundation depth and surface
        q_at_interface = (sl_0.unit_dry_weight * gwl) + (sl_0.unit_bouy_weight * (h0 - gwl))
        q_d = (sl_0.unit_dry_weight * gwl) + (sl_0.unit_bouy_weight * (fd.depth - gwl))
        unit_eff_weight_0_at_fd_depth = q_d / fd.depth
        unit_eff_weight_0_at_interface = sl_0.unit_bouy_weight
        unit_eff_weight_1_below_foundation = sl_1.unit_bouy_weight

    elif fd.depth < gwl <= fd.width + fd.depth:
        if gwl < h0:  # Case 3: GWL at between foundation depth and foundation depth plus width, and GWL < layer 1 depth

            average_unit_bouy_weight = sl_0.unit_bouy_weight + (
            ((gwl - fd.depth) / fd.width) * (sl_0.unit_dry_weight - sl_0.unit_bouy_weight))

            q_at_interface = (sl_0.unit_dry_weight * gwl) + (sl_0.unit_bouy_weight * (h0 - gwl))
            unit_eff_weight_0_at_fd_depth = sl_0.unit_dry_weight
            unit_eff_weight_0_at_interface = average_unit_bouy_weight
            unit_eff_weight_1_below_foundation = sl_1.unit_bouy_weight

        else:  # Case 4: GWL at between foundation depth and foundation depth plus width, and GWL > layer 1 depth
            average_unit_bouy_weight = sl_1.unit_bouy_weight + (
            ((gwl - h0) / fd.width) * (sl_1.unit_dry_weight - sl_1.unit_bouy_weight))

            q_at_interface = sl_0.unit_dry_weight * h0
            unit_eff_weight_0_at_fd_depth = sl_0.unit_dry_weight
            unit_eff_weight_0_at_interface = sl_0.unit_dry_weight
            unit_eff_weight_1_below_foundation = average_unit_bouy_weight

    elif gwl > fd.depth + fd.width:  # Case 5: GWL beyond foundation depth plus width
        q_at_interface = sl_0.unit_dry_weight * h0
        unit_eff_weight_0_at_fd_depth = sl_0.unit_dry_weight
        unit_eff_weight_0_at_interface = sl_0.unit_dry_weight
        unit_eff_weight_1_below_foundation = sl_1.unit_dry_weight
    else:
        raise ValueError("Could not interpret inputs")  # never reached

    # maximum value (qu <= qt)
    q_ult6 = q_at_interface - unit_eff_weight_0_at_fd_depth * fd.depth
    q_0 = (sl_0.cohesion * sl_0.nc_factor_0) + (0.5 * unit_eff_weight_0_at_interface * fd.width * sl_0.ng_factor_0)
    q_b2 = (q_at_interface * sl_1.nq_factor_1 * sl_1.s_q_1)
    q_1 = (sl_1.cohesion * sl_1.nc_factor_1) + (0.5 * unit_eff_weight_1_below_foundation * fd.width * sl_1.ng_factor_1)
    q_b3 = (unit_eff_weight_1_below_foundation * fd.width * sl_1.ng_factor_1 * sl_1.s_g_1 / 2)
    q_ult5 = r * (unit_eff_weight_0_at_interface * ((h0 - fd.depth) ** 2)) * (1 + (2 * fd.depth / (h0 - fd.depth))) * (
        np.tan(np.deg2rad(sl_0.phi)) / fd.width) * s
    q_t2 = (unit_eff_weight_0_at_fd_depth * fd.depth * sl_0.nq_factor_0 * sl_0.s_q_0)
    q_t3 = (unit_eff_weight_0_at_interface * fd.width * sl_0.ng_factor_0 * sl_0.s_g_0 / 2)

    # qb
    q_b1 = (sl_1.cohesion * sl_1.nc_factor_1 * sl_1.s_c_1)
    q_b = q_b1 + q_b2 + q_b3

    q1_q0 = q_1 / q_0

    # calculate the ca factor
    # if sl_0.cohesion == 0:
    #     c1_c0 = 0
    # else:
    #     c1_c0 = sl_1.cohesion / sl_0.cohesion
    x = np.array([0.000, 0.082, 0.206, 0.298, 0.404, 0.509, 0.598, 0.685, 0.772])
    y = np.array([0.627, 0.700, 0.794, 0.855, 0.912, 0.948, 0.968, 0.983, 0.997])

    # raise Warning("ca should be interpolated using q1/q2 not cohesion, see Figure 4 in MH1978")
    ca_c0 = np.interp(q1_q0, x, y)
    ca = ca_c0 * sl_0.cohesion

    # ks
    x_0 = np.array([0, 20.08, 22.42, 25.08, 27.58, 30.08, 32.58, 34.92, 37.83, 40.00, 42.67, 45.00, 47.00, 49.75])
    y_0 = np.array([0.93, 0.93, 0.93, 0.93, 1.01, 1.17, 1.32, 1.56, 1.87, 2.26, 2.72, 3.35, 3.81, 4.82])
    x_2 = np.array([0, 20.08, 22.50, 25.08, 27.58, 30.08, 32.50, 35.00, 37.67, 40.17, 42.67, 45.00, 47.50, 50.00])
    y_2 = np.array([1.55, 1.55, 1.71, 1.86, 2.10, 2.33, 2.72, 3.11, 3.81, 4.43, 5.28, 6.14, 7.46, 9.24])
    x_4 = np.array([0, 20.00, 22.51, 25.10, 27.69, 30.11, 32.45, 35.04, 37.88, 40.14, 42.65, 45.07, 47.33, 50.08])
    y_4 = np.array([2.49, 2.49, 2.64, 2.87, 3.34, 3.81, 4.43, 5.20, 6.29, 7.38, 9.01, 11.11, 14.29, 19.34])
    x_10 = np.array([0, 20.00, 22.50, 25.08, 28.00, 30.00, 32.50, 34.92, 37.50, 40.17, 42.42, 45.00, 47.17, 50.08])
    y_10 = np.array([3.27, 3.27, 3.74, 4.44, 5.37, 6.07, 7.16, 8.33, 10.04, 12.30, 15.95, 21.17, 27.47, 40.00])
    x_int = sl_0.phi

    if sl_0.phi < 1:
        fd.ks = 0
    else:

        if q1_q0 == 0:
            fd.ks = np.interp(x_int, x_0, y_0)

        elif q1_q0 == 0.2:
            fd.ks = np.interp(x_int, x_2, y_2)

        elif q1_q0 == 0.4:
            fd.ks = np.interp(x_int, x_4, y_4)

        elif q1_q0 == 1.0:
            fd.ks = np.interp(x_int, x_10, y_10)

        elif 0 < q1_q0 < 0.2:
            ks_1 = np.interp(x_int, x_0, y_0)
            ks_2 = np.interp(x_int, x_2, y_2)
            fd.ks = (((ks_2 - ks_1) * q1_q0) / 0.2) + ks_1

        elif 0.2 < q1_q0 < 0.4:
            ks_1 = np.interp(x_int, x_2, y_2)
            ks_2 = np.interp(x_int, x_4, y_4)
            fd.ks = (((ks_2 - ks_1) * (q1_q0 - 0.2)) / 0.2) + ks_1

        elif 0.4 < q1_q0 < 1.0:
            ks_1 = np.interp(x_int, x_4, y_4)
            ks_2 = np.interp(x_int, x_10, y_10)
            fd.ks = (((ks_2 - ks_1) * (q1_q0 - 0.4)) / 0.6) + ks_1
        else:
            raise DesignError(
                "Cannot compute 'ks', bearing ratio out-of-range (q1_q0 = %.3f) required: 0-1." % q1_q0)

    # qu
    q_ult4 = (r * (2 * ca * (h0 - fd.depth) / fd.width) * a)
    q_ult5_ks = q_ult5 * fd.ks
    q_ult = q_b + q_ult4 + q_ult5_ks - q_ult6

    q_t1 = (sl_0.cohesion * sl_0.nc_factor_0 * sl_0.s_c_0)
    q_t = q_t1 + q_t2 + q_t3

    if q_ult > q_t:
        if h0 > fd.width/2:
            fd.q_ult = q_t

        else:
            vert_eff_stress_interface = sp.vertical_effective_stress(h0)
            vert_eff_stress_lowest = sp.vertical_effective_stress(fd.width+fd.depth)
            average_eff_stress = (vert_eff_stress_interface + vert_eff_stress_lowest) / 2

            c_2_eff = sl_1.cohesion + average_eff_stress * np.tan(np.radians(sl_1.phi))

            if sl_0.cohesion > c_2_eff:
                fd.q_ult = q_t

            else:
                # vd = {}
                # vd[1] =[1, 1, 1, 1, 1]
                # vd[0.667] = [1, 1.033, 1.064, 1.088, 1.109]
                # vd[0.5] = [1, 1.056, 1.107, 1.152, 1.193]
                # vd[0.333] = [1, 1.088, 1.167, 1.241, 1.311]
                # vd[0.25] = [1, 1.107, 1.208, 1.302, 1.389]
                # vd[0.2] = [1, 1.121, 1.235, 1.342, 1.444]
                # vd[0.1] = [1, 1.154, 1.302, 1.446, 1.584]

                h_over_b = (h0 - fd.depth)/fd.width
                c1_over_c2 =sl_0.cohesion/c_2_eff

                c_1_over_c_2 = [0.1, 0.2, 0.25, 0.333, 0.5, 0.667, 1.]
                m_1 = [1.584, 1.444, 1.389, 1.311, 1.193, 1.109, 1.]
                m_125 = [1.446, 1.342, 1.302, 1.241, 1.152, 1.088, 1.]
                m_167 = [1.302, 1.235, 1.208, 1.167, 1.107, 1.064, 1.]
                m_25 = [1.154, 1.121, 1.107, 1.088, 1.056, 1.033, 1.]
                m_5 = [1, 1, 1, 1, 1, 1, 1]

                if h_over_b == 0.1:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_1)
                elif h_over_b == 0.125:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_125)
                elif h_over_b == 0.167:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_167)
                elif h_over_b == 0.250:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_25)
                elif h_over_b >= 0.5:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_5)
                elif 0.1 < h_over_b < 0.125:
                    m_a = np.interp(c1_over_c2, c_1_over_c_2, m_1)
                    m_b = np.interp(c1_over_c2, c_1_over_c_2, m_125)
                    m = np.interp(h_over_b, [0.1,0.125], [m_a,m_b])
                elif 0.125 < h_over_b < 0.167:
                    m_a = np.interp(c1_over_c2, c_1_over_c_2, m_125)
                    m_b = np.interp(c1_over_c2, c_1_over_c_2, m_167)
                    m = np.interp(h_over_b, [0.125, 0.167], [m_a, m_b])
                elif 0.167 < h_over_b < 0.25:
                    m_a = np.interp(c1_over_c2, c_1_over_c_2, m_167)
                    m_b = np.interp(c1_over_c2, c_1_over_c_2, m_25)
                    m = np.interp(h_over_b, [0.167, 0.250], [m_a, m_b])
                elif 0.25 < h_over_b < 0.5:
                    m_a = np.interp(c1_over_c2, c_1_over_c_2, m_25)
                    m_b = np.interp(c1_over_c2, c_1_over_c_2, m_5)
                    m = np.interp(h_over_b, [0.250, 0.500], [m_a, m_b])

                fd.q_ult = (sl_0.cohesion * m * sl_0.nc_factor_0) + (unit_eff_weight_0_at_fd_depth * fd.depth)

    else:
        fd.q_ult = q_ult

    return fd.q_ult