% ========================================================================
% Linearized Navier-Stokes equations with kx = 0. 
% Find the eigenvalues of TTs using chebfun and integral formulation.
%
% additional m-file needed: svdfr
%
% Original equations:
%
%   x1'''' - (2*kz*kz + ii*om*R)*x1'' + kz*kz*(kz*kz + ii*om*R)*x1 = ...
%       R*kz*kz*d2 + R*ii*kz*d3';
%
%   x2'' - (kz*kz + ii*om*R)*x2 - R*ii*kz*U'*x1 = -R*ii*kz*d1;
%
%   u = -(ii/kz)*x2;
%   v = x1;
%   w = (ii/kz)*x1';
%
%   BCs:
%   x1(-1) = x1(+1) = 0; x1'(-1) = x1'(+1) = 0
%   x2(-1) = x2(+1) = 0
%
% written by: Binh Lieu
% ========================================================================

%% system parameters
baseflow = 1; % baseflow

% System parameters:
R = 2000; % Reynolds number
 
kzval = linspace(0.1,5,30); % spanwise wave-number
kzgrd = length(kzval);
 
om = 0; % temporal frequency
 
dom = domain(-1,1);     % domain of your function
y = chebfun('y',dom);
fone = chebfun(1,dom);      % fone(y) = 1
fzero = chebfun(0,dom);     % fzero(y) = 0
 
%% based flow
if baseflow
    % % matrix representation of Umean, Uy, Uyy
    % in Poiseuille flow
    U = diag(1 - y.^2);
    Uy = diag(-2*y);
    Uyy = diag(-2*fone);
else
    % in Couette flow
    U = diag(y);
    Uy = diag(fone);
    Uyy = diag(fzero);
end
 
%% Boundary conditions
Wa0{1} = [1, 0, 0, 0; 0, 1, 0, 0];
Wa0{2} = [1, 0];
Wb0{1} = [1, 0, 0, 0; 0, 1, 0, 0];
Wb0{2} = [1, 0];
 
%% Looping over kz
Smax = zeros(kzgrd,1); A0 = cell(2,2); B0 = cell(2,3); C0 = cell(3,2);
for indz = 1:kzgrd
    
    kz = kzval(indz); kz2 = kz*kz; kz4 = kz2*kz2;
    
    % coefficients of the operator A0
    % kz^2*(kz^2/R + i*om)*phi1 + 0*phi1' - (2*kz^2*R + i*om)*phi1'' + 
    % 0*phi1''' + (1/R)*phi1'''' + 0*phi2 + 0*phi2' + 0*phi2''
    A11 = [kz2*(kz2/R + 1i*om)*fone, fzero, ...
        -(2*kz2/R + 1i*om)*fone, fzero, (1/R)*fone];
    A12 = [fzero, fzero, fzero];
    % -i*kz*U'*phi1 + 0*phi1' + 0*phi1'' + 0*phi1''' + 0*phi1'''' -
    % (kz^2/R + i*om)*phi2 + 0*phi2' + (1/R)*phi2''
    A21 = [-1i*kz*Uy*fone, fzero, fzero, fzero, fzero];
    A22 = [-(kz2/R + 1i*om)*fone, fzero, (1/R)*fone];
        
    A0 = {A11, A12; A21, A22};
    
    % coefficients of the operator B0
    % 0*d1 + 0*d1' - i*kz*d2 + 0*d2' + 0*d3 + 1*d3'
    B11 = [fzero, fzero]; B12 = [-1i*kz*fone, fzero]; 
    B13 = [fzero, fone];
    % -1*d1 + 0*d1' + 0*d2 + 0*d2' + 0*d3 + 0*d3'
    B21 = [-fone, fzero]; B22 = [fzero, fzero]; B23 = [fzero, fzero];
        
    B0 = {B11, B12, B13; B21, B22, B23};
    
    % coefficients of the operator C0
    % u = 0*phi1 + 0*phi1' + 1*phi2 + 0*phi2'
    C11 = [fzero, fzero]; C12 = [fone, fzero];
    % v = i*kz*phi1 + 0*phi1' + 0*phi2 + 0*phi2'
    C21 = [1i*kz*fone, fzero]; C22 = [fzero, fzero];
    % w = 0*phi1 - 1*phi1' + 0*phi2 + 0*phi2'
    C31 = [fzero, -fone]; C32 = [fzero, fzero];
        
    C0 = {C11, C12; C21, C22; C31, C32};
        
    % solving for the left principal singular pair
    [Sfun, Sval] = svdfr(A0,B0,C0,Wa0,Wb0,1,1);

    % saving the largest singular value for each value of kz        
    Smax(indz) = Sval(1);
    
    [indz]

end
 
% Plotting the largest singular value
plot(kzval,Smax,'-','LineWidth',1.1);
xlab = xlabel('k_z', 'interpreter', 'tex');
set(xlab, 'FontName', 'cmmi10', 'FontSize', 20);
h = get(gcf,'CurrentAxes');
set(h,'FontName','cmr10','FontSize',15,'xscale','log','yscale','lin');
