from numpy import *

def simpson(x, f, F=[]):
    """
    Use composite Simpsons Rule quadrature based on the subdivision xi
    Input:
    xi: location of nodes, defining the subdivision of the interval
    (including the endpoints)
    f:  function to integrate
    F:  antiderivative (for computing errors)
    
    
    Output:
    
    I = the numerical approximation to the integral on each subinterval
    E = error (absolute value) on each subinterval (only if F is known)
    
    """
    n = len(x)         # number of nodes
    assert(n % 2 == 1) # should be an odd number
    h = x[1]-x[0]
    # according to the book's definition, 1 panel in Simpson = 2h
    # number of panels
    m = int((n-1)/2)
    I = zeros(m)
    E = zeros(m)
    for i in range(m):
        I_simpson = h/3*(f(x[2*i])+4*f(x[2*i+1])+f(x[2*i+2]))
        if F:
            I_exact = F(x[2*i+2])-F(x[2*i])
            E[i] = abs(I_simpson-I_exact)
        I[i] = I_simpson
    return I, E

# Test Simpsons Rule
if __name__ == "__main__":
    # Test
    N = 201
    a = 0
    b = 1
    x = linspace(a,b,N)
    f = lambda x: sqrt(x)
    F = lambda x: x*sqrt(x)*2./3.
    I, E = simpson(x, f, F)
    Iq= sum(I)
    print("N = %6d, numerical = %e, error = %e" % (N,Iq,F(b)-F(a)-Iq))


    m = 32
    N = 2*32+1

    # (a)
    a,b=0,1
    x = linspace(a,b,N)
    f = lambda x: x**3
    fprim = lambda x: 3*x*x
    g = lambda x: sqrt(1+fprim(x)**2)
    I,_=simpson(x, g)
    print('(a):  arclength ~ %e' % sum(I))
    
    # (b)
    a,b=0,pi/4
    x = linspace(a,b,N)
    f = lambda x: tan(x)
    fprim = lambda x: 1/cos(x)**2
    g = lambda x: sqrt(1+fprim(x)**2)
    I,_=simpson(x, g)
    print('(b):  arclength ~ %e' % sum(I))

    # (c)
    a,b=0,1
    x = linspace(a,b,N)
    f = lambda x: atan(x)
    fprim = lambda x: 1/(1+x*x)
    g = lambda x: sqrt(1+fprim(x)**2)
    I,_=simpson(x, g)
    print('(c):  arclength ~ %e' % sum(I))