% ========================================================================
% Orr-Sommerfeld equation (streamfunction formulation)
% Find the eigenvalues of TTs using chebfun and integral formulation.
%
% additional m-file needed: svdfr
%
% Original equations:
%
%   (1/R)*x1'''' + a2*x1'' + a0*x1 = -d1' + 1i*kx*d2;
%
%   u = x1';
%   v = -1i*kx*x1;
%
%   a2(y) = -(2*(1/R)*kx^2 + 1i*kx*U(y) + 1i*om)
%   a0(y) = (1/R)*kx^4 + 1i*kx^3*U(y) + 1i*kx*U''(y) + 1i*om*kx^2
%
%   BCs:
%   x1(-1) = x1(+1) = 0; x1'(-1) = x1'(+1) = 0
%
% written by: Binh Lieu
% ========================================================================

%% system parameters
kxval = logspace(-2,1,30);     % streamwise wavenumber
% kxval = 0.5;
kxgrd = length(kxval);

% omval = -logspace(-4,1,50);     % temporal frequency
omval = -0.3;
omgrd = length(omval);

R = 2000; % Reynolds number

baseflow = 1; % baseflow

%% domain of interest
dom = domain(-1,1);
y = chebfun('y',dom);

% one and zero functions
fone = chebfun(1,[-1 1]);
fzero = chebfun(0,[-1 1]);

%% 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];
Wb0{1} = [1, 0, 0, 0; 0, 1, 0, 0];

%% Looping over kx and om
A0 = cell(1,1); B0 = cell(1,2); C0 = cell(2,1);
Smax = zeros(kxgrd,omgrd);
for indx = 1:kxgrd
    kx = kxval(indx); kx2 = kx*kx; kx4 = kx2*kx2;
    for indom = 1:omgrd
        om = omval(indom);
        
        % coefficients of the A0 operator
        % a0*phi + 0*phi? + a2*phi?? + 0*phi??? + 1/R*phi????
        a2 = -(2*kx2*(1/R)*fone + 1i*kx*U*fone + 1i*om*fone);
        a0 = (1/R)*kx4*fone + 1i*kx*kx2*U*fone + 1i*kx*Uyy*fone + ...
            1i*om*kx2*fone;
        A0{1} = [a0, fzero, a2, fzero, (1/R)*fone];
        
        % coefficients of the B0 operator
        % 0*d1 - d1? + i*kx*d2 + 0*d2?
        B0{1,1} = [fzero, -fone]; B0{1,2} = [1i*kx*fone, fzero];
        
        % coefficients of the C0 operator
        % u = 0*phi + 1*phi?; v = -i*kx*phi + 0*phi?
        C0{1,1} = [fzero, fone]; C0{2,1} = [-1i*kx*fone, fzero];
        
        % solving for the left principal singular pair
        [Sfun, Sval] = svdfr(A0,B0,C0,Wa0,Wb0,1,1);
        
        % saving the largest singular value for each pair of kx and om
        Smax(indx,indom) = Sval(1);
        
        [indx,indom]
    end
end

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