%   acous2.m
%   Acoustique numerique  
%   Logiciel 'HaWAY2' avec Octave-Matlab.
%   fd    haway1  1 fevrier 2001,      tp1     16 octobre 2001.
%         haway2  06   fevrier  2001,  tp2     20  novembre 2001.
%         acous2  02   fevrier  2006, 25 fevrier 2025.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              Haway2. Test du mode d ordre un.
%   nombre de mailles ne 
%   nombre de pas de temps nt
%   nombre de sommets
%   pas d espace
    ne    = 50
    niter = 100
ns    = ne + 1 ; 
dx    = 1. / ne ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   condition initiale  
%   jinit   parametre relatif a la condition initiale
%   kmode   numero du mode             
%   jinit = 0  champ de pression nul
%         = 1  vitesse nulle,  pression selon un mode en cosinus
%         = 2  vitesse nulle,  pression selon un mode en sinus
jinit = 2 ; 
kmode = 2 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   conditions aux limites
%   jclgau  type de condition limite a gauche 
%   jcldro  type de condition limite a droite 
%   si    jcl*** = 1  : Dirichlet homogene en vitesse
%                = 2  : Condition d impedance ; pp = ZZ * uu 
%   zzgau valeur de l impedence a gauche   
%   zzdro valeur de l impedence a droite   
jclgau = 2 ; 
jcldro = 2 ; 
zzgau  = -1 ; 
zzdro  = 0 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   jdt type de calcul du pas de temps 
%   jdt  = 1   calcul automatique avec condition cfl
%        = 2   pas de temps impose
jdt  = 1  ;
    if (jdt == 1) ; 
        cfl  = 0.6 
    dt   = cfl * dx ; 
    else ;
        dt   = 0.1
        cfl  = dt / dx ; 
    end ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                  abscisses des sommets
xx        = 0 ; 
x         = zeros(ns,1) ; 
for     j = 1:ns ; 
x(j:j,1)  = xx ; 
xx        = xx + dx ; 
end ;
%                                  milieux des elements
xx        = 0.5 * dx ; 
y         = zeros(ne,1) ; 
for     k = 1:ne ; 
y(k:k,1)  = xx ; 
xx        = xx + dx ; 
end ;
%                                  initialisation du champ de vitesse
uzero     = zeros(ns,1) ; 
%                                  initialisation du champ de pression
xx           = 0.5 * dx ; 
if (jinit == 0) 
pzero        = zeros(ne,1) ; 
else
for        k = 1:ne ; 
if (jinit == 1) 
pzero(k:k,1) = cos (pi*kmode*xx) ;  
end ; 
if (jinit == 2) 
pzero(k:k,1) = sin (0.5*pi*kmode*xx) ;  
end ; 
xx           = xx + dx ; 
end ;
end  ;
%
%
plot (x,uzero,'o3;vitesse initiale;',  y,pzero,'*k;pression initiale;' )
title('conditions initiales')
axis ([ 0  1.0   -1.1   1.1 ]) 
pause (1.5) 
%                                                   iterations en temps 
time           = 0 ; 
uu             = uzero ; 
pp             = pzero ; 
for       itps = 1 :  niter
   if (jclgau == 2) ; 
   xi        = zzgau * dt / dx ; 
   uu(1:1,1) = ( (1.d0 + xi)* uu(1:1,1)  - 2.d0 * dt * pp(1:1,1) / dx ) ...
               / (1.d0 - xi) ;
   end ;    
   for   j     = 2 :  ns-1 ; 
   uu(j:j,1)   = uu(j:j,1) - dt * (pp(j:j,1)-pp(j-1:j-1,1)) / dx ; 
   end ;
   if (jcldro == 2) ; 
   xi          = zzdro * dt / dx ; 
   uu(ns:ns,1) = ( (1.d0 - xi)* uu(ns:ns,1) + 2.d0 * dt* pp(ne:ne,1) / dx ) ...
              /    (1.d0 + xi) ; 
   end ;
   for   k     = 1 :  ne ;  
   pp(k:k,1)   = pp(k:k,1)  - dt * (uu(k+1:k+1,1) - uu(k:k,1)) / dx ; 
   end ;

plot (x,uu,'o3;vitesse;', y,pp,'*k;pression;' ) ; 
axis ([ 0  1.0   -1.1   1.1 ])
pause(.01) 
end ;

plot (x,uu,'o3;vitesse;', y,pp,'*k;pression;' ) ; 
axis ([ 0  1.0   -1.1   1.1 ])  

title('acoustique    condition limite d impedance')
print -dpdf acous-2-final.pdf
print -djpg acous-2-final.jpg



