clear all
 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=sqrt(k/M);% pusation de pilonnement
 omega2=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=linspace(1,2,NPQ);
 yq=linspace(1,2,NPQ);
 x=linspace(1,2,NPV);
 x2=linspace(1,2,NPV);
 y=linspace(1,2,NPV);
 omegar1=0.*linspace(1,2,NPV);
 amor1=0.*linspace(1,2,NPV);
 omegar2=0.*linspace(1,2,NPV);
 amor2=0.*linspace(1,2,NPV);
 %=== boucle sur Q=partialcm/partialalpha 
 %=== pour trouver les zones ou petit delta>0
 for ii=1:NPQ
   Q=Qmin+(ii-1)*Qmax/(NPQ-1);
   xq(ii)=Q;
   c1=(omega1^2)*(1-xi);
   c2=a*R*(omega1^2*xi-omega2^2);
   c3=a^2*R^2*omega2^2;
   yq(ii)=Q^2*c1+Q*c2+c3;
 end
 %== boucle sur les vitesses pour le calcul du discriminant grand delta,
 Q=QQ;
 for ii=1:NPV
     V=Vmin+(ii-1)*Vmax/(NPV-1);
     x(ii)=V;
     x2(ii)=V^2;
     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;
     y(ii)=A*V^4+2*B*V^2+C;
     % calcul des racines en lambda carre
         FA=1;
         FB=-((omega1^2)*xi+(omega2^2)+x2(ii)*(a*R-Q)/JG);
         FC=(omega1^2)*(omega2^2-Q * x2(ii)/JG);
         omegar1(ii)=real((-FB+sqrt(y(ii)))/(2*FA));
         amor1(ii)=imag((-FB+sqrt(y(ii)))/(2*FA));
         omegar2(ii)=real((-FB-sqrt(y(ii)))/(2*FA)); 
         amor2(ii)=imag((-FB-sqrt(y(ii)))/(2*FA));
 end
figure 
subplot(2,2,1)
plot(xq,yq)
xlabel('coefficient q');ylabel('petit delta')
title('Variation du signe de petit delta')
subplot(2,2,2)
plot(x2,y)
xlabel('Vitesse V au carre');ylabel('Discriminant Deta')
title('Evolution du discriminant')
subplot(2,2,3)
plot(x,omegar1,x,omegar2)
xlabel('Vitesse');ylabel('Evolution des deux pulsations propres')
title('Courbe annoncant le flutter')
subplot(2,2,4)
plot(x,amor1,x,amor2)
xlabel('Vitesse');ylabel('Evolution des amortissements aero')
title('Courbe des amortissements')