function [cp, cg, evf] = amodes (Z, c, frq, nm, zr)
%
% function [cp, cg, evf] = amodes (Z, c, frq, nm, zr)
%
% Calculation of Acoustic Modes
%
% Calculates acoustic modes modes, P, from
% Pzz+ [w^2/c^2 - ev^2]*P=0, subject to P(-D)=0 and P(0)=0
% (I am at the moment a little unsure of the boundary conditions!)
% (And how to implement them in this matrix formalism.)
%
% Uses Numerov's method following the suggestion of Jensen, Kuperman,
% Porter, and Schmitt, Computational Ocean Acoustics, p. 282-286.
% Googlebooks for a preview:
% http://books.google.com/books?id=QHtx4zYPbzMC
%
% The approach here first calculates the components of the matrices required
% for Numerov formulation of the differential equation. The eigenvalue
% equation that results is: A*phi=(k^2)*B*phi, which is a standard eigenvalue
% equation. The eigenvectors phi are found with the call to eigs(A,B,...).
%
% Roughly, B is tridiagonal with 10/12 on the main diagonal and 1/12 on
% the +/-1 diagonals.
%
% INPUT:
%
% Z,c are input depth and sound speed. These must already be interpolated
% to an evenly-spaced grid to achieve a desired accuracy. 5-m spacing seems
% to work o.k.; this depends on frequency, of course.
%
% frq and nm are frequency and number of modes desired.
%
% zr is a vector of depths at which the mode shapes are desired.
% This might correspond to the depths of hydrophones in a simulated
% reception. amodes interpolates the modes it finds to the zr depths.
%
% Some timing on a 2.4 GHz Q6600 cpu:
% At 100 Hz, using
% 1000 sound speeds, 50 modes in 4.3 s; d(cp)=6e-4, d(cg)=1e-3
% 2500 sound speeds, 50 modes in 38.2 s; d(cp)=2e-5, d(cg)=1e-3
% 5000 sound speeds, 50 modes in 259.0 s; d(cp)=2e-5, d(cg)=1e-3
% d(cp) is difference in phase speed relative to chebmode,
% Same for d(cg) and group speed. These taken to be some measure
% of the accuracy of the calculations.
%
% OUTPUT:
%
% cp and cg are mode phase and group speed.
% evf is a matrix of nm modes at the depths zr.
%
% Uses eigs to calculate only nm modes. The use of eig is possible, but
% this calculates length(Z) modes, most of which were tossed out!
%
% B. Dushaw February-April 2009
% allocate for evf
evf=zeros(length(zr),nm,length(frq));
% Grid spacing
h=abs(Z(2)-Z(1));
% Frequency to radians/s
w=2*pi*frq;
% Calculate slowness.
sp=1./c;
% Main loop over frequency
for kkk=1:length(frq),
sp2=(w(kkk)^2)*sp.*sp; % frequency over sound speed squared
% Diagonal elements
dj=(-2/(h^2) + (10/12)*sp2)/(h*1024);
% Upper off diagonal
eja=(1/(h^2) + (1/12)*sp2(2:end))/(h*1024);
% Lower off diagonal
ejb=(1/(h^2) + (1/12)*sp2(1:(end-1)))/(h*1024);
% Construct the matrix A - A is not symmetric.
A1=diag(dj);
A2=diag(eja,1);
A3=diag(ejb,-1);
A=A1+A2+A3;
A(1,:)=[]; % Unsure about the boundary conditions
A(:,1)=[]; % In their simpler example Jensen et al say to null out the first row and column.
% Seems to work better with this as well - better accuracy, and modes are zero
% at the surface and bottom like they are supposed to be.
A(end,:)=[];
A(:,end)=[];
clear A1 A2 A3
% Construct the matrix B - B is symmetric.
B1=(10/12)*eye(length(A));
B2=(1/12)*diag(ones(1,length(A)-1),1);
B3=(1/12)*diag(ones(1,length(A)-1),-1);
B=(B1+B2+B3)./(h*1024);
clear B1 B2 B3
if 1==1,
% Use eigs to just get the nm largest real eigenvalues
% Cuts computation time by a factor of 3-4 over eig.
% When only a small number of modes are requested, one needs to have
% more Lanczos vectors for convergence.
P=max(2*nm,75);
P=min(P,length(A)); % P can be at most length(A)
% tol is by default "eps" which is 2.2e-16 on my computer. Increasing this
% tolerance makes things run much faster, without noticeably affecting the
% accuracy of the modes or eigenvalues that I can tell.
tol=1e-10;
%tol=eps
% Options for eigs:
% 'tol'= numerical tolerance,
% 'maxit'=maximum number of iterations,
% 'p'=number of Lanczos vectors,
% 'disp'=verbosity of calculation.
% (there are other options)
opts=struct('tol',tol,'maxit',500,'p',P,'disp',0);
% All of the compute time is in this call.
[V,kr,flag]=eigs(sparse(A),sparse(B),nm,'LR',opts);
% Convergence=flag % If convergence is not zero, you have problems.
kr=diag(kr)'; % eigenvalues - KH in the literature.
% Arrange evalues and evectors in descending order
% These may be in order already.
[kr,I]=sort(kr,'descend');
V=V(:,I);
else
% Use eig to get all eigenvalues.
% Solves for ALL modes; very computationally expensive.
[V,kr]=eig(A,B);
kr=diag(kr)'; % eigenvalues - KH in the literature.
[kr,I]=sort(diag(kr),'descend');
V=V(:,I);
kr=kr(1:nm);
V=V(:,1:nm);
end
% kr are the wavenumbers
kr=sqrt(kr);
[n m]=size(V);
V=[zeros(1,m);V;zeros(1,m)]; % mode is zero at the surface and bottom.
% I don't understand how this works, but it
% seems to.
if 1==2,
% calculate local vertical wavenumber
for kk=1:nm,
KV=sqrt(w(kkk)*w(kkk)*sp.*sp - kr(kk)*kr(kk));
LV=2*pi./KV;
plot(LV/5,Z)
title(['Mode Number: ' num2str(kk)])
grid on
axis([0 35 -5000 0])
pause
end
end
if 1==2,
% Now the test...how well is the eigenvalue problem solved?
V1=V; V1(1,:)=[];
for kk=1:4,
Q=A*V1-kr(kk)*kr(kk)*B*V1;
plot(Q(:,kk),Z(2:end))
hold on
end
pause
% In my tests, Q came out to be at most 2e-16
end
% Interpolate modes to 5-m spacing for better integration accuracy.
s=sign(Z(end));
hs=5*s; % 5-m intervals seems to work fine; 10-m works ok as well, looks like.
ZP=Z(1):hs:Z(end);
hs=abs(hs);
% Normalize the modes
F=interp1(Z,V,ZP,'spline'); % always spline, never cubic!
F=F.*F;
N=hs*trapz(F);
N=sqrt(N);
N=ones(size(Z))*N;
V=V./N;
% Set the phase of modes at the top
if 1==2,
V=V.*(ones(length(V),1)*sign(V(3,:)));
else
for kook=1:nm,
I=find(abs(V(:,kook))>1e-3);
I=I(1);
V(:,kook)=V(:,kook)*sign(V(I,kook));
end
end
% Phase speed
cp(:,kkk)=w(kkk)./kr;
% Group speed per Jensen et al, or Munk et al 1995 ....
% V's are normalized just above, so we're good to go here.
Q=sp*ones(1,nm);
Q=V.*Q;
Q=interp1(Z,Q,ZP,'spline'); % always spline, never cubic!
Q=Q.*Q;
icg=hs*(w(kkk)./kr).*trapz(Q);
cg(:,kkk)=1./icg;
% Interpolate the modes to the requested instrument depths.
evf(:,:,kkk)=interp1(Z,V,zr,'spline');
end % end loop for frequencies.