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.


