Method to calculate the surface of the intersection of a rectangle (aligned with axes) and another surface
    defined by two functions f_lower and f_upper.

    Args:
        rectangle:
            Rectangle defined as : ((x1, x2), (y1, y2)).
        f_lower:
            Function defining the lower bound of the surface.
        f_upper:
            Function defining the upper bound of the surface.
        bounds_lower:
            Interval in which the f_lower function is defined.
        bounds_upper:
            Interval in which the f_upper function is defined.
        check:
            Whether to check if f_lower is always lower than f_upper.
        numpoints_check:
            Number of points used to check whether f_lower is always lower than f_upper

    Returns:
        The surface of the intersection of the rectangle and the surface defined by f_lower and f_upper.

def rectangle_surface_intersection(rectangle, f_lower, f_upper,
                                   bounds_lower=None, bounds_upper=None,
                                   check=True, numpoints_check=500):
    """
    Method to calculate the surface of the intersection of a rectangle (aligned with axes) and another surface
    defined by two functions f_lower and f_upper.

    Args:
        rectangle:
            Rectangle defined as : ((x1, x2), (y1, y2)).
        f_lower:
            Function defining the lower bound of the surface.
        f_upper:
            Function defining the upper bound of the surface.
        bounds_lower:
            Interval in which the f_lower function is defined.
        bounds_upper:
            Interval in which the f_upper function is defined.
        check:
            Whether to check if f_lower is always lower than f_upper.
        numpoints_check:
            Number of points used to check whether f_lower is always lower than f_upper

    Returns:
        The surface of the intersection of the rectangle and the surface defined by f_lower and f_upper.
    """
    x1 = np.min(rectangle[0])
    x2 = np.max(rectangle[0])
    y1 = np.min(rectangle[1])
    y2 = np.max(rectangle[1])
    # Check that f_lower is allways lower than f_upper between x1 and x2 if no bounds are given or between the bounds
    #  of the f_lower and f_upper functions if they are given.
    if check:
        if bounds_lower is not None:
            if bounds_upper is not None:
                if not all(np.array(bounds_lower) == np.array(bounds_upper)):
                    raise ValueError('Bounds should be identical for both f_lower and f_upper')
                if not '<' in function_comparison(f1=f_lower, f2=f_upper,
                                                  x1=bounds_lower[0], x2=bounds_lower[1],
                                                  numpoints_check=numpoints_check):
                    raise RuntimeError('Function f_lower is not allways lower or equal to function f_upper within '
                                       'the domain defined by the functions bounds.')
            else:
                raise ValueError('Bounds are given for f_lower but not for f_upper')
        elif bounds_upper is not None:
            if bounds_lower is None:
                raise ValueError('Bounds are given for f_upper but not for f_lower')
            else:
                if not '<' in function_comparison(f1=f_lower, f2=f_upper,
                                                  x1=bounds_lower[0], x2=bounds_lower[1],
                                                  numpoints_check=numpoints_check):
                    raise RuntimeError('Function f_lower is not allways lower or equal to function f_upper within '
                                       'the domain defined by the functions bounds.')
        else:
            if not '<' in function_comparison(f1=f_lower, f2=f_upper,
                                              x1=x1, x2=x2, numpoints_check=numpoints_check):
                raise RuntimeError('Function f_lower is not allways lower or equal to function f_upper within '
                                   'the domain defined by x1 and x2.')
    if bounds_lower is None:
        raise NotImplementedError('Bounds should be given right now ...')
    else:
        if x2 < bounds_lower[0] or x1 > bounds_lower[1]:
            return (0.0, 0.0)
        if x1 < bounds_lower[0]:
            xmin = bounds_lower[0]
        else:
            xmin = x1
        if x2 > bounds_lower[1]:
            xmax = bounds_lower[1]
        else:
            xmax = x2
        def diff(x):
            flwx = f_lower(x)
            fupx = f_upper(x)
            minup = np.min([fupx, y2 * np.ones_like(fupx)], axis=0)
            maxlw = np.max([flwx, y1 * np.ones_like(flwx)], axis=0)
            zeros = np.zeros_like(fupx)
            upper = np.where(y2 >= flwx, np.where(y1 <= fupx, minup, zeros), zeros)
            lower = np.where(y1 <= fupx, np.where(y2 >= flwx, maxlw, zeros), zeros)
            return upper-lower
        return quad(diff, xmin, xmax)