% Copyright (c) 2015, Lennon O'Naraigh
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
% 
% 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
% 
% 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer 
% in the documentation and/or other materials provided with the distribution.
% 
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, 
% INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, 
% OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 
% OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 
% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, 
% EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

%  30/01/2015 Solves the **OSS** equation for a Poiseuille flow in a channel.

function [lambda,L_aug,M_aug,L_os,M_os] =OSS_solver(alpha,beta,T4_1,T2_1,T0_1,T0_U2_1,T2_U0_1,...
    T0_U0_1,T2_1_sq,T0_1_sq,T0_U0_1_sq,cpl_1)

global im
global Re_global
global N_1_global

[Re,N_1]=fix_all_parameters();

Re_global=Re;
N_1_global=N_1;

im=sqrt(-1);

%**************************************************************************

cpl_1=cpl_1*beta;

%**************************************************************************
%  Define the eigenvalue problem

ksquare=alpha*alpha+beta*beta;

L_1_os= (1/Re)*(T4_1 -2*ksquare*T2_1+ (ksquare^2)*T0_1 )  + im*alpha*T0_U2_1  - im*alpha*(T2_U0_1 -ksquare*T0_U0_1 );
M_1_os =  T2_1 -  ksquare*T0_1;
L_1_sq=(1/Re)*(T2_1_sq -ksquare*T0_1_sq )-im*alpha*T0_U0_1_sq;
M_1_sq =T0_1_sq;

%**************************************************************************

[L_aug,M_aug,L_os,M_os] = make_augmented_matrices(L_1_os,L_1_sq,M_1_os,M_1_sq,cpl_1);


lambda=eig(L_aug,M_aug);

for i=1:length(lambda)
    if(abs(lambda(i))>10000)
        lambda(i)=-10000;
    end
    
    B=isnan(lambda(i));
    if(B==1)
         lambda(i)=-10000;
    end

end

%**************************************************************************

end

%**************************************************************************
%**************************************************************************
%**************************************************************************
%**************************************************************************

function [L_aug,M_aug,L_os,M_os] = make_augmented_matrices(L_1_os,L_1_sq,M_1_os,M_1_sq,cpl_1)

global N_1_global

N_1=N_1_global;

N_aug = 2*(N_1+1);

L_aug = zeros(N_aug,N_aug);
M_aug = zeros(N_aug,N_aug);

%**************************************************************************
% Wall-normal velocity block

for j=1:N_1+1
   %  BC 1
   L_aug(1,j)=T(-1,j-1); 
   %  BC 2
   L_aug(2,j)=d1T(-1,j-1);
end

for i=1:(N_1-3)
    for j=1:(N_1+1)
        L_aug(i+2,j)=L_1_os(i,j);
        M_aug(i+2,j)=M_1_os(i,j);
    end
end

for j=1:N_1+1
   %  BC 4
   L_aug(N_1,j)=T(1,j-1);
   %  BC 4
   L_aug(N_1+1,j)=d1T(1,j-1);
end

%**************************************************************************
% Wall-normal vorticity block

nn=N_1+1;

for j=1:N_1+1
   L_aug(nn+1,nn+j)=T(-1,j-1); 
end

for i=1:(N_1-1)
    for j=1:(N_1+1)
        L_aug(nn+i+1,nn+j)=L_1_sq(i,j);
        M_aug(nn+i+1,nn+j)=M_1_sq(i,j);
    end
end

for j=1:N_1+1
   L_aug(2*nn,nn+j)=T(1,j-1); 
end

%**************************************************************************
% Coupling:

for i=1:(N_1-1)
    for j=1:(N_1+1)
        L_aug(nn+1+i,j)=cpl_1(i,j);
    end
end

%**************************************************************************

L_os=L_aug(1:nn,1:nn);
M_os=M_aug(1:nn,1:nn);

end

%**************************************************************************
%**************************************************************************
%**************************************************************************
%**************************************************************************

function t = T(z,k)

t=cos(k*acos(z));

end

%**************************************************************************

function dt = d1T(z,k)

% specifiy dt at boundary

if(abs(z)==1)
    if(z==1)
        dt=k*k;
    else
        dt=((-1)^(k+1))*k*k;
    end
else
    dt = (k/sqrt(1-z*z))*sin(k*acos(z));
end

end

%**************************************************************************

function d2t = d2T(z,k)

t=cos(k*acos(z));

if(abs(z)==1)
    if(z==1)
        d2t=(1/3)*(k^4-k^2);
    else
        d2t=(1/3)*((-1)^k)*(k^4-k^2);
    end
else
    dt = (k/sqrt(1-z*z))*sin(k*acos(z));
    d2t=(z/(1-z*z))*dt - (k*k/(1-z*z))*t;
end

end

%**************************************************************************

% % function d3t = d3T(z,k)
% % 
% % t=cos(k*acos(z));
% % 
% % if(abs(z)==1)
% %     if(z==1)
% %         d3t=(1/(3*5))*k*k*(k*k-1)*(k*k-2*2);
% %     else
% %         d3t=(1/(3*5))*((-1)^(k+3))*k*k*(k*k-1)*(k*k-2*2);
% %     end
% % else
% %     dt = (k/sqrt(1-z*z))*sin(k*acos(z));
% %     d2t=(z/(1-z*z))*dt - (k*k/(1-z*z))*t;
% %     d3t = (3*z/(1-z*z))*d2t + ((1-k*k)/(1-z*z))*dt;
% % end
% % 
% % end

%**************************************************************************

function d4t = d4T(z,k)

t=cos(k*acos(z)); 
dt = (k/sqrt(1-z*z))*sin(k*acos(z));                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           
d2t=(z/(1-z*z))*dt - (k*k/(1-z*z))*t;
d3t = (3*z/(1-z*z))*d2t + ((1-k*k)/(1-z*z))*dt;

d4t = (5*z/(1-z*z))*d3t + ((4-k*k)/(1-z*z))*d2t;

end

%**************************************************************************
%**************************************************************************


