% ******************************************************
% Modes
% ******************************************************
clear all

% *** set up Munk SSP ***

zs  = 1000.0;			% source depth
f   = 50; omega = 2 * pi * f	% frequency in Hertz
eps = 0.00737; c0  = 1500;

d   = 5000;			% bottom depth
nz  = 5000;			% number of finite-difference points
h   = d / nz; h2  = h * h;	% mesh spacing
z   = linspace( 0, d, nz );	% grid coordinates

x = 2 * ( z - 1300 ) / 1300;  
c = c0 * ( 1 + eps * ( x - 1 + exp( -x ) ) );

plot( z, c ); view( 90, 90 );
xlabel( 'Depth (m)' ); ylabel( 'Sound Speed (m/s)' )

% ******************************************************
% set up finite difference matrix
% ******************************************************

v = omega * omega ./ ( c .* c );

D = -2*ones( nz, 1 ) / h2 + v';
E =    ones( nz, 1 ) / h2;
A = spdiags([E D E], -1:1, nz, nz);

% ******************************************************
% Find eigenvalues and eigenvectors
% ******************************************************

m = 102 %50			% number of modes to use

k2 = eig( A );
k2 = flipud( sort( k2 ) );
k2 = k2( 1:m );
k = sqrt( k2 );

% inverse iteration for eigenfunctions

for ii = 1:m
   psi(:,ii) = ( A - k2( ii ) * speye( size( A ) ) ) \ ones( nz, 1 );
%  psi(:,ii) = ( A - k2( ii ) * speye( size( A ) ) ) \ psi(:,ii);	% extra iteration
   psi(:,ii) = psi(:,ii) / norm( psi(:,ii) );				% normalize
end

% ******************************************************
% plot the modes
% ******************************************************

pcolor( 1:m, z, psi ); shading interp; colormap(jet); colorbar; view( 0, -90 );
title('Mode shapes for the Munk profile')
xlabel( 'Mode' ); ylabel( 'Depth (m)' )

plot( z, psi( :, 1:25:m ) ); view( 90, 90 );
xlabel( 'Depth (m)' ); ylabel( 'Mode shape' )

% --- weigh modes by mode excitation

isd = zs / d * nz;		% index of source depth
a = psi( isd, : );
psi = psi * diag( a, 0 );	% scale modes by a

% --- show which modes are excited

pcolor( 1:m, z, psi ); shading interp; colormap(jet); colorbar;
title('Mode shapes for the Munk profile')
xlabel( 'Mode' ); ylabel( 'Depth (m)' )
view( 0, -90 );	% flip plot so that z-axis is pointing down

% ******************************************************
% form pressure field
% ******************************************************

nr = 201;
r = linspace( 1000.0, 100000.0, nr ); % nr ranges to 100 km
phase = diag( 1.0 ./ sqrt( k ) ) * exp( i * k * r ) * diag( 1.0 ./ sqrt( r ) );

p = psi * phase;
tl = 20.0 * log10( abs( p ) );

pcolor( r, z(1:10:nz), tl( 1:10:nz,:) ); caxis( [ -100 -60 ] ); shading interp; colormap( jet ); colorbar; view( 0, -90 );
xlabel( 'Range (m)' ); ylabel( 'Depth (m)' );
title( 'Normal mode intensity' )