API Reference - Solve Module

def force_equilibrium(slice_df, theta_list, fs_guess=1.5, tol=1e-6, max_iter=50, debug=False, right_facing=False):
    """
    Limit‐equilibrium by force equilibrium in X & Y with variable interslice angles.

    Parameters:
        slice_df (pd.DataFrame): must contain columns
            'alpha' (slice base inclination, degrees),
            'phi'   (slice friction angle, degrees),
            'c'     (cohesion),
            'dl'    (slice base length),
            'w'     (slice weight),
            'u'     (pore force per unit length),
            'dload' (distributed load),
            'beta'  (distributed load inclination, degrees),
            'kw'    (seismic force),
            't'     (tension crack water force),
            'p'     (reinforcement force),
            'h_pile' (pile force, optional), 'theta_p' (pile inclination, RADIANS, optional)
        theta_list (array-like): slice‐boundary force inclinations (degrees),
                                 length must be n+1 if there are n slices
        fs_guess (float): initial guess for factor of safety
        tol (float): convergence tolerance on residual
        max_iter (int): maximum number of Newton (secant) iterations
        debug (bool): print residuals during iteration

    Returns:
        (bool, dict or str):
           - If converged: (True, {'method':'force_equilibrium','FS':<value>})
           - If failed:   (False, "error message")
    """
    import numpy as np

    n = len(slice_df)
    if len(theta_list) != n+1:
        return False, f"theta_list length ({len(theta_list)}) must be n+1 ({n+1})"

    # extract and convert to radians
    alpha   = np.radians(slice_df['alpha'].values)
    phi     = np.radians(slice_df['phi'].values)
    c       = _c_eff(slice_df)
    w       = slice_df['w'].values
    u       = slice_df['u'].values
    dl      = slice_df['dl'].values
    D       = slice_df['dload'].values
    beta    = np.radians(slice_df['beta'].values)
    kw      = slice_df['kw'].values
    T       = slice_df['t'].values
    P       = slice_df['p'].values
    theta   = np.radians(np.asarray(theta_list))

    # Pile forces (backward compatible — zeros if no pile data)
    H_pile  = slice_df['h_pile'].values  if 'h_pile'  in slice_df.columns else np.zeros(n)
    theta_p = slice_df['theta_p'].values if 'theta_p' in slice_df.columns else np.zeros(n)

    # v12 support/load terms (all zero when the columns are absent)
    def _c12(name):
        return slice_df[name].values.astype(float) if name in slice_df.columns else np.zeros(n)
    P_pt  = _c12('p_pt')
    pa_cx = _c12('pa_cx'); pa_cy = _c12('pa_cy')
    pp_cx = _c12('pp_cx'); pp_cy = _c12('pp_cy')
    H_pas = _c12('h_pile_pas')
    LL    = _c12('lload')
    ll_b  = np.radians(_c12('ll_beta'))
    # The march mixes the orientation-normalized frame (alpha, theta) with the
    # stored axial components (normalized to cos(psi) > 0). The horizontal
    # component enters as a resisting magnitude (facing-agnostic), but the real
    # tension on a right-facing slope is the negation of the stored direction,
    # so its VERTICAL component flips (same correction as janbu()).
    if right_facing:
        pa_cy = -pa_cy
        pp_cy = -pp_cy

    N = np.zeros(n)  # normal forces on slice bases (filled by each march)
    Z = np.zeros(n+1)  # interslice forces, Z[0] = 0 by definition (no force entering leftmost slice)

    def residual(FS):
        """Right-side interslice force Z[n] for a given FS (force-closure residual).

        Delegates the per-slice march to the shared `_equilibrium_march`, then
        copies its output into the persistent N/Z arrays (via slice assignment) so
        the admissibility check and the slice_df writes after the root-find see the
        values at the converged FS. Pure extraction — arithmetic is unchanged.
        """
        N_march, Z_march = _equilibrium_march(
            alpha, phi, c, w, u, dl, D, beta, kw, T, P, H_pile, theta_p, theta, FS,
            P_pt=P_pt, pa_cx=pa_cx, pa_cy=pa_cy, pp_cx=pp_cx, pp_cy=pp_cy,
            H_pas=H_pas, LL=LL, ll_b=ll_b)
        N[:] = N_march
        Z[:] = Z_march
        return Z[n]

    if debug:
        r0 = residual(fs_guess)
        print(f"FS_guess={fs_guess:.6f} → residual={r0:.4g}")

    # Root-find FS such that the right-end interslice force Z[n] = 0, by Newton-
    # secant from fs_guess. (A bracketed brentq solver was tried to also catch the
    # genuinely-low-FS surfaces where the single-guess secant diverges, but it
    # resurrects non-physical over-strength roots near FS->0 on surfaces where the
    # secant correctly fails, which the admissibility check below does not reliably
    # separate from legitimate low-FS solutions. Deferred.)
    try:
        FS_opt = newton(residual, fs_guess, tol=tol, maxiter=max_iter)
    except Exception as e:
        return False, f"force_equilibrium failed to converge: {e}"

    # A non-positive factor of safety is unphysical: the secant has converged onto a
    # spurious root (it happens on some steep right-facing surfaces). Reject it.
    if not np.isfinite(FS_opt) or FS_opt <= 0:
        return False, f"force_equilibrium: non-physical factor of safety ({FS_opt})"

    # Re-evaluate at the root so N and Z reflect FS_opt (newton's last call may differ).
    residual(FS_opt)

    # Admissibility guard. A free search can drive the force-equilibrium solver onto
    # grossly non-physical surfaces (a large fraction of slices in base tension, or
    # pervasive interslice tension) and report a spurious low FS. Reject those by
    # EXTENT only — a few negative base normals or a non-monotonic thrust line occur
    # in valid solutions and are NOT rejected (valid benchmark criticals run ~0-4%
    # negative normals, <=20% interslice tension).
    frac_N_neg = float(np.mean(N < 0)) if n else 0.0
    # Interslice "tension" sign is convention-dependent: on right-facing slopes the
    # caller negates theta_list, which flips the sign of Z, so there tension is Z>0.
    Z_int = Z[1:n]
    z_tension = (Z_int > 0) if right_facing else (Z_int < 0)
    frac_Z_tension = float(np.mean(z_tension)) if n > 1 else 0.0
    if frac_N_neg > 0.5 or frac_Z_tension > 0.5:
        return False, (
            "force_equilibrium: inadmissible solution "
            f"({100*frac_N_neg:.0f}% of base normals in tension, "
            f"{100*frac_Z_tension:.0f}% interslice tension)")

    slice_df['n_eff'] = N  # store effective normal forces in slice_df
    slice_df['z'] = Z[:-1]  # store interslice forces in slice_df, adjust length to n slices

    # Report-only admissibility screen (corps/lowe inherit it). The parallel-force
    # march exposes no line of thrust, so only the base-tension and interslice-
    # tension signatures apply; the third is skipped. On right-facing slopes the
    # caller negates theta_list, flipping Z's sign, so tension there is Z>0 — pass
    # the physical-convention Z (tension < 0) the helper expects.
    Z_phys = -Z if right_facing else Z
    warns = _admissibility_warnings(c, N, Z_phys)

    if debug:
        r_opt = residual(FS_opt)
        print(f" Converged FS = {FS_opt:.6f}, residual = {r_opt:.4g}")

    return True, {'FS': FS_opt, 'warnings': warns}

