clear all;
clc;
close all;
Nsteps=30;
f=1825; 	%ka=10
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.3;	% piston radius
x=linspace(2/100,2,100);
y=linspace(2/100,2,100);

for m=1:100
   for l=1:100
      
      ux(l,m)=0;
      uy(l,m)=0;
      uxsum=0;
      uysum=0;
      for nrad=1:36
         sgma=nrad*a/36;
            for nang=1:36
               psi=nang*2*pi/36;
      			zs=sgma*cos(psi);
      			ys=sgma*sin(psi);
      			rp=sqrt((zs)^2+(y(l)-ys)^2+x(m)^2);
      			dp=exp(-i*k*rp)/rp*sgma;
               dur=dp*(1-i/(k*rp));
               uxsum=uxsum+dur*x(m)/rp;
               uysum=uysum+dur*(y(l)-ys)/rp;
   			end
		 end
       ux(l,m)=uxsum;
       uy(l,m)=uysum;
	end
end

for k1=1:2:100
   l1=1+(k1-1)/2;
   for k2=1:2:100
      l2=1+(k2-1)/2;
      amp=sqrt(ux(k1,k2)^2+uy(k1,k2)^2);
      ux3(l1,l2)=sqrt(amp)*ux(k1,k2)/amp;
      uy3(l1,l2)=sqrt(amp)*uy(k1,k2)/amp;
   end
   x2(l1)=x(k1);
   y2(l1)=x(k1);
end

figure(1)
ux2=real(ux3);
uy2=real(uy3);
H=quiver(x2,y2,ux2,uy2,1.5);
axis equal
axis([0 2 0 2])

%M=moviein(Nsteps); 

for kk=1:Nsteps
   ux2=real(ux3*exp(i*w*t(kk)));
   uy2=real(uy3*exp(i*w*t(kk)));
   H=quiver(x2,y2,ux2,uy2,1.5);
   axis equal
	xlabel('x');
   ylabel('y');
   title('Velocity')
   axis([0 2 0 2])
   
%   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)