clear all;
clc;
close all;
f=2000; 
w=2*pi*f;
c=344;
rho=1.15;
rhoc=rho*c;
k=w/c;
a=0.05;	% sphere radius
x=linspace(-1,1,25);
y=linspace(-1,1,25);

for m=1:25
   for l=1:25
      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)))/rhoc;
         pt(l,m)=pt(l,m)+exp(-i*k*x(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))/rhoc;
         Ix(l,m)=0.5*real(pt(l,m)*conj(ux(l,m)));
         Iy(l,m)=0.5*real(pt(l,m)*conj(uy(l,m)));
         amp=sqrt(Ix(l,m)^2+Iy(l,m)^2);
         Ix(l,m)=sqrt(sqrt(amp))*Ix(l,m)/amp;
         Iy(l,m)=sqrt(sqrt(amp))*Iy(l,m)/amp;
       else
          Ix(l,m)=0;
          Iy(l,m)=0;
       end
    end
end

figure(1)

H=quiver(x,y,Ix,Iy,1.5);
axis equal
xlabel('x');
ylabel('y');
title('Time-averaged Intensity')
axis([-1 1 -1 1])

% to store the images for web posting (bitmaps)
%filnme='spher';
%eval(['print -dbitmap ' filnme ' -f']);
