%  test7.m    demo3 pour l X        demo3_2 pour octave 6.4 
%  advection avec conditions limites periodiques
%  schema explicite  decentre ;  profil regulier  , 
%  fd,  4 janvier 2K, 21 fevrier 2K, 15 janvier 2002, 17 janvier 2005.
%      05 fevrier 2006, fevrier 2025 .
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cfl=.8;  tmax=1. ;  n1=101;   dx=1/(n1-1);
%%%%   [chcfl,msg] = sprintf (' nombre de Courant : %5.2f',cfl);
%%%%  [chn1,msg] = sprintf  ('        nombre de mailles : %3i',n1-1) ; 
%%%%  [chtps,msg] = sprintf ('        temps final  : %5.2f',tmax);chtps ;
%%%%  titx = [chcfl , chn1] ; titx ;
%%%%  titx1 = [chcfl , chn1 , chtps] ; titx1 ;
x=0:dx:1 ;   n=length(x) ; 
%              profil regulier de classe C2 comme condition initiale
pdebut = [16*314  -4*794  539  0 0 0 ] ; 
pfin   = [16*350   -4*878 589  0 0 0 ] ; 
p      = [1 -1 0 1 ] ;
u0=zeros(n,1);      u0(2:n1/2-1,1)=1;  
for j=1:n-1 ; xx = (j-1)*dx ;
if xx<0.25     ; u0(j:j,1)= polyval(pdebut,xx);
elseif xx>0.75 ; u0(j:j,1)= polyval(pfin,1-xx); 
else  u0(j:j,1)= polyval(p,xx-0.5) ;  end ;  end  ; u0(n:n,1)=u0(1:1,1) ;
%%% xgset yrange[-0.2:1.2] ; gset xrange[-0: 1.0] ; 
%%%  axis ([0 1 -.2 1.2]) 
%%%  title (' Advection d un profil regulier      Conditions limites periodiques   \
%%% Schema  explicite  decentre ') ;   
%%% xlabel (titx ) ;
plot (x,u0,'*;condition initiale;') 
pause (1.5)
%%% xlabel (titx1 ) ;
dt=cfl*dx;    t=0:dt:tmax;   m=length(t)-1 ;
ue=zeros(n,1) ; ue(1:n,1)=u0(1:n) ; 
V=ue  ; U=V ;  W=V ;  
%  debut de la boucle en temps
for i=1:m ; 
W(2:n,1)   = U(2:n,1) - cfl * (U(2:n,1)-U(1:n-1,1)) ;
W(1:1,1)   = W(n:n,1) ; V=W ;   U=V ;
%   solution exacte 
    for j=1:n ; zz =  (j*dx- i*dt) / dx ; mm = round (zz) ;
    if mm<=0 ;    p=mm+n-1; mu=mm ; 
    while p<=0 ; pp=mm ;  p=pp+n-1 ; mm=p ; end ; mm=p; 
    zz = zz + mm-mu ; end ; 
    if mm==1 ;  mmm1=n-1; mmp1 = 2 ; else ; if mm==n; mmm1=n-1; mmp1=2 ;
    else mmm1=mm-1 ; mmp1 = mm+1 ; end ; end ; 
    if mm <= zz ; ue(j:j,1) = (1-zz+mm)*u0(mm)+(zz-mm)*u0(mmp1);
    else  ue(j:j,1) = (mm-zz)*u0(mmm1)+ (1+zz-mm)*u0(mm) ; end ;  end;  
%   fin de la solution exacte 
%%         plot(x,ue(:,1),"ob;;" , x,U(:,1),'*;;' )
end ;

plot(x,ue(:,1),"b;solution exacte;", x,U(:,1),'*;solution approchee;',x,U(:,1),"k;;") ; 
title('schema  explicite  decentre')
axis ([0 1 -.01  1.01])
print -djpg figure-3.jpg
print -dpdf figure-3.pdf

%
%   gset terminal postscript ;  gset output "octave/explicite.ps"; 
%  plot(x,ue(:,1),"b;solution exacte;", x,U(:,1),'*;solution approchee;',x,U(:,1),"r;;")
