def linear_stress_block(segments, q1, q1_y, q2, q2_y):
    """
    A function to calculate P,Mx, and My by line
    integral along given line segments for a linearly
    varing stress distribution

    Parameters
    ----------
    segments: List of two tuples/lists of two floats
                each segment should be of the form [[x1,y1],[x2,y2]]
                each x,y should be a float. y1 or y2 should exactly match
                q1_y or q2_y.
                example input:
                [[[x11,y11],[x21,y21]],...,[[x1i,y1i],[x2i,y2i]]]

    q1: float
        stress value 1 for region, usually the start stress
    q1_y: float
        y-coordinate/elevation where q1 applies,
        used to swap q1 and q2 when needed.
    q2: float
        stress value 2 for region, usually the end stress
    q2_y: float
        y-coordinate/elevation where q2 applies,
        used to swap q1 and q2 when needed.

    Returns:
    ---------
    P: float
        Sum of Axial force from all segments
    Mx: float
        Sum of Moments about the x-axis from all segments
    My: float
        Sum of Moments about the y-axis from all segments
    x: float
        x centroid coordinate of P action
    y: float
        y centroid coordinate of P action
    details: list of floats
            list of P,Mx,My values per segment
    """

    P = 0
    Mx = 0
    My = 0
    x = 0
    y = 0
    details = []

    for s in segments:
        x1 = s[0][0]
        y1 = s[0][1]
        x2 = s[1][0]
        y2 = s[1][1]

        # set start and end stress based on y coordinate of segment end point
        if y2 == q2_y:
            qs = q1
            qe = q2
        else:
            qs = q2
            qe = q1

        axial = (1/6.0)*(y2-y1)*((qs*((2*x1)+x2))+(qe*(x1+(2*x2))))
        P += axial

        momentx = ((1/12.0)*(y2-y1)
                    * (
                        qs*x1*((3*y1)+y2)
                        + qs*x2*(y1+y2)
                        + qe*x1*(y1+y2)
                        + qe*x2*(y1+(3*y2))
                        )
                    )
        Mx += momentx

        momenty = ((-1/24.0)*(y1-y2)
                    * (
                        x1*x1*((3*qs)+qe)
                        + 2*x1*x2*(qs+qe)
                        + x2*x2*(qs+(3*qe))
                        )
                    )
        My += momenty

        details.append([axial,momentx,momenty])

    x = My/P
    y = Mx/P

    return P,Mx,My,[x,y],details