%clear all
V=210;
J0=.1;
L=1;
S=1;
a=.02;
ro=1.2;
T=5;
NPT=100*256;
alpha=linspace(1,2,NPT);
alphadot=linspace(1,2,NPT);
temps=linspace(1,T+1,NPT)-1;
dt=T/NPT;
q=ro*S*L/(2*J0);
omega=2*pi*20;
omega2=omega^2;
amor=.00009*omega;
alpha0=.6;%2*pi*(20/180);
alpha1=0.;
alpha(1)=alpha0;
alpha(2)=alpha(1)+dt*alpha1;
alphadot(1)=alpha1;
alphadot(2)=alphadot(1);
for i=3:NPT
    
    vap2=V^2-2*V*alphadot(i-1)*a*sin(alpha(i-1))+(a*alphadot(i-1))^2;
    alphap=atan(tan(alpha(i-1))-a*alphadot(i-1)/(V*cos(alpha(i-1))));
    %
    alpha(i)=2*alpha(i-1)-alpha(i-2)-amor*alphadot(i-1)*dt-(dt^2)*omega2*(alpha(i-1)-alpha0)+dt^2*(q*vap2)*cm0raf(alphap);
    alphadot(i)=(alpha(i)-alpha(i-1))/dt;
end
figure
plot(alpha,alphadot)
title('Solution dans le plan des phases')
xlabel('alpha')
ylabel('dot alpha')