clear all
close all
% mouvement d'un cil par equation d'advection 
%d'apres these Chatelin 
% mouvement du cil dans (x,y)

global T
T=1/10;

nbx=100;
nbt=36;

dx=1/nbx;
dt=T/nbt;
t=(1:nbt)*T/nbt;
xi=linspace(0,1,nbx);

Px=zeros(nbx,nbt);
Py=zeros(nbx,nbt);
Lx=zeros(nbx,nbt);
Ly=zeros(nbx,nbt);
B=zeros(nbx,1);
lcil=zeros(nbt,1);

% hfx=figure, hold on;
% hfy=figure, hold on;
hf=figure, hold on;

for k=1:nbt
    modelas=elas(t(k));  %fonction qui calcule nu(t)
    A=diag(-modelas*ones(nbx-1,1)*dt/dx,-1)... %diagonale inf a cause de -1
        +diag(1+modelas*ones(nbx,1)*dt/dx,0);  %diagonale a cause du zero a la fin
    A(1,1)=1;A(1,2)=0;
    InvA=inv(A);
    
    %pour la position x du cil
    if (k>1)
        B=Lx(:,k-1);
    end
    B(1,1)=2*cos(2*pi*t(k)/T);
    Lx(:,k)=InvA*B;   
    Px(:,k)=cumsum(Lx(:,k))*dx;
%     figure(hfx)
%     plot(xi,Px(:,k))
    
    %pour la la position y du cil
    if (k>1)
        B=Ly(:,k-1);
    else
        B=ones(nbx,1);
    end
    B(1,1)=1;
    Ly(:,k)=InvA*B;   
    Py(:,k)=cumsum(Ly(:,k))*dx;
%     figure(hfy)
%     plot(xi,Py(:,k))   
    
    %calcul de la longueur du cil pour renormalisation de la longueur a 1
    lcil(k)=sum((Lx(:,k).^2+Ly(:,k).^2).^0.5)*dx;
end

figure(hf);hold on
%config de reference 
hp=plot(zeros(nbx,1),xi,'b');
axis equal, grid on
axis([-1 1 0 1.1])
hx=xlabel('x');
hy=ylabel('y');
set(hp,'LineWidth',2)
F(1)=getframe(gcf);
for k=1:nbt
    hp=plot(Px(:,k)/lcil(k),Py(:,k)/lcil(k),'b');
    axis equal, grid on
    axis([-1 1 0 1.1])
    hx=xlabel('x');
    hy=ylabel('y');
    set(hp,'LineWidth',2)
    F(k+1)=getframe(gcf);
end

movie(F,10,2)