def solve_all(slice_df, rapid=False):
    """
    Executes all available limit equilibrium solution methods sequentially.

    Runs six different limit equilibrium methods on the provided slice dataframe
    and displays the factor of safety for each method. This is useful for comparing
    results across multiple solution approaches.

    Parameters
    ----------
    slice_df : pandas.DataFrame
        Slice dataframe containing all required columns for all methods.
        Must include: 'alpha', 'phi', 'c', 'w', 'u', 'dl', 'dload', 'd_x', 'd_y',
        'beta', 'kw', 't', 'y_t', 'p', 'x_c', 'y_cg', and additional columns
        required for specific methods (e.g., 'r', 'xo', 'yo' for circular methods).

    Returns
    -------
    None
        Results are printed to console for each method.

    Notes
    -----
    Methods executed in order:
    1. Ordinary Method of Slices (OMS)
    2. Bishop's Simplified Method
    3. Janbu's Simplified Method
    4. Corps of Engineers Method
    5. Lowe & Karafiath Method
    6. Spencer's Method
    7. Morgenstern-Price Method

    If any method fails, an error message is displayed but execution continues
    with the remaining methods.
    """
    solve_selected('oms', slice_df, rapid=rapid)
    solve_selected('bishop', slice_df, rapid=rapid)
    solve_selected('janbu', slice_df, rapid=rapid)
    solve_selected('corps', slice_df, rapid=rapid)
    solve_selected('lowe', slice_df, rapid=rapid)
    solve_selected('spencer', slice_df, rapid=rapid)
    solve_selected('mprice', slice_df, rapid=rapid)

def solve_selected(method_name, slice_df, rapid=False):
    """
    Executes a specified limit equilibrium solution method and displays results.

    Parameters
    ----------
    method_name : str
        Name of the solution method function to call. Must be one of:
        'oms', 'bishop', 'janbu', 'spencer', 'corps', 'lowe',
        'mprice'
    slice_df : pandas.DataFrame
        Slice dataframe containing all required columns for the specified method
        (see individual method documentation for column requirements)
    rapid : bool, optional
        If True, performs rapid drawdown analysis using the specified method.
        Default is False.

    Returns
    -------
    dict or str
        If successful: dictionary containing method results (includes 'FS' and method-specific parameters)
        If failed: error message string

    Notes
    -----
    This function automatically prints the factor of safety and method-specific
    parameters to the console. For methods with additional parameters:
    - Spencer: displays theta (interslice force angle)
    - Janbu: displays fo (correction factor)
    - Corps of Engineers: displays theta
    """

    func = globals()[method_name]

    if rapid:
        success, result = rapid_drawdown(slice_df, method_name)
    else:
        success, result = func(slice_df)
    if not success:
        print(f'Error: {result}')
        return result

    if func == oms:
        print(f'OMS: FS={result["FS"]:.3f}')
    elif func == bishop:
        print(f'Bishop: FS={result["FS"]:.3f}')
    elif func == spencer:
        print(f'Spencer: FS={result["FS"]:.3f}, theta={result["theta"]:.2f}')
        for w in result.get('warnings', []):
            print(f'  Spencer admissibility warning: {w}')
    elif func == janbu:
        print(f'Janbu Corrected FS={result["FS"]:.3f}, fo={result["fo"]:.2f}')
    elif func == corps:
        print(f'Corps Engineers: FS={result["FS"]:.3f}, theta={result["theta"]:.2f}')
        for w in result.get('warnings', []):
            print(f'  Corps admissibility warning: {w}')
    elif func == lowe:
        print(f'Lowe & Karafiath: FS={result["FS"]:.3f}')
        for w in result.get('warnings', []):
            print(f'  Lowe & Karafiath admissibility warning: {w}')
    elif func == mprice:
        print(f'Morgenstern-Price ({result["f_type"]}): FS={result["FS"]:.3f}, '
              f'lambda={result["lambda"]:.3f}')
        for w in result.get('warnings', []):
            print(f'  Morgenstern-Price admissibility warning: {w}')

    return result