% Program to test the accuracy of amodes against chebmode.
% Uses cssprofile by default; edit to change.
clear all
figure(2)
clf
figure(1)
clf

path(path,'/home/dushaw/matlab/aogmode')
% We need the cssprofile program and we are comparing to chebmode.

frq=28;   % frequency
nm=50;    % number of modes

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 1==2
% Get sound speed from eigenray/mpiram formatted data file.
   C=load('test.ssp');
% temp.ssp is a file of z c at various ranges.
   I=C(:,1)==-1.;
   I=find(I);
   N1=I(2)-1;
   z=C(2:N1,1);
   c=C(:,2);
   N2=length(c)/N1;
   c=reshape(c,N1,N2);
   r=c(1,:);
   c(1,:)=[];

% Use just first profile for now.
   z=-z;
   c=c(:,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Resample c onto an evenly-spaced grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Number of samples in depth for the grid.
% This can be adjusted for more accuracy/higher number of modes.
% Looks like NS=1000 over 5000 m depth, or 5-m spacing seems to 
% work about right.
  NS=5500;   

% Depth grid
  h=(max(z)-min(z))/NS;
  Z=(0:-h:min(z))';

  c=interp1(z,c,Z,'cubic');
  z=Z;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

z=(0:10:5000)';
c=cssprofile(z);
z=-z;

zr=-(0:10:5000)';

tic
[CP, CG, EVF] = amodes (z, c, frq, nm, zr);
toc

MN=34:36;
MN=1:5;
MN=25:35;
subplot(121)
plot(EVF(:,MN),zr/1000)
ylabel('Depth (km)')
hold on

tic
[cp, cg, evf]=chebmode(-z,c,frq,nm,-zr);
toc

plot(evf(:,MN),zr/1000,'*')
subplot(122)
DM=EVF-evf;
plot(DM(:,MN),zr/1000)
ylabel('Depth (km)')
xlabel('Mode Difference')

figure(2)
subplot(211)
plot(CP-cp','*')
hold on
plot(CP-cp','-')
title('Phase Speed Difference:  amode vs. chebmode')
ylabel('Speed (m/s)')

subplot(212)
plot(CG-cg','*')
hold on
plot(CG-cg','-')
title('Group Speed Difference:  amode vs. chebmode')
ylabel('Speed (m/s)')
xlabel('Mode Number')