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,20);
y=linspace(-1,1,20);

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

figure(1)
ux2=real(ux);
uy2=real(uy);
H=quiver(x,y,ux2,uy2,1.5);
axis equal
axis([-1 1 -1 1])

M=moviein(Nsteps); 

for kk=1:Nsteps
   ux2=real(ux*exp(i*w*t(kk)));
   uy2=real(uy*exp(i*w*t(kk)));
   H=quiver(x,y,ux2,uy2,1.5);
   axis equal
	xlabel('x');
   ylabel('y');
   title('Velocity')
   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)