# -*- coding: utf-8 -*-
# import basic modules
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
import pylab as pyl
import math

NPT=700
NPampli=300 # nombre de points pour T et a0
a0=.9999
omega=2.*math.pi
T=4. # valeur max de a0, pulsation et periode max
x=np.linspace(1,2,NPT)-1
y=np.linspace(1,2,NPampli)-1 #%discetisation en T et a0
z=np.zeros((NPT,NPampli),float) # matrice des resultats
for j in range(NPampli):# boucle sur  les valeurs de l'amplitude
   e=math.sqrt(1-a0*y[j])#% coefficient d'attenuation
   omegap=omega*e #% fequence perturbee
   for i in range(NPT):  #boucle sur les periodes
       xi=.5*(omega/omegap+omegap/omega)
       phase1=omega*T*x[i]/2
       phase2=omegap*T*x[i]/2
       q=np.max(np.abs(math.cos(phase2)*math.cos(phase1)-xi*math.sin(phase1)*math.sin(phase2))-1,0)
       if (q>0):
                z[i,j]=1 # cas instable
       else:
              z[i,j]=0 # cas stable
xx=T*x #axe des periodes
yy=a0*y # axe des amplitudes

fig = plt.figure()
ax = fig.gca()
Y, X = pyl.meshgrid(yy, xx)

pyl.contourf(X, Y, z, 8, alpha=.75, cmap='jet')
#C = pyl.contour(X, Y, z, 8, colors='black', linewidth=.5)
ax.set_ylabel('Amplitude de la perturbation')
ax.set_xlabel('T')


plt.title(u'Zone d\u0027instabilites de l\u0027équation de Hill')

plt.show()