import numpy as np
import matplotlib.pylab as plt

g = 9.81

def integrate(f, a, b, N=1000):
    xn = np.linspace(a, b, N)
    summe = 0.0
    dx = abs((a-b)/(N-1))
    for i in range(1,N):
        area = ( f(xn[i]) + f(xn[i-1]) ) / 2.0 * dx
        summe += area
    return summe

def T(theta, l):
    K = lambda x: 1.0 / np.sqrt(1.0 - np.sin(theta/2.0)**2 * x)
    return 4.0 * np.sqrt(l/g)*integrate(K, 0.0, np.pi/2.0)
    
def calc_T(l=1.0, theta_max=90.0, N=50):
    angles = np.linspace(0.0, theta_max, N)
    time = []
    for theta in (angles/180.0*np.pi):
        time.append(T(theta, l))
    return angles, time

l = 1.0
theta_max = 90.0
N = 50

angles, numeric = calc_T(l, theta_max, N)
analytic = np.ones(N)*2.0*np.pi*np.sqrt(l/g)

plt.figure()
plt.plot(angles, numeric, 'o-')
plt.plot(angles, analytic, '-')
plt.legend(('numerisch', 'analytisch'))
plt.show()