% Animates quadrupole - last modified 10/29/99
clear all;
n=30;
M=moviein(n);		% for Matlab movie animation
f=500;				% frequency
w=2*pi*f;
t=linspace(pi/15/w,2*pi/w,30);	%30 time steps per period
c=344;
k=w/c;
a=0.04;			% sphere radius
x=linspace(-2,2,100);	% 4 m x 4 m area, with 10000 points
y=linspace(-2,2,100);

for m=1:100
   for l=1:100
      r(m,l)=sqrt((x(m)^2)+(y(l)^2));	%source at origin
      r2(m,l)=sqrt((x(m)^2)+((y(l)-0.1)^2)); %kd=0.91
      r3(m,l)=sqrt(((x(m)-0.1)^2)+(y(l)^2));
      r4(m,l)=sqrt(((x(m)-0.1)^2)+((y(l)-0.1)^2));
      if r(m,l)>=a
         pt(m,l)=exp(i*k*(r(m,l)-a))/r(m,l);
      else
         pt(m,l)=1/a;	% for visualization purposes
      end
      if r2(m,l)>=a
         pt2(m,l)=-exp(i*k*(r2(m,l)-a))/r2(m,l);
      else
         pt2(m,l)=-1/a;
      end
      if r3(m,l)>=a
         pt3(m,l)=-exp(i*k*(r3(m,l)-a))/r3(m,l);
      else
         pt3(m,l)=-1/a;
      end
		if r4(m,l)>=a
         pt4(m,l)=exp(i*k*(r4(m,l)-a))/r4(m,l);
      else
         pt4(m,l)=1/a;
      end
    end
end

pt=pt+pt2+pt3+pt4;
pmag=sqrt(sqrt(abs(pt)));	% 1/4 power of amplitude for better contrast
pphase=angle(pt);
clim=max(max(pmag));
V=[-clim clim];		% sets color scheme
pt=pmag.*exp(i*pphase);	% red: maximum positive
				% blue: minimum negative pressure
figure(1)
axis equal		% to eliminate distortion
axis([-2 2 -2 2])
for l=1:30
   pt2=real(pt*exp(-i*w*t(l)));
   h=pcolor(x,y,pt2);
   axis equal
	xlabel('x');
	ylabel('y');
	shading interp;
   caxis(V);
   axis([-2 2 -2 2])
   M(:,l)=getframe;
% to store the images for web posting (bitmaps)
%filnme=strcat('spher',int2str(l));
%eval(['print -dbitmap ' filnme ' -f']);
end