% ffp pekeris wave guide
% Author: Peter Gerstoft, Jie Li
% Date: December 9, 2016

clear all;
% close all;

%% parameters

f=150;                   % frequency
omega=2*pi*f;            % angular frequency
zs=35;                   % source depth      
zr=25;                   % receiver depth
cw=1500;                 % sound speed of water
Nlayer=30;               % number of layers
dz=5;                    % layer thickness
iSint=(zs/dz)+1;         % source interface
IRint=(zr/dz)+1;         % receiver interface

for il=1:Nlayer-1   
    c(il)=cw;  z(il)=(il-1)*dz;  rho(il)=1;
end

c(Nlayer)=1800;          % bottom sound speed
rho(Nlayer)=1.8;         % density
z(Nlayer)=Nlayer*dz;     % depth of last layer

thick=diff(z);           % thickness of each layer
kp=omega./c;       
mu=rho*omega^2;
sw=-4*pi/(mu(iSint));    % source strength

kmax=omega/real(cw)*2; % maximum of kr  
Nfft=2048;               % sampling points of kr
dkr=kmax/(Nfft-1);       % sampling interval of kr
kr=(0:dkr:kmax);         
dr=2*pi/kmax; 
R=(0:(Nfft-1))*dr;       % range
R(1)=R(1)+10^(-10);      % to avoid divide by zero
eps = 3/(log10(exp(1))*max(R));              % contour offset, 60dB
kr=kr-1i*eps;            

%% wavenumber integration 

for ik=1:Nfft
    
    alpha=sqrt(kp.^2-kr(ik).^2)*1i;      % 1i*k_z
    vv=[ -sw/(2*pi)  0].';                   % vector v at the source interface; (Eq. 4.81)
    v=[1 -mu(Nlayer)/alpha(Nlayer)].';       % vector v at the lowest interface (Eq. 4.185 and Eq. 4.84);
    
    % bottom to receiver layer
    for il=(Nlayer-1):-1:1
        cmat=[-alpha(il) alpha(il);  -mu(il)   -mu(il)];
        P=cmat*diag([exp(alpha(il)*thick(il)) exp(-alpha(il)*thick(il))])/cmat;  % (Eq. 4.78 and Eq. 4.184)
        v=P*v;                               % Eq. 4.77
        vr(:,il)=v;                          % save for the receiver
        if (il<iSint)                        % source
            vv=P*vv;      vs(:,il)=vv;
        else
            vs(:,il)=[0 0]';
        end
    end                               % done layer loop
    
    w=-vv(2)/v(2);                    % pressure-release boundary condtion, see Eq. 4.186
    phif(ik,:)=vr(2,:)*w+vs(2,:);     % Eq. 4.82 add inhomogeneous term
    
end                                   % wavenumber loop

phi=phif(:,IRint).';

for irr=IRint
    phi=phif(:,irr).';

    for nr = 1:length(R)
        p3(nr) = sum(phi.*besselj(0,kr*R(nr)).*kr)*dkr;
    end
    
    p3=p3.*exp(eps*R);
    fieldh(irr,:)=-20*log10(abs(p3));
    
end

%% plot

figure
plot(real(kr), abs(phi), 'r');hold on;
xlabel('Horizontal wavenumber'); 
ylabel('Amplitude');                 

figure;
plot(R(1:200)/10^3,fieldh(IRint,1:200),'r');hold on;
xlabel('Range (km)'); 
ylabel('Transmission loss (dB)');   
set(gca,'YDir','reverse');
