import math


def buckling_factor_m_n(
        D11, D12, D22, D66, mode_m=1, mode_n=1,
        N_x=0, N_y=0, length_x=1, length_y=1):
    """
    returns the buckling factor of a simply-supported orthotropic laminate
    for a specific buckling mode, based on out-of-plane stiffnesses

    INPUT

    - D11, D12, D22 and D66: out-of-plane stifnesses
    - length_x: plate dimensions (x-direction)
    - length_y: plate dimensions (y-direction)
    - N_x: compressive loading intensity in the x-direction
    - N_y: compressive loading intensity in the x-direction
    - mode_m: buckling mode (number of half-waves) in the x direction
    - mode_n: buckling mode (number of half-waves) in the y direction
    """
    alpha = (mode_m/length_x)**2
    beta = (mode_n/length_y)**2
    return (math.pi**2)*(D11*alpha**2 \
        + 2*(D12 + 2*D66)*alpha*beta \
        + D22*beta**2)/(alpha*abs(N_x) + beta*abs(N_y))