%   acous5.m
%   Acoustique numerique  
%   Logiciel 'HaWAY5' 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.
%         haway5  26   fevrier  2001,    tp5     26  decembre 2001.
%         acous5  02   fevrier  2006, 25 fevrier 2025.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
c0 =  sprintf (' HaWAY5. 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
%     haway5     source ponctuelle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   nombre de mailles             ne 
%   nombre de pas de temps        nt
    ne    = 100
    niter = 150
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   calcul du nombre de sommets   ns       et du   pas d espace    dx
  ns    = ne + 1  ;
dx    = 1. / ne ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   condition initiale de pression  jinitpr 
%   jinitpr  = 0   champ de pression nul
%            = 1   pression en cosinus
%            = 2   pression en sinus
%   zmodepr : numero du mode en pression
jinitpr  = 0 ; 
zmodepr  = 5 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   condition initiale de vitesse   jinitvi
%   jinitvi  = 0    champ de vitesse  nul
%            = 1    vitesse  en cosinus
%            = 2    vitesse  en sinus
%            = 3    corde pincee : vitesse  affine par morceaux
jinitvi  = 0 ; 
zmodevi  = 0.1 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   condition limite a gauche 
%   type de condition limite a gauche :  jclgau       1
%   jclgau = 1  Dirichlet en vitesse ;   uu(x=0) = vvgau * sin(omegvv * t)
%          = 2  Condition d impedance ;  pp(x=0) = ppgau - zzgau  * uu(x=0)
%   vvgau       vitesse de reference  a gauche 
%   omegvv      pulsation de reference a gauche
%   ppgau       valeur de la pression a gauche   
%   zzgau       valeur de l impedence a gauche   
jclgau = 1 ; 
vvgau  = 0 ; 
omegvv = 35 ;  
ppgau  = 0 ; 
zzgau  = 1 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   condition limite a droite 
%   jcldro  type de condition limite a droite 
%   jcldro = 1  Dirichlet homogene en vitesse :     uu(x=1) = 0
%          = 2  Condition d impedance a droite ;    pp(x=1) = zzdro * uu(x=1)
%   zzdro       valeur de l impedence a droite   
jcldro = 2 ; 
zzdro  = 0 ;  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   source pour l equation de la masse 
%   source ponctuelle pour l equation de la masse de la forme suivante :   
%   q=delta(x-xsource)*[flux + ffcc * cos(omega*t) + ffsi * sin(omega*t)]
%   xsource    position de la source de masse  
%   flux       valeur du flux constant de la source de masse 
%   ffcc       valeur du flux en cosinus de la source de masse 
%   ffsi       valeur du flux en sinus de la source de masse 
%   omega      pulsation de la source de masse
xsource = .292 ; 
flux    = 0 ; 
ffco    = 0 ; 
ffsi    = 0.5 ; 
omega   = 16.4995063 ;  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   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.4 ; 
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
      zz      = xsource / dx  ;
      for   k = 1:ne ;  
      if ((k <= zz) && (zz < k+1))  jsource = k ;   
      end ;   end ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          pretraitement du pincement du champ de vitesse
      if (jinitvi == 3) ; 
      xxvi    = zmodevi 
      xx      = 0 
      for   k = 1:ne 
      if ((xx <= xxvi) & (xxvi < xx+dx))  jxxvi = k ,  end ;
      xx      = xx + dx ; end ;  end ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%                                  initialisation du champ de vitesse  
if (jinitvi == 0) ;     uzero        = zeros(ns,1) ; else ; 
   xx      = 0 ;   for   j = 1:ns ; 
   if (jinitvi == 1) ;  uzero(k:k,1) = cos (pi*zmodevi*xx) ; end ; 
   if (jinitvi == 2) ;  uzero(k:k,1) = sin (pi*zmodevi*xx) ; end ; 
   if (jinitvi == 3) ;  if (j <= jxxvi) ; 
                        uzero(j:j,1) = (j-1) / ( jxxvi-1) ; else 
                        uzero(j:j,1) = (ns-j) / ( ns-jxxvi) ; end ;  end ; 
   xx = xx + dx ;  end ;  end ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%                                  initialisation du champ de pression
xx           = 0.5 * dx ; 
if (jinitpr == 0) ; pzero        = zeros(ne,1) ; 
else
for        k = 1:ne ; 
if (jinitpr == 1) ; pzero(k:k,1) = cos (pi*zmodepr*xx) ;   end ; 
if (jinitpr == 2) ; pzero(k:k,1) = sin (pi*zmodepr*xx) ;   end ; 
xx           = xx + dx ; 
end ;
end  ;
%
% gset yrange[-1.1:1.1] ; gset xrange[0: 1.0] ; 
axis ([ 0  1.0   -1.75   1.75 ]) ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                               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   -.55   .55]) ;
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 +  ffco * cos(omega*time) ...
                              +  ffsi * sin(omega*time))  / dx  ;  end ;
   end ;



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

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






