clear all;
clc;
close all;
Nsteps=30;
f=2000; 
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.05;	% sphere radius
x=linspace(-1,1,200);
y=linspace(-1,1,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(l,m)=exp(-i*k*(r(m,l)-a))/r(m,l);
         pt(l,m)=pt(l,m)+exp(-i*k*x(m));
       else
         pt(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);	

figure(1)
pt2=real(pt);
H=pcolor(x,y,pt2);
axis equal
axis([-1 1 -1 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([-1 1 -1 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)