%clear all;
NPT=1000; NPampli=1000;% nombre de points pour T et a0
a0=.999;omega=2.*pi;T=4.; % valeur max de a0, pulsation et periode max
x=linspace(1,2,NPT)-1;y=linspace(1,2,NPampli)-1;%discetisation en T et a0
z=zeros(NPT,NPampli);% matrice des resultats
  for j=1:NPampli% boucle sur  les valeurs de l'amplitude
    e=sqrt(1-a0*y(j));% coefficient d'attenuation
         omegap=omega*e;% fequence perturbee
    for i=1:NPT  %boucle sur les periodes
        xi=.5*(omega/omegap+omegap/omega);
         phase1=omega*T*x(i)/2;
           phase2=omegap*T*x(i)/2;
        q=max(abs(cos(phase2)*cos(phase1)-xi*sin(phase1)*sin(phase2))-1,0);
        if(q >0) 
            z(i,j)=1;% cas instable
        else
            z(i,j)=0;% cas stable
        end
    end
  end
xx=T*x; %axe des periodes
yy=a0*y;% axe des amplitudes
figure
surfc(yy,xx,z);
shading interp
view([90,-90])
ylabel('T');
    xlabel('amplitude de la perturbation');
    zlabel('zones instables en jaunes');
    title('Zone d instabilites de l equation de Hill')