def EI(sections, normal=None):  # {{{
    """Calculate the bending stiffnes of a cross-section.

    The cross-section is composed out of rectangular nonoverlapping sections
    that can have different Young's moduli.

    Each section is represented by a 4-tuple (width, height, offset, E).
    The offset is the distance from the top of the section to the top of the
    highest section. This should always be a positive value.
    E is the Young's modulus of the material of this section.

    Arguments:
        sections: Iterable of section properties.
        normal: The Young's modulus to which the total cross-section will be
            normalized. (Not used anymore, retained for compatibility.)

    Returns:
        Tuple of EI, top and bottom. Top and bottom are with respect to the
        neutral line.

    Examples:
        >>> E = 210000
        >>> B = 100
        >>> H = 20
        >>> sections = ((B, H, 0, E),)
        >>> EI(sections)
        (14000000000.0, 10.0, -10.0)

        >>> B = 100
        >>> h = 18
        >>> t = 1
        >>> H = h + 2 * t
        >>> E = 210000
        >>> sections = ((B, t, 0, E), (B, t, h+t, E))
        >>> EI(sections)
        (3794000000.0, 10.0, -10.0)

        >>> E1, E2 = 200000, 71000
        >>> t1, t2 = 1.5, 2.5
        >>> H = 31
        >>> B = 100
        >>> sections = ((B, t1, 0, E1), (B, t2, H-t2, E2))
        >>> EI(sections)
        (9393560891.143106, 11.530104712041885, -19.469895287958117)
    """
    normal = sections[0][-1]
    normalized = tuple((w * E / normal, h, offs) for w, h, offs, E in sections)
    A = sum(w * h for w, h, _ in normalized)
    S = sum(w * h * (offs + h / 2) for w, h, offs in normalized)
    yn = S / A
    # Find any geometry that straddles yn.
    to_split = tuple(g for g in sections if g[2] < yn and g[1] + g[2] > yn)
    geom = tuple(g for g in sections if g not in to_split)
    # split that geometry.
    # The new tuple has the format (width, height, top, bottom)
    new_geom = []
    for w, h, offs, E in to_split:
        h1 = yn - offs
        h2 = h - h1
        new_geom.append((w, h1, h1, 0, E))
        new_geom.append((w, h2, 0, -h2, E))
    # Convert the remaining geometry to reference yn.
    for w, h, offs, E in geom:
        new_geom.append((w, h, yn - offs, yn - offs - h, E))
    EI = sum(E * w * (top ** 3 - bot ** 3) / 3 for w, h, top, bot, E in new_geom)
    top = max(g[-3] for g in new_geom)
    bot = min(g[-2] for g in new_geom)
    return EI, top, bot