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

#Romberg integration code adapted from Wikipedia
# http://en.wikipedia.org/wiki/Romberg's_method

def Romberg(f, N, xmin, xmax, PRINT = False):
    """Approximate the definite integral of f from a to b by Romberg's method.
    eps is the desired accuracy."""
    R = [[0.5 * (xmax - xmin) * (f(xmin) + f(xmax))]]  # R[0][0]
    if PRINT:
        print ' '.join('%11.4f' % x for x in R[0])
    for i in range(1,N+1):
        h = float(xmax-xmin)/2**i
        R.append((i+1)*[None])  # Add an empty row.
        R[i][0] = 0.5*R[i-1][0] + \
                  h*sum(f(xmin+(2*k-1)*h) for k in range(1, 2**(i-1)+1))
        for m in range(1, i+1):
            R[i][m] = R[i][m-1] + (R[i][m-1] - R[i-1][m-1]) / (4**m - 1)
        if PRINT:
            print ' '.join('%11.4f' % x for x in R[i])
        last = R[i][i]
    return last

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 = Romberg(f,N,xmin,xmax)
        h = (xmax-xmin)*1.0/(2**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__":

    Romberg(f,9,0,2*pi,True)

    Nlist = range(1,10)
    
    tableprint(f,Nlist,0,2*pi,3*pi/4)

    Elist = []
    hlist=[]
    for N in Nlist:
        Elist.append( error(Romberg(f,N,0,2*pi),3*pi/4 ) )
        hlist.append(2*pi/(2**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()