def hs_bounds(K1, K2, M1, M2, f1, f2):
    """
    Hashin Shtrickman bounds for a 2 mineral mixture

    Usage:
        (K_hsu, K_hsl, M_hsu, M_hsl) = hs_bounds(K1, K2, M1, M2, f1, f2)

    Inputs:
        K1 =   bulk modulus mineral 1
        K2 =   bulk modulus mineral 2
        M1 =  shear modulus mineral 1
        M2 =  shear modulus mineral 2
        f1 = volume fraction mineral 1
        f2 = volume fraction mineral 2

    Outputs:
        K_hsu =  Bulk Modulus Hashin-Shtrickman upper bound
        K_hsl =  Bulk Modulus Hashin-Shtrickman lower bound
        M_hsu = Shear Modulus Hashin-Shtrickman upper bound
        M_hsl = Shear Modulus Hashin-Shtrickman lower bound
    """

    K_hsu = K1 + f2 * ( (K2-K1)**-1 + f1*(K1+4*M1/3)**-1 )**-1
    K_hsl = K2 + f1 * ( (K1-K2)**-1 + f2*(K2+4*M2/3)**-1 )**-1
    M_hsu = M1 + f2/( (M2-M1)**-1 + (2*f1*(K1+2*M1))/(5*M1*(K1+4*M1/3)) )
    M_hsl = M2 + f1/( (M1-M2)**-1 + (2*f2*(K2+2*M2))/(5*M2*(K2+4*M2/3)) )

    return K_hsu, K_hsl, M_hsu, M_hsl