# -*- coding: utf-8 -*-
"""
Created on Wed Jan  1 14:30:33 2015
ok
@author: joseorellana
"""

# 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 cmath 

a=-.05
JG=1
M=1000.
k=900.
c=16.
J0=JG+(a**2)*M 
xi=J0/JG
R=10.
Qmin=0. # plage de recherche de Q borne inf
Qmax=5. #plage de recherche de Q  borne sup
omega1=np.sqrt(k/M) # pusation de pilonnement
omega2=np.sqrt(c/JG) # pulsation de torsion
NPV=200
NPQ=100 # choix du nombre de points pour les courbes
QQ=.1 # choix d'une valeur de Q
Vmin=0. #vitesse minimale
Vmax=10. # vitesse maximale
#===== declaration des vecteurs de travail
xq=np.linspace(1,2,NPQ)
yq=np.linspace(1,2,NPQ)
x=np.linspace(1,2,NPV)
x2=np.linspace(1,2,NPV)
y=np.linspace(1,2,NPV)
omegar1=0.*np.linspace(1,2,NPV)
amor1=0.*np.linspace(1,2,NPV)
omegar2=0.*np.linspace(1,2,NPV)
amor2=0.*np.linspace(1,2,NPV)
c1=(omega1**2)*(1-xi)
c2=a*R*(omega1**2*xi-omega2**2)
c3=a**2*R**2*omega2**2
 #=== boucle sur Q=partialcm/partialalpha 
 #=== pour trouver les zones ou petit delta>0
for ii in range(NPQ):
   Q=Qmin+ii*Qmax/(NPQ-1)
   xq[ii]=Q
   yq[ii]=Q**2*c1+Q*c2+c3
#== boucle sur les vitesses pour le calcul du discriminant grand delta,
Q=QQ;
A=((a*R-Q)**2)/(JG**2)
B=2*(omega1**2)*Q/JG+(xi*omega1**2+omega2**2)*(a*R-Q)/JG
C=(xi*omega1**2+omega2**2)**2-4*omega1**2*omega2**2
FA=1
for ii in range(NPV):
     V=Vmin+ii*Vmax/(NPV-1)
     x[ii]=V
     x2[ii]=V**2
     y[ii]=A*V**4+2*B*V**2+C
     # calcul des racines en lambda carre
     FB=-((omega1**2)*xi+(omega2**2)+x2[ii]*(a*R-Q)/JG)
     FC=(omega1**2)*(omega2**2-Q * x2[ii]/JG)
     omegar1[ii]=((-FB+cmath.sqrt(y[ii]))/(2*FA)).real
     amor1[ii]=((-FB+cmath.sqrt(y[ii]))/(2*FA)).imag
     omegar2[ii]=((-FB-cmath.sqrt(y[ii]))/(2*FA)).real
     amor2[ii]=((-FB-cmath.sqrt(y[ii]))/(2*FA)).imag




plt.subplot(2, 2, 1)
plt.plot(xq, yq, '-')
plt.xlabel('coefficient q')
plt.ylabel('petit delta')
plt.title('Variation du signe de petit delta')

plt.subplot(2, 2, 2)
plt.plot(x2, y, '-')
plt.xlabel('Vitesse V au carre')
plt.ylabel('Discriminant Deta')
plt.title('Evolution du discriminant')

plt.subplot(2,2,3)
plt.plot(x,omegar1,x,omegar2)
plt.xlabel('Vitesse')
plt.ylabel('Evolution des deux pulsations propres')
plt.title('Courbe annoncant le flutter')

plt.subplot(2,2,4)
plt.plot(x,amor1,x,amor2)
plt.xlabel('Vitesse')
plt.ylabel('Evolution des amortissements aero')
plt.title('Courbe des amortissements')


plt.show()