%   acous4.m
%   Acoustique numerique  
%   Logiciel 'HaWAY4' 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.
%         haway4  12   fevrier  2001,    tp4     26  decembre 2001.
%         acous4  02   fevrier  2006, 25 fevrier 2025.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
c0 =  sprintf (' HaWAY4. 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
%     haway4     source ponctuelle sinusoidale
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   nombre de mailles             ne 
%   nombre de pas de temps        nt
%   nombre de sommets             ns
%   pas d espace                  dx
    ne    = 300
    niter = 400
ns    = ne + 1 ; 
dx    = 1. / ne ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   condition initiale  
%   jinit   parametre relatif a la condition initiale
%   kmode   numero du mode             
%   jinit = 1  champ nul
%         = 2  vitesse nulle,  pression en cosinus
%         = 3  vitesse nulle,  pression en sinus
jinit = 3 ; 
kmode = 5.5 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   condition limite a gauche 
%   jclgau  type de condition limite a gauche 
%   jclgau = 1  Dirichlet homogene ou non 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 = 1 ; 
ppgau  = 0 ; 
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  = 10 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   Elements instationnaires du champ de vitesse a gauche 
%   u(gauche) =  vvgau * sin (omegvv * t) ;
vvgau  = 0.3 ; 
omegvv = 35 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   source pour l equation de la masse du type 
%          flux + ffgg*sin(omega*time))   
%   situee a la position xsource
%   xsource    position de la source de masse  
%   flux       valeur du flux de la source de masse 
xsource  = .1 ; 
flux     = 0 ; 
ffgg     = .5 ; 
omega    = 35  ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   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 == 1) 
pzero        = zeros(ne,1) ; 
else
for        k = 1:ne ; 
if (jinit == 2) 
pzero(k:k,1) = cos (pi*kmode*xx) ;  
end ; 
if (jinit == 3) 
pzero(k:k,1) = sin (pi*(kmode+0.5)*xx) ;  
end ; 
xx           = xx + dx ; 
end ;
end  ;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                               conditions initiales
time           = 0 ;
timini         = time ;
uu             = uzero ; 
pp             = pzero ;

%
plot (x,uzero,'o3;vitesse initiale;',  y,pzero,'*k;pression initiale;' )
title('conditions initiales')
axis ([ 0  1.0   -1.8   1.8]) ; 
pause (1.5) 
%   


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                iterations en temps 
for       itps = 1 :  niter
   time        = time + dt ;
%                               condition de vitesse imposee  a gauche
   if (jclgau == 1) ;
   uu(1:1,1) =  vvgau * sin (omegvv * time) ;
   end ;    
%                                     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 + ffgg*sin(omega*time))  / dx  ;
   end ;
   end ;

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

plot (x,uu,'o3;vitesse;', y,pp,'*k;pression;' ) ;
axis ([ 0  1.0   -1.8   1.8]) ; 
title('acoustique     source ponctuelle sinusoidale')
print -dpdf acous-4-final.pdf
print -djpg acous-4-final.jpg


