% 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.

%  January 2015.  Solves the OS/Squire equations in a CHANNEL Yecko (IJMF,
%  2007)

function [lambda,L_aug,M_aug,L_os,M_os,L_sq,M_sq] =laminar_channel_solver_3deta(alpha,beta,...
                    T4_1,T2_1,T0_1,T0_U2_1,T2_U0_1,T0_U0_1,...
                    T4_21,T2_21,T0_21,T0_U2_21,T2_U0_21,T0_U0_21,...
                    T2_1_sq,T0_1_sq,T0_U0_1_sq,cpl_1,...
                    T2_21_sq,T0_21_sq,T0_U0_21_sq,cpl_21)

global n_height
global Re
global m
global r
global We
global Fr

global A1
global A2
global a1
global a2

global dedz_1
global dedz_21

global N_1
global N_21

global im

%  The complex unit

im=sqrt(-1);

%**************************************************************************
%  Physical parameters

[n_height,Re,m,r,N_1,N_21,We,Fr]=fix_all_parameters();
[~,~,~,~,~,~,A1,A2,a1,a2]=get_u_laminar();

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

dedz_1 = 2;
dedz_21 = 2/n_height;

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

cpl_1=cpl_1*beta;
cpl_21=cpl_21*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 );
L_21_os=(m/(r*Re))*(T4_21-2*ksquare*T2_21+(ksquare^2)*T0_21)  + im*alpha*T0_U2_21 - im*alpha*(T2_U0_21-ksquare*T0_U0_21);

M_1_os =  T2_1 -  ksquare*T0_1;
M_21_os = T2_21 - ksquare*T0_21;

L_1_sq= (1/Re)    *(T2_1_sq -ksquare*T0_1_sq )-im*alpha*T0_U0_1_sq;
L_21_sq=(m/(r*Re))*(T2_21_sq-ksquare*T0_21_sq)-im*alpha*T0_U0_21_sq;

M_1_sq =T0_1_sq;
M_21_sq=T0_21_sq;

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

[L_aug,M_aug,L_os,M_os,L_sq,M_sq] = make_augmented_matrices(L_1_os,L_21_os,L_1_sq,L_21_sq,M_1_os,M_21_os,M_1_sq,M_21_sq,cpl_1,cpl_21,alpha,beta);

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,L_sq,M_sq] = make_augmented_matrices(L_1,L_21,L_1_sq,L_21_sq,M_1,M_21,M_1_sq,M_21_sq,cpl_1,cpl_21,alpha,beta)

global m
global r
global Re
global We

global N_1
global N_21

global dedz_1
global dedz_21

global im

N_aug = 2*((N_1+1) + (N_21+1))+1;
N_os=(N_1+1) + (N_21+1);
ksquare=alpha*alpha+beta*beta;

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

% L_os=zeros(N_os,N_os);
% M_os=zeros(N_os,N_os);

%**************************************************************************
% Liquid 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(i,j);
        M_aug(i+2,j)=M_1(i,j);
    end
end

%**************************************************************************
% Liquid - gas interfacial conditions

u0 = U_1(0);
du_1 = dU_1(0);
du_21 = dU_21(0);

for j=1:(N_1+1)
    %  IC 1
    L_aug(N_1,j)=T(1,j-1);
    
    %  IC 2
    L_aug(N_1+1,j)=dedz_1*d1T(1,j-1);
   
    %  IC 3
    L_aug(N_1+2,j)=(  (dedz_1^2)*d2T(1,j-1)+ksquare*T(1,j-1)   ); %(ddu_21-m_0i*ddu_1)*T(1,j-1);
    % M_aug(N_1+2,j)=(im/alpha)*1*((dedz_1^2)*d2T(1,j-1)+ksquare*T(1,j-1));

    %  IC 4
    L_aug(N_1+3,j)=-im*alpha*Re*1*u0*dedz_1*d1T(1,j-1)+ 1*((dedz_1^3)*d3T(1,j-1) - 3*ksquare*dedz_1*d1T(1,j-1)) + im*alpha*1*Re*du_1*T(1,j-1)+ im*(ksquare/alpha)*((ksquare*(Re/We))/(du_1-du_21))*dedz_1*d1T(1,j-1);
    M_aug(N_1+3,j)=1*Re*dedz_1*d1T(1,j-1);  
