%Angular Spectrum Propagation Code
%
%Source condition file should set the following variables
% Usxy  - two-d array with amplitude and phase of source velocity
% k - wavenumber of sound
% alpha - attenuation of sound
% xvec - vector with x in space
% lenx - length of x vector
% yvec - vector with y in space
% leny - length of y vec

%Written by Robin Cleveland, Boston University
%ME720, Acoustics 2
%This version March 2008

%%%%%%%%%%%
%USER INPUT
%%%%%%%%%%%

%Select z locations were output is desired
zvec=[0.01:0.01:1]*zrayl;

%Set to 1 to create plots of Phis and Psxy at each z
PLOT=0;   

%Selection locations for output in x and y planes
xidx=round(lenx/2+1);yidx=round(leny/2+1);  %Midpoints of data sets

%%%%%%%%%%%%%%%%%%%
%CALCULATIONS START
%%%%%%%%%%%%%%%%%%%

lenz=length(zvec);
%Create empty arrays for figures of axial pressure, xz, and yz
paxial=0*zvec;
pxz=zeros(lenx,lenz);
pyz=zeros(leny,lenz);
subplotlenx=ceil(sqrt(lenz));
subplotleny=ceil(lenz/subplotlenx);

%Create u,v spatial transform vectors
dx=(xvec(end)-xvec(1))/(lenx-1);
dy=(yvec(end)-yvec(1))/(leny-1);
u=[-lenx/2:lenx/2-1]'/(dx*lenx);
v=[-leny/2:leny/2-1]/(dy*leny);
u=fftshift(u); %reorder to match matlab FFT format
v=fftshift(v); %reorder to match matlab FFT format

kz=sqrt((k+i*alpha)^2-(2*pi)^2*((u.*u)*ones(size(v))+ones(size(u))*(v.*v)));
    
%Transform source condition to k-space
Usuv=ifft2(Usxy);
Phis=Usuv./(i*kz);  %Convert to velocity potential source

figure(2);
imagesc(fftshift(abs(Usuv)));
title ('Source in transform-space');

%Step through all z output locations
for zcnt=1:lenz,
	z=zvec(zcnt);
    
%Calculate P in k-space at z
    Phiz(:,:)=Phis(:,:).*exp(i*kz*z);   
%Transform back to spatial domain
    Pzxy(:,:)=fft2(i*k*Phiz(:,:));
	
    if PLOT>0,
    %plot k-space data
      figure(3);
      subplot(subplotlenx,subplotleny,zcnt);	
      imagesc(fftshift(abs(Phiz(:,:))))
      axis('equal');axis('tight');
      title(['Angular Spectrum  z=',num2str(z)]);
    %plot xy-data
      figure(4);
      subplot(subplotlenx,subplotleny,zcnt);	
	  imagesc(xvec*1e3,yvec*1e3,abs(Pzxy(:,:)));
	  axis('equal');axis('tight');
      title(['Amplitude z=',num2str(z)]);
	  ylabel('y (mm)');
      xlabel('x (mm)');
    end
    
%Save data for other plots
     paxial(zcnt)=Pzxy(xidx,yidx);
     pxz(:,zcnt)=Pzxy(:,yidx);
     pyz(:,zcnt)=Pzxy(xidx,:)';
	
end   %for z loop

%Create image of x-z data
figure(5)
imagesc(zvec*1e3,xvec*1e3,abs(pxz));
axis([0 zvec(end) -1.5*a 1.5*a]*1e3);
title(['y =' num2str(1000*yvec(yidx)) ' mm']);
ylabel('x (mm)');xlabel('z (mm)');

%Create image of y-z data
figure(6)
imagesc(zvec*1e3,yvec*1e3,abs(pyz));
axis([0 zvec(end) -1.5*a 1.5*a]*1e3);
title(['x =' num2str(1000*xvec(xidx)) ' mm']);
ylabel('y (mm)');xlabel('z (mm)');

%Plot axial pressure amplitude
figure(7);
plot(zvec*1e3,abs(paxial));
title('Axial pressure');
xlabel('z (mm)');ylabel('|p|/p0');