import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

g = 9.81
l = 0.5


def f(y, t):
    x, vx = y
    return [vx, (-1) * g / l * (x - (x ** 3) / 6)]


t = np.linspace(0, 50, 10000)

X0 = [0, 0]
x, vx = odeint(f, X0, t, full_output=False).T


Amplituda=[]
Chastota=[]


'''
def AH(nach, konec):
    N = 100
    dx = (konec - nach) / N
    fi = nach
    while fi < konec:
        t = np.linspace(0, 50, 10000)
        X0 = [fi, 0]
        x, vx = odeint(f, X0, t, full_output=False).T
        j = 0
        t_posled = 0
        Amplituda1=[]
        Chastota1=[]
        for i in xrange(0, 10000):
            if i > 0 and x[i - 1] * x[i] < 0:
                j += 1
                t_posled = 50. / 10000 * i
        np.append.Amplituda1 = fi
        np.append.Chastota1 = t_posled / j * 2
        fi += dx
        return Amplituda1, Chastota1
'''
N = 100
konec=math.pi / 2
nach=math.pi / 100
dx = (konec - nach) / N
fi = nach
Amplituda1 = []
Chastota1 = []
while fi < konec:
    X0 = [fi, 0]
    x, vx = odeint(f, X0, t, full_output=False).T
    j = 0
    t_posled = 0
    for i in xrange(0, 10000):
        if i > 0 and x[i - 1] * x[i] < 0:
            j += 1
            t_posled = 50. / 10000 * i
    np.append.Amplituda1 = fi
    np.append.Chastota1 = j / t_posled / 2
    fi += dx


#Amplituda, Chastota = AH(math.pi / 100, math.pi / 2)

fig = plt.figure(facecolor='white')
plt.plot(Chastota1, Amplituda1)
plt.grid(True)
plt.show()