import sympy


def airystressint(XA, XB, YA, YB, P, a, b, E, nu):
    """
    Calculate integrals of strain energy values based on an Airy stress function
    """
    x1, x2, xa, xb, ya, yb = sympy.symbols("x_1, x_2 x_a x_b y_a y_b")

    sigma = sympy.Matrix([
        [3 * P / (2 * a**3 * b) * x1 * x2],
        [0],
        [3 * P / (4 * a * b) * (1 - (x2**2) / (a**2))]
    ])

    C = 3 * P / (4 * E * a**3 * b)
    epsilon = sympy.Matrix([
        [2 * C * x1 * x2],
        [-2 * nu * C * x1 * x2],
        [2 * (1 + nu) * C * (a**2 - x2**2)]
    ])

    SE = sigma.dot(epsilon)
    SEi = sympy.lambdify((xa, xb, ya, yb), SE.integrate((x1, xa, xb), (x2, ya, yb)))

    return SEi(XA, XB, YA, YB)