% Sound field calculation, the sea bottom is described by two parameters: P and Q
% Version 2     2012
% Main function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The output of the main function:
% complex horizontal wavenumber of mode N:        kr(N)+i*ki(N)
% eigen function of mode N in water layer:        Pm(N,:) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all; close all; clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Input parameters by user
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%frequency (Hz)
f=200;            
%% Parameters of the bottom
P=11;  Q=0.4;     %parameters of complex reflection coefficients
%% Parameters in water volume
H=100;                               %depth of water(m)
den0=1.0;                            %density of water
depth(1)=0; ssp_in(1)=1490.0;        %sound speed profile. Format: depth;sound speed.
depth(2)=100.0;ssp_in(2)=1490.0;

N=1000;          %numbers of layers 
 
%% Parameters for calculate sound field
sd=95.0;            %depth of source
rd=95.0;            %depth of receiver, only used for TL 
r=500:10:20000;    %range: start:step:stop(m);  for the calculation of sound intensity

%% limit the range of phase speed for searching  
cmin=1400;
cmax=2000;

alpha=pi/P;      %critical angle, modified 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W=2*pi*f;
h=1/N; 
%% SSP interpolation
dep=0:H/N:H;
ssp=interp1(depth,ssp_in,dep);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cmax=min(cmax,ssp(N+1)/cos(alpha));   
cmin=min(ssp)+1.0e-5;
%% shooting %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jj=1; 
dc=0.1;   % searching step of phase speed
for cn=cmin:dc:cmax                                 
    kn=W/cn;                                 
    Z0=cal_z0(W,P,kn,ssp(N+1),H);
    Z1=cal_z1(Z0,W,P,kn,ssp(N+1),ssp(N),H,h);       
    Z=cal_z(Z0,Z1,W,kn,ssp,H,N);                    
    Zsurf(jj)=real(Z(N+1));                           
    jj=jj+1;
end

%% find the cx
x0=Zsurf(1);                                                                         
kk=0;
for  ii=2:jj-1;               % modified 
    if (Zsurf(ii)*x0<=0)    
        kk=kk+1;
        if abs(Zsurf(ii))>abs(Zsurf(ii-1))  % finding the accurate solution
            Zsurf_ok(kk)=Zsurf(ii-1);
        else
            Zsurf_ok(kk)=Zsurf(ii);
        end
    XI(kk)=(ii-1)*dc;
        x0=Zsurf(ii);      
    end
end

figure,
plot(cmin:dc:cmax,Zsurf,'*-',cmin+XI,Zsurf_ok,'rs');
xlabel('modal phase velocity(m/s)')
ylabel('value at surface')
title('Shooting');

Nmode=kk;  % Nmode is the max mode number
Pm=zeros(Nmode,N+1);     

for ii=1:Nmode
    cx(ii)=Search_c(W,P,H,ssp,N,h,cmin+XI(ii)-dc,cmin+XI(ii)+dc); 
    %%%%%%%%%%%%
    kr(ii)=W/cx(ii);
    Z0=cal_z0(W,P,kr(ii),ssp(N+1),H);
    Z1=cal_z1(Z0,W,P,kr(ii),ssp(N+1),ssp(N),H,h);
    Pm(ii,N+1:-1:1)=cal_z(Z0,Z1,W,kr(ii),ssp,H,N);
    %%%%%%%%%%%%%
    Sn=Pm(ii,:)*Pm(ii,:)'*h;
    kh=W/ssp(1);
    x=sqrt(1-(kr(ii)/kh).^2);
    ki(ii)=Q*x/(kr(ii)*H*H*Sn+P/kh/x);
    %%%%%%%%%%%%%
    disp(sprintf('wavenumber of mode %d (real part and imaginary part): %e,%e',ii,kr(ii),ki(ii)))  
    Nn(ii)=cal_Nn(W,Pm(ii,:),den0,kr(ii),ssp,H/N,P);          
    Pm(ii,:)=Pm(ii,:)*sqrt(Nn(ii));                                  %mode depth functions
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%output   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%   To display all the eigen functions
for n=1:Nmode
dep1=0:2:H;    
fy(n,:)=interp1(dep,real(Pm(n,:)),dep1);;
end

figure,
subplot(1,n,1);
plot(fy(1,:)',dep1);
set(gca,'ydir','reverse');
title('mode1');
ylabel('depth(m)');

for m=2:n
subplot(1,n,m);
plot(fy(m,:)',dep1);

set(gca,'ydir','reverse');
title(['mode' int2str(m)]);
end

%% the TL 
Nrange=length(r);


for ii=1:Nmode
    Pm_z(ii)=interp1(dep,Pm(ii,:),rd);
    Pm_z0(ii)=interp1(dep,Pm(ii,:),sd);
end

for ii=1:Nrange   
    for nn=1:Nmode
    K(nn,ii)=exp(i*kr(nn)*r(ii)-ki(nn).*r(ii)+i*(pi/4))./sqrt(kr(nn));
    P(nn,ii)=sqrt(1/(8*pi*r(ii))).*Pm_z(nn).*Pm_z0(nn).*K(nn,ii);
    end
    p(ii)=sum(P(:,ii));      % sound pressure
    A(ii)=angle(p(ii));      % phase of sound pressure
    S(ii)=abs(p(ii));        % amplitude of sound pressure
   
end

I=-20*log10(S);      % Transmission Lose
figure
plot(r,I);
set(gca,'ydir','reverse');
xlabel('range(m)');
ylabel('Transmission Lose(dB)');


