"""Functions to approximate integrals.

mid_riemann       Calculate an integral using the middle Riemann sum.
left_riemann      Calculate an integral using the left Riemann sum.
right_riemann     Calculate an integral using the right Riemann sum.
trapezoidal_rule  Calculate an integral using the trapezoidal rule.
simpsons_rule     Calculate an integral using Simpson's rule.

"""
def mid_riemann(fn, a, b, n):
    """mid_riemann(fn, a, b, n) -> sum

    Calculate the middle Riemann sum of a function from a to b.  fn
    should be a function that takes a value between a and b and
    returns a number.  n is the number of divisions.

    """
    if b < a:
        raise ValueError, "a should be less than b"
    delta = float(b - a) / n        # dx
    x = a
    sum = 0.0
    while(x < b):
        midpoint = x + delta/2
        height = fn(midpoint)       # f(x)
        area = delta * height       # f(x)dx
        
        if(area >= 0):
            sum = sum + area
        else:
            sum = sum + -area

        x = x + delta
    return(sum)

def left_riemann(fn, a, b, n):
    """left_riemann(fn, a, b, n) -> sum

    Calculate the left Riemann sum of a function from a to b.  fn
    should be a function that takes a value between a and b and
    returns a number.  n is the number of divisions.

    """
    if(b < a):
        raise ValueError, "a should be less than b"
    delta = float(b - a) / n        # dx

    sum = mid_riemann(fn, a-delta/2, b-delta/2, n)
    return(sum)

def right_riemann(fn, a, b, n):
    """right_riemann(fn, a, b, n) -> sum

    Calculate the right Riemann sum of a function from a to b.  fn
    should be a function that takes a value between a and b and
    returns a number.  n is the number of divisions.

    """
    if(b < a):
        raise error, "a should be less than b"
    delta = float(b - a) / n        # dx

    sum = mid_riemann(fn, a+delta/2, b+delta/2, n)
    return(sum)


def trapezoidal_rule(fn, a, b, n):
    """trapezoidal_rule(fn, a, b, n) -> sum

    Approximate an integral using the Trapezoidal Rule.  fn should be
    a function that takes a value between a and b and returns a
    number.  n is the number of divisions.

    """
    # Trapezoidal rule is:
    # delta / 2 * (f(x0) + 2f(x1) + 2f(x2) + ... + f(xN))

    if(b < a):
        raise error, "a should be less than b"
    delta = float(b - a) / n        # dx
    sum = 0.0
    x = a
    for i in range(n+1):
        height = fn(x)      # f(x)
        if((i == 0) or (i == n)):
            sum = sum + height
        else:
            sum = sum + 2*height
        x = x + delta
    area = sum * delta / 2
    return(area)

def simpsons_rule(fn, a, b, n):
    """simpsons_rule(fn, a, b, n) -> sum

    Approximate an integral using the Simpson's Rule.  fn should be a
    function that takes a value between a and b and returns a number.
    n is the number of divisions.

    """
    # Simpson's rule is:
    # delta / 3 * (f(x0) + 4f(x1) + 2f(x2) + ... + f(xN))

    if(b < a):
        raise error, "a should be less than b"
    delta = float(b - a) / n        # dx
    sum = 0.0
    x = a
    for i in range(n+1):
        height = fn(x)      # f(x)
        if((i == 0) or (i == n)):
            sum = sum + height
        elif(_even(i)):
            sum = sum + 2*height
        else:
            sum = sum + 4*height
        x = x + delta
    area = sum * delta / 3
    return(area)

def integral_2d(fxy, x1, x2, nx, y1, y2, ny):
    """integral_2d(fxy, x1, x2, nx, y1, y2, ny) -> sum

    Calculate a 2-dimensional integral.  This must be in the form:
    INTx INTy f(x,y)dxdy
    nx is the number of divisions for x.  ny is the number of divisions
    for y.

    """
    if(x1 > x2):
        raise error, "x1 must be less than x2"
    if(y1 > y2):
        raise error, "y1 must be less than y2"

    # Do a quick Riemann sum
    dx = float(x2-x1)/nx
    dy = float(y2-y1)/ny
    volume = 0.0
    x = x1
    for i in range(nx):
        y = y1
        for j in range(ny):
            height = fxy(x, y)
            # Bug: Assumes only positive numbers
            volume = volume + height * dx * dy
            y = y + dy
        x = x + dx
    return(volume)

def _even(x):
    return(x&1 == 0)