%   acous3.m
%   Acoustique numerique  
%   Logiciel 'HaWAY3' avec Octave-Matlab.
%   fd    haway1  1 fevrier 2001,        tp1     16 octobre 2001.
%         haway2  06   fevrier  2001,    tp2     20 novembre 2001.
%         haway3  10   fevrier  2001,    tp3     14, 26  decembre 2001.
%         acous3  02   fevrier  2006, 25 fevrier 2025.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
c0 =  sprintf (' HaWAY3. Resolution numerique de l equation des ondes ')
%     haway1     premier test avec le mode le plus bas
%                conditions de Dirichlet homogenes en vitesse
%     haway2     condition limite d impedance
%                condition initiale modale       
%     haway3     terme de source ponctuelle pour la masse
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   nombre de mailles             ne 
%   nombre de pas de temps        nt
%   nombre de sommets             ns
%   pas d espace                  dx
    ne    = 100
    niter = 200
ns    = ne + 1 ; 
dx    = 1. / ne ; 
%          xsource    position de la source de masse 
%          flux       flux de la source de masse 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   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 = 0 ; 
kmode = 0 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   condition limite a gauche 
%   jclgau  type de condition limite a gauche 
%   jclgau = 1  Dirichlet homogene en vitesse
%          = 2  Condition d impedance ;   pp = ppgau - zzgau  * uu
%   ppgau       valeur de la pression a gauche   
%   zzgau       valeur de l impedence a gauche   
jclgau = 2 ; 
ppgau  = 1 ; 
zzgau  = 1 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   condition limite a droite 
%   jcldro  type de condition limite a droite 
%   jcldro = 1  Dirichlet homogene en vitesse
%          = 2  Condition d impedance ;    pp = zzdro * uu
%   zzdro       valeur de l impedence a droite   
jcldro = 2 ; 
zzdro  = 2 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   source pour l equation de la masse
%   xsource    position de la source de masse  
%   flux       valeur du flux de la source de masse 
xsource  = .1 ; 
flux     = .4 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   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 ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                 pretraitement de la source de masse
      for   k = 1:ne ; 
      zz      = xsource / dx  ;
      if ((k <= zz) && (zz < k+1))  jsource = k ; 
      end ;
      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  ;
%
% gset yrange[-1.1:1.1] ; gset xrange[0: 1.0] ; 
axis ([ 0  1.0   -1.1   1.1 ]) ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                               conditions initiales
time           = 0 ;
timini         = time ;
uu             = uzero ; 
pp             = pzero ; 
[chn1,msg]  = sprintf (' nombre de  mailles : %4i',ne) ; 
%%%%%  if jdt==1 ;  [chcfl,msg] = sprintf ('            cfl : %5.2f',cfl) ; end ;
%%%%%  if jdt==2 ;  [chcfl,msg] = sprintf ('         dt : %9.2e',dt) ; end ;
%%%%%  [chtps,msg] = sprintf ('           temps   : %7.4f',timini); 
%%%%%  [chtps,msg] = sprintf ('           temps   : %7.4f',timini); 
%%%%%  titx = [chn1 , chcfl , chtps] ;    title (c0)  , xlabel (titx) ;
%%%%%  gset yrange[-1.1:1.1] ; gset xrange[0: 1.0] ; 
%
plot (x,uzero,'o3;vitesse initiale;',  y,pzero,'*k;pression initiale;' )
title('conditions initiales')
axis ([ 0  1.0   -.4   .9 ])  
pause (1.5) 
%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                iterations en temps 
for       itps = 1 :  niter
   time        = time + dt ;
%                                     condition d impedance a gauche
   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)-ppgau) / 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 ; 
   if (k == jsource) pp(k:k,1) =  pp(k:k,1) - dt * flux / dx ;   end ;
   end ;


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

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

title('acoustique     source ponctuelle pour la masse')
print -dpdf acous-3-final.pdf
print -djpg acous-3-final.jpg


