% This is to calculate the complex pressure only
% To animate, use same m-file as, say, hw10

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
      
      pt(l,m)=0;
      psum=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;
      			psum=psum+dp;
   			end
		 end
		pt(l,m)=psum;
	end
m
end

pmag=sqrt(abs(pt));
pphase=angle(pt);
clim=max(max(pmag));
V=[-clim clim];							
pt=pmag.*exp(i*pphase);	
