%K&F Figure 8.22
clear all;
clc;
close all;
Nsteps=30;
f=4380; % kd=8
w=2*pi*f;
t=linspace(2*pi/Nsteps/w,2*pi/w,Nsteps);
c=344;
rho=1.15;
rhoc=rho*c;
k=w/c;
a=0.01;	% piston radius
x=linspace(2/200,1,100);
y=linspace(2/200,1,100);

for m=1:100
   for l=1:100
      r=sqrt((x(m)^2)+(y(l)^2));
      r2=sqrt((y(l)^2)+((x(m)-.1)^2));
      r3=sqrt((y(l)^2)+((x(m)+.1)^2));
      r4=sqrt((y(l)^2)+((x(m)-.2)^2));
      r5=sqrt((y(l)^2)+((x(m)+.2)^2));
         p1=exp(-i*k*(r-a))/r;
         p2=exp(-i*k*(r2-a))/r2;
         p3=exp(-i*k*(r3-a))/r3;
         p4=exp(-i*k*(r4-a))/r4;
         p5=exp(-i*k*(r5-a))/r5;
         pt(l,m)=p1+p2+p3+p4+p5;
     end
end

pmag=sqrt(sqrt(abs(pt)));
pphase=angle(pt);
clim=max(max(pmag));
V=[-clim clim];							
pt=pmag.*exp(i*pphase);	

figure(1)
pt2=real(pt);
H=pcolor(x,y,pt2);
axis equal
axis([0 1 0 1])

%M=moviein(Nsteps); 

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

%movie(M,20)