clear all;
clc;
close all;
Nsteps=30;
f=30; 
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.2;	% sphere radius
x=linspace(0,10,200);
y=linspace(0,10,200);

for m=1:200
   for l=1:200
      r(m,l)=sqrt((x(m)^2)+(y(l)^2));
      if r(m,l)>=a
         pt(m,l)=exp(i*k*(r(m,l)-a))/r(m,l);
         ut(m,l)=pt(m,l)*(1-j/(k*r(m,l)));
      else
         pt(m,l)=0;
         ut(m,l)=0;
      end
    end
end

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

umag=sqrt(sqrt(abs(ut)));
pphase=angle(ut);
clim2=max(max(umag));
if clim2>clim
   V=[-clim2 clim2];
end
ut=umag.*exp(i*pphase);

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

subplot(1,2,2)
ut2=real(ut);
H=pcolor(x,y,ut2);
axis equal
axis([0 10 0 10])

M=moviein(Nsteps); 

for kk=1:Nsteps
   pt2=real(pt*exp(-i*w*t(kk)));
   subplot(1,2,1)
   h=pcolor(x,y,pt2);
   axis equal
	xlabel('x');
   ylabel('y');
   title('Pressure')
   shading interp;
   caxis(V);
   axis([0 10 0 10])
   
   ut2=real(ut*exp(-i*w*t(kk)));
   subplot(1,2,2)
   h=pcolor(x,y,ut2);
   axis equal
	xlabel('x');
   ylabel('y');
   title('Velocity')
	shading interp;
   caxis(V);
   axis([0 10 0 10])

   M(:,kk)=getframe;
% to store the images for web posting (bitmaps)
filnme=strcat('spher',int2str(kk));
eval(['print -dbitmap ' filnme ' -f']);
end
