#!/net/python/bin/python
from numpy import *
import pylab as p

def TrapInt(f,N,xmin,xmax):
    h = (xmax-xmin)*1.0/N
    x=xmin
    sum = -0.5*( f(xmin)+f(xmax) )
    for i in arange(N+1):
        sum += f(x+h*i)

    return h*sum

def error(measured,exact):
    return abs( (measured-exact)/exact )

def tableprint(f,Nlist,xmin,xmax,exact):
    print "N\tStepSize\tI_num\t\terror"
    for N in Nlist:
        num = TrapInt(f,N,xmin,xmax)
        h = (xmax-xmin)*1.0/N
        print '%i\t' % N,
        print '%s\t\t' % round(h),
        print '%s\t\t' % round(num),
        print round(error(num,exact))

def f(x):
    return ( sin(8*x) )**4

def round(x, n=2):
    """
    rounds x to n significant digits, returns a string representation
    in scientific notation
    """
    if n < 1:
        raise ValueError("number of significant digits must be >= 1")
    return '%.*g' % (n,x)

if __name__ == "__main__":

    """
    x=arange(0,2*pi,0.01)
    y=f(x)
    p.plot(x,y)
    p.title('sin(8x)^4')
    p.xlabel('x')
    p.ylabel('f(x)')
    p.show()
    """

    Nlist = [1.0]
    while True:
        if Nlist[-1]*2 > 512:
            break
        else:
            Nlist.append(Nlist[-1]*2)
    
    tableprint(f,Nlist,0,2*pi,3*pi/4)

    Elist = []
    hlist=[]
    for N in range(1,501):
        Elist.append( error(TrapInt(f,N,0,2*pi),3*pi/4 ) )
        hlist.append(2*pi/N)

    Elist = array(Elist)
    hlist = array(hlist)
    p.loglog(hlist,Elist)
    p.xlabel('h')
    p.ylabel('Error')
    p.title('Error in integration of ( sin(8x) )^8 for step-size h')
    p.show()