% Duct with enlargement - last modified 11/03/99
clear all;
clc;
close all;
Nsteps=30;
f=86*2;	% frequency
rho=1.15;
c=344.;
rhoc=rho*c;
l=1;	%length of duct
a1=0.0762/2;
S1=pi*a1^2;
a2=0.1524/2;
S2=pi*a2^2;
m=S2/S1;
w=2*pi*f;
k=w/c;
t=linspace(2*pi/Nsteps/w,2*pi/w,Nsteps);	%Nsteps time steps per period
alpha1=2.93e-5/a1.*sqrt(f);
alpha2=2.93e-5/a2.*sqrt(f);
khat=k-i*alpha1;
kl=khat*l;
U0=0.01;
R=(1-m)/(m+1);
A=rho*c*U0./(exp(i*kl)-R*exp(-i*kl));
B=R*A;
C=A+B;

x=linspace(0,2,100);	% computational domain
y=linspace(-0.2,0.2,32);

for m=1:50
   for k=8:24
      pt(k,m)=A*exp(i*khat*(l-x(m)))-B*exp(-i*khat*(l-x(m)));
   end
end

for m=51:100
   for k=1:32
      pt(k,m)=C*exp(i*khat*(l-x(m)));
   end
end

pmag=abs(pt);				
pphase=angle(pt);
clim=max(max(pmag));
V=[-clim clim];	% sets color scheme
pt=pmag.*exp(i*pphase);	% red: maximum positive
								% blue: minimum negative pressure
figure(1)
pt2=real(pt);
H=pcolor(x,y,pt2);
axis equal % to eliminate distortion
axis([0 2 -.2 .2])

M=moviein(Nsteps);	% for Matlab movie animation

for k=1:Nsteps
   pt2=real(pt*exp(i*w*t(k)));
   H=pcolor(x,y,pt2);
   axis equal
	xlabel('x');
   ylabel('y');
   title('pressure')
	shading interp;
   caxis(V);
   axis([0 2 -.2 .2])
   M(:,k)=getframe;
% to store the images for web posting (bitmaps)
%   filnme=strcat('spher',int2str(k));
%   eval(['print -dbitmap ' filnme ' -f']);
end

movie(M,20)
