This is used to calculate an elliptical sector area for a numerical analysis problem...

area problema de la elipse

from math import sqrt, cos

MAX = 1.52097701e8 #km
MIN = 1.47098074e8 #km
A = (MAX + MIN)/2.
C = (MAX - MIN)/2.;
E = C/A
B = sqrt(A**2-C**2)
#theta = var('theta')
def R(theta): return A*(1-E**2)/(1 - E*cos(theta))

def trapezoid(f,a,b,Iold,k):
    if k == 1: Inew = (f(a) + f(b))*(b - a)/2.0
    else:
        n = 2**(k -2 ) # Number of new points
        h = (b - a)/n # Spacing of new points
        x = a + h/2.0
        sum = 0.0
        for i in range(n):
            sum = sum + f(x)
            x = x + h
        Inew = (Iold + h*sum)/2.0
    return Inew

## module romberg
''' I,nPanels = romberg(f,a,b,tol=1.0e-15).
    Romberg intergration of f(x) from x = a to b.
    Returns the integral and the number of panels used.
'''
from numpy import zeros

def romberg(f,a,b,tol=1.0e-6):
    def richardson(r,k):
        for j in range(k-1,0,-1):
            const = 4.0**(k-j)
            r[j] = (const*r[j+1] - r[j])/(const - 1.0)
        return r
    r = zeros(21)
    r[1] = trapezoid(f,a,b,0.0,1)
    r_old = r[1]
    for k in range(2,21):
        r[k] = trapezoid(f,a,b,r[k-1],k)
        r = richardson(r,k)
        
        if abs(r[1]-r_old) < tol*max(abs(r[1]),1.0):
            return r[1],2**(k-1)
        r_old = r[1]
    print "Romberg quadrature did not converge"

def f(x): return R(x)**2
print romberg(f,0.,1.)[0]