end

for j=1:(N_21+1)
    %  IC 1
    L_aug(N_1,N_1+1+j)=-T(-1,j-1);
    
    %  IC 2
    L_aug(N_1+1,N_1+1+j)= -dedz_21*d1T(-1,j-1);
    
    %  IC 3
    L_aug(N_1+2,N_1+1+j)=-m*(  (dedz_21^2)*d2T(-1,j-1)+ksquare*T(-1,j-1)  );
    % M_aug(N_1+2,N_1+1+j)=-(im/alpha)*m*((dedz_21^2)*d2T(-1,j-1)+ksquare*T(-1,j-1));

    
    %  IC 4
    temp_var = -im*alpha*Re*r*u0*dedz_21*d1T(-1,j-1)+ m*((dedz_21^3)*d3T(-1,j-1) - 3*ksquare*dedz_21*d1T(-1,j-1)) + im*alpha*r*Re*du_21*T(-1,j-1)+im*(ksquare/alpha)*((ksquare*(Re/We))/(du_1-du_21))*dedz_21*d1T(-1,j-1); 
    L_aug(N_1+3,N_1+1+j)=-temp_var;
    M_aug(N_1+3,N_1+1+j)=-r*Re*dedz_21*d1T(-1,j-1);      
end

L_aug(N_1+1,N_aug)=-im*alpha*(du_1-du_21);

%**************************************************************************
% Gas block

for i=1:(N_21-3)
    for j=1:(N_21+1)
        L_aug(i+N_1+3,j+N_1+1)=L_21(i,j);
        M_aug(i+N_1+3,j+N_1+1)=M_21(i,j);
    end
end

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

% Wall-gas BCs

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


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


%**************************************************************************
%**************************************************************************
% Squire part now:

nn=N_1+N_21+2;

% Liquid - Bottom wall BC

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

% Liquid block

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

% Liquid - gas interfacial conditions
% Liquid side

for j=1:N_1+1
   L_aug(nn+N_1+1,nn+j)= T(1,j);
   % M_aug(nn+N_1+1,nn+j)= 0?;
   
   % Coupling (liquid side only)
   % L_aug(nn+N_1+1,j)=im*beta*(du_1-du_21)*T(1,j);
   
   L_aug(nn+N_1+2,nn+j)=1*dedz_1*d1T(1,j);
   % M_aug(nn+N_1+2,nn+j)=0;
end

% Liquid - gas interfacial conditions
% Gas  side

for j=1:N_21+1
   L_aug(nn+N_1+1,nn+N_1+1+j)= -T(-1,j);
   % M_aug(nn+N_1+1,nn+N_1+1+j)=  T(-1,j);
   
   L_aug(nn+N_1+2,nn+N_1+1+j)=m*dedz_21*d1T(-1,j);
   % M_aug(nn+N_1+2,nn+N_1+1+j)=0;
end

L_aug(nn+N_1+1,N_aug)=-im*beta*(du_21-du_1);

% Gas block

for i=1:(N_21-1)
    for j=1:(N_21+1)
        L_aug(nn+N_1+2+i,nn+N_1+1+j)=L_21_sq(i,j);
        M_aug(nn+N_1+2+i,nn+N_1+1+j)=M_21_sq(i,j);
    end
end

% Gas - top wall BC

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

% Kinematic condition:

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

L_aug(N_aug,N_aug)=-im*alpha*u0;
M_aug(N_aug,N_aug)=1;


%**************************************************************************
%**************************************************************************
% 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

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


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

N_sq=(N_1+1) + (N_21+1);
L_sq=L_aug(N_os+(1:N_sq),N_os+(1:N_sq));
M_sq=M_aug(N_os+(1:N_sq),N_os+(1:N_sq));


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

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

function u = U_1(z1)

global A1
global a1

u=A1*(z1^2)+a1*z1+1;

end

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

function du = dU_1(z1)

global A1
global a1

du=2*A1*z1+a1;

end

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

function d2u = d2U_1(~)

global A1

d2u=2*A1;

end

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

function u = U_21(z2)

global A2
global a2

u=A2*(z2^2)+a2*z2+1;
    
end

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

function du = dU_21(z2)

global A2
global a2

du=2*A2*z2+a2;
    
end

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

function d2u = d2U_21(~)

global A2

d2u=2*A2;
    
end