%
% POUR FAIRE UNE APPLICATION CONSOLE, transformer ce script en fonction:
% function oscillateur()
% puis supprimer toute référence au fichier AVI et aux frames (en mettant avi à 0, ligne 7)
% et enfin compiler avec:
% mcc -m -B sgl oscillateur oscill_modele oscill_entree
avi=0;
%
% POUR FAIRE UNE ANIMATION, rétablir avi à 1 (ligne 7)
% (création du fichier, ajout de frame, fermeture du fichier)
% puis transformer le fichier AVI en Animated GIF, avec Corel par exemple.
% Scenari
global scenario
% scenario est à valeurs dans { 'amorti' , 'libre' , 'entretenu' , 'stable' ,'filtre' }
scenario='test';
% Constantes physique par défaut: peuvent çetre écrasées suivant le scenario.
n=5;        % Nombre de masses
omega0=2;   % sqrt( Force de rappel / Masse )
b=0;        % Raideur / Masse: système conservatif
tmax=100;   % durée de la simulation
ns=5;       % nombre de spires (pour le dessin)
% les conditions initiales sont stockées dans x0 et dépendent de la simulation
% (les n vitesses puis les n positions, la condition initiale statique est donc [0,0,...,0,1,2,...,n])
% La durée de la simulation est dans tmax

%
% Constantes du modèle
% --------------------
x0=[1:n,zeros(1,n)]';
switch lower(scenario),
case 'test',
    disp('Oscillateur harmonique, mode test'),
    n=3;
    x0=[1.5,2,2.5,zeros(1,n)]';
    omega0=1/2;
case 'amorti',
    disp('Oscillateur harmonique simple'),
    n=1;
    omega0=1;
    b=1/8;
    x0=[2,0];
    ns=20;
    Amp=(x0(1)-1)*sqrt(1+(b/2/omega0)^2);
case 'libre',
    disp('Oscillateur harmonique libre'),
    x0(n+1)=1;
    x0(2*n)=-1;
    tmax=50;
case 'entretenu',
    disp('Oscillateur harmonique entretenu'),
    tmax=200;
case 'stable',
    disp('Oscillateur harmonique stabilisé'),
    x0(n+1:n+n)=x0(n+1:n+n)+randn(n,1)*0.8;
case 'filtre',
    disp('Oscillateur harmonique comme filtre'),
    tmax=200;
    x0=[0.9574 1.9047 2.8814 3.8938 4.9378 0.1489 0.2178 0.2346 0.1961 0.1114];
    b=1/10;
otherwise
    disp('Scenario inconnu'),
end,
% A = Matrice du système
A=-2*diag(ones(1,n)) + diag(ones(n-1,1),1) + diag(ones(n-1,1),-1);
A=[zeros(n,n),eye(n);omega0^2*A,-b*eye(n)];
B=zeros(2*n,1);
B(n+1)=1;
C=zeros(2*n,1);
C(n+n)=omega0^2*(n+1);
%
% Simulation
% ----------
%options = odeset('OutputFcn',@sorties);
[t,x] = ode23('oscill_modele',[0,tmax],x0,[],A,B,C);
% recalcul de l'entrée car je n'arrive pas à faire marcher OutputFcn !!!
u=zeros(length(t),1);
for k=1:length(t),
    u(k)=oscill_entree(t(k),x(k,:));
%    if (k>10) & (abs(u(k)-1/5)<1e-5), t(k),u(k),x(length(t),:), end,
end;
x=x(:,1:n);
%
% Animation, le résultat est une animation au format AVI
% ------------------------------------------------------
fig=figure(1);
set(fig,'DoubleBuffer','on');
ct=tmax/50;
cx=(n+1)/50;
cu=cx/4;
color='rgbymc';
if avi,
    aviobj = avifile('oscillateur.avi','compression','none','quality',100),
end,
lk=length(t);
dk=round(length(t)/100);
for k=1:dk:lk,
    set(gca,'nextplot','replace')
    plot(t(1:k),u(1:k),'k'),
    axis([0 tmax -1 n+1]),
    set(gca,'nextplot','add')
    for i=1:n,
        plot(t(1:k),x(1:k,i),color(1+mod(i-1,6))),
    end,
    tc=t(k);
    tp=[tc-ct,tc+ct,tc+ct,tc-ct,tc-ct];
    up=[u(k)-cu,u(k)-cu,u(k)+cu,u(k)+cu,u(k)-cu];
    patch(tp,up,'k');    
    T=tc+ct/2*(-1).^[1:2*ns+2];
    xc=[x(k,:),n+1+cx];
    d=(xc(1)-u(k)-cx-cu)/(2*ns+1);
    X=u(k)+cu+[0:2*ns+1]*d;
    plot(T,X,'k-');hold on,
    for i=1:n,
        xp=[xc(i)-cx,xc(i)-cx,xc(i)+cx,xc(i)+cx,xc(i)-cx];
        patch(tp,xp,color(1+mod(i-1,6)));
        d=(xc(i+1)-xc(i)-2*cx)/(2*ns+1);
        X=xc(i)+cx+[0:2*ns+1]*d;
        plot(T,X,'k-');
    end,
    if n==1,
        plot(t(1:k),1+Amp*exp(-(b/2)*t(1:k)),'c:'),
        plot(t(1:k),1-Amp*exp(-(b/2)*t(1:k)),'c:'),
    end,
    if avi,
        frame = getframe(gca);
        aviobj = addframe(aviobj,frame);
    end,
    drawnow,
end,
if avi,
    aviobj = close(aviobj);
end,
