function [PHI C PHI2]=MODES(dz,N2,Nm,flag_0)
% USAGE: [PHI C PHI2]=MODES(dz,N2,Nm)
% Obtain vertical modes for arbitrary stratification
%
% INPUTS:
% dz  [1 x 1 ]   vertical spacing (must be uniform)
% N2  [Nz x 1]   stratification 
% Nm  [1 x 1 ]   number of modes to extract 
% flag_0 [1 x 1] Set to 1 (default) to include mode-0, set to 0 to exclude mode-0
%
% OUTPUTS:
% PHI  [Nm x Nz]  pressure and velocity structure eigenfunctions 
% C    [Nm x 1 ]  eigenspeeds 
% PHI2 [Nm x Nz]  vertical velocity structure eigenfunctions 
%
% Sam Kelly, 21 AUG 2015 (smkelly@d.umn.edu)

% Default number of modes
if nargin<4
	flag_0=1;
	if nargin<3
		Nm=length(N2)-1;
	end
end

if flag_0==1
	% Subtract barotropic mode from modes to solve
	Nm=Nm-1;
end

N=length(N2);
H=dz*N;

% Second derivative matrix
D2=zeros(N-1,N-1);
for i=2:N-2
    D2(i,i-1)=1/dz^2;
    D2(i,i)=-2/dz^2;
    D2(i,i+1)=1/dz^2;
end

% Rigid lid
D2(1,1)=-2/dz^2; 
D2(1,2)=1/dz^2;

% Flat bottom
D2(N-1,N-1)=-2/dz^2;
D2(N-1,N-2)=1/dz^2;

% "A" Matrix
N2_tmp=(N2(1:end-1)+N2(2:end))/2;
A=diag(-N2_tmp);% Hydrostatic

% Solve generalized eigenvalue problem
opts.issym=1;
opts.isreal=1;
[PHI k2]=eigs(D2,A,Nm,'SM',opts);

% Sort modes by eigenspeed
k2=diag(k2);
k2(k2<0)=Inf;
C=1./sqrt(k2); % eigenspeed
[C,ind]=sort(C,1,'descend');
PHI=PHI(:,ind);


% Add in boundary values
PHI=[zeros([1 Nm]); PHI; zeros([1 Nm])]; 

% Average between points to get W structure
PHI2=(PHI(1:end-1,:)+PHI(2:end,:))/2;
 
% Take derivative to get U and P structure
%PHI=-diff(PHI)./dz;
PHI=diff(PHI)./(-dz);

% Normalize 
A=repmat(nansum(PHI.^2.*dz,1)./H,[N 1]).^(1/2);
A(A==0)=Inf;
PHI=PHI./A;

A=repmat(nansum(PHI2.^2.*repmat(N2,[1 Nm])*dz,1)./(H*C.^2).',[N 1]).^(1/2);
A(A==0)=Inf;
PHI2=PHI2./A;
%PHI2(:,PHI2(1,:)<0)=-PHI2(:,PHI2(1,:)<0);

% Set surface PHI to positive value
%PHI2(:,PHI2(1,:)<0)=-PHI2(:,PHI2(1,:)<0);
PHI2(:,PHI(1,:)<0)=-PHI2(:,PHI(1,:)<0);
PHI(:,PHI(1,:)<0)=-PHI(:,PHI(1,:)<0);

if flag_0==1
	% Add in zeroth mode
	PHI=[ones([N 1]) PHI];
	PHI2=[zeros([N 1]) PHI2];
	C=[sqrt(9.81*H); C];
end

return



