function [sig,w,b]=SSF(z,U,Bz,k,l,nu,kappa,bcw,bcb,imode)
%
% USAGE: [sig,w,b]=SSF(z,U,Bz,k,l,nu,kappa,bcw,bcb,imode)
%
% Stability analysis for a viscous, diffusive, stratified, parallel shear flow
% INPUTS:
% z = vertical coordinate vector (evenly spaced)
% U = velocity profile
% Bz = buoyancy gradient profile (Bz=squared BV frequency)
% k,l = wave vector 
% nu, kappa = viscosity, diffusivity
% bcw = boundary conditions on w at (1) z=z(0) and (2) z=z(N+1)
%       r=rigid, f=frictionless (default)
% bcb = boundary conditions on b at (1) z=z(0) and (2) z=z(N+1)
%       f=insulating, c=fixed-buoyancy (default)
% imode = mode choice (default imode=1)
%         imode=0: output all modes, sorted by growth rate
%
% OUTPUTS:
% sig = growth rate of FGM
% w = vertical velocity eigenfunction
% b = buoyancy eigenfunction
%
% CALLS:
% ddz2, ddz4
%
% CITATION:
%   Instability in Geophysical Flows
%   W.D. Smyth & J.R. Carpenter
%   Cambridge University Press, 2019
%
% W. Smyth, Oregon State University, February 2016

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stage 1: Preliminaries
%

% check for equal spacing
if abs(std(diff(z))/mean(diff(z)))>.000001
    disp(['SSF: values not evenly spaced!'])
    sig=NaN;
    return
end

% defaults
if(~exist('bcw'))
    % bcs for w are both frictionless
    bcw='ff'; 
end
if(~exist('bcb'))
    %boundary conditions for b are both constant-buoyancy
    bcb='cc'; 
end
if ~exist('imode');imode=1;end

% define constants
ii=complex(0.,1.);
del=mean(diff(z));N=length(z);
kt=sqrt(k^2+l^2);

% differentate U
Uzz=ddz2(z)*U;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stage 2: Derivative matrices and BCs

% D2: 2nd derivative with impermeable boundary
D2=ddz2(z);  
D2(1,:)=0;D2(1,1:2)=[-2 1]/del^2;
D2(N,:)=0;D2(N,N-1:N)=[1 -2]/del^2;

% Fourth derivative with frictionless or rigid boundaries
D4=ddz4(z);
if bcw(1)=='f';    % frictionless
    D4(1,:)=0;D4(1,1:3)=[5 -4 1]/del^4;
elseif bcw(1)=='r';     % rigid
    D4(1,:)=0;D4(1,1:3)=[7 -4 1]/del^4; 
else
    display(['SSF: Boundary condition 1 on w not understood'])
    sig=nan;
    return
end
D4(2,:)=0;D4(2,1:4)=[-4 6 -4 1]/del^4;
if bcw(2)=='f'; % frictionless
    D4(N,:)=0;D4(N,N-2:N)=[1 -4 5]/del^4;
elseif bcw(2)=='r'
    D4(N,:)=0;D4(N,N-2:N)=[1 -4 7]/del^4;
else
    display(['SSF: Boundary condition 2 on w not understood'])
    sig=nan;
    return
end
D4(N-1,:)=0;D4(N-1,N-3:N)=[1 -4 6 -4]/del^4;

% 
% 2nd derivative matrix for buoyancy 

D2b=ddz2(z);
if bcb(1)=='c'
    % Dirichlet bcs
    D2b(1,:)=0;D2b(1,1:2)=[-2 1]/del^2;
elseif bcb(1)=='f'  
    % Neumann bcs
    D2b(1,:)=0;D2b(1,1:2)=[-1 1]*2/3/del^2;
else
    display(['ddz2: Boundary condition 1 on buoyancy not understood'])
    d=NaN;
    return
end

if bcb(2)=='c'
    % Dirichlet bcs
    D2b(N,:)=0;D2b(N,N-1:N)=[1 -2]/del^2;
elseif bcb(2)=='f'  
    % Neumann bcs
    D2b(N,:)=0;D2b(N,N-1:N)=[1 -1]*2/3/del^2;
else
    display(['ddz2: Boundary condition 2 not understood'])
    d=NaN;
    return
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stage 3: Assemble stability matrices
%
% Laplacian and squared Laplacian matrices
Id=eye(N);
L=D2-kt^2*Id;
Lb=D2b-kt^2*Id;
LL=D4-2*kt^2*D2+kt^4*Id;

N2=2*N;
NP=N+1;
A=zeros(N2,N2);B=zeros(N2,N2);

% assemble matrix A
A=[L Id*0 ; Id*0 Id];

% Compute submatrices of B using Levi's syntax
b11=-ii*k*diag(U)*L+ii*k*diag(Uzz)+nu*LL;
b21=-diag(Bz);
b12=-kt^2*Id;
b22=-ii*k*diag(U)+kappa*Lb;

% assemble matrix B
B=[b11 b12 ; b21 b22];
                        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stage 4: Solve eigenvalue problem and extract results
%
% Solve generalized eigenvalue problem
[v,e]=eig(B,A);sigma=diag(e);

% Sort eigvals
[sr,ind]=sort(real(sigma),1,'descend');
sigma=sigma(ind);
v=v(:,ind);


% Extract the selected mode(s)
if imode==0
    sig=sigma;
    w=v(1:N,:);
    b=v(NP:N2,:);
elseif imode>0
    sig=sigma(imode);
    w=v(1:N,imode);
    b=v(NP:N2,imode);
end

% sig=sigma(imode);
% w=v(1:N,imode);
% b=v(NP:N2,imode);

return