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

function [nrm,mysing0]=xrun_script_spatiotemporal_rapid_time(t_vec)

global A1
global A2
global a1
global a2

im=sqrt(-1);

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

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

%**************************************************************************
%  Collocation grids

kk_1_os=1:(N_1-3);
eta_1_os=cos((pi/(N_1-2))*kk_1_os);

kk_21_os=1:(N_21-3);
eta_21_os=cos((pi/(N_21-2))*kk_21_os);

kk_1_sq=1:(N_1-1);
eta_1_sq=cos((pi/(N_1-0))*kk_1_sq);

kk_21_sq=1:(N_21-1);
eta_21_sq=cos((pi/(N_21-0))*kk_21_sq);

dedz_1 = 2;
dedz_21 = 2/n_height;

%**************************************************************************
% derive OS matrices for the liquid

zeros_1_os=zeros(N_1-3,N_1+1);
T4_1=zeros_1_os;
T2_1=zeros_1_os;
T0_1=zeros_1_os;
T0_U2_1=zeros_1_os;
T2_U0_1=zeros_1_os;
T0_U0_1=zeros_1_os;

for i=1:(N_1-3)
    z=(eta_1_os(i)-1)*(1/2);
    u_val=U_1(z);
    ddu_val=ddU_1(z);
    for k=0:N_1       
        
        T4_1(i,k+1)=(dedz_1^4)*d4T(eta_1_os(i),k);
        T2_1(i,k+1)=(dedz_1^2)*d2T(eta_1_os(i),k);
        T0_1(i,k+1)=T(eta_1_os(i),k);
         
        T0_U2_1(i,k+1)=T(eta_1_os(i),k)*ddu_val;
        T2_U0_1(i,k+1)=(dedz_1^2)*d2T(eta_1_os(i),k)*u_val;
        T0_U0_1(i,k+1)=T(eta_1_os(i),k)*u_val;
          
    end
end

% derive OS matrices for the gas

zeros_21_os=zeros(N_21-3,N_21+1);
T4_21=zeros_21_os;
T2_21=zeros_21_os;
T0_21=zeros_21_os;
T0_U2_21=zeros_21_os;
T2_U0_21=zeros_21_os;
T0_U0_21=zeros_21_os;

for i=1:(N_21-3)
    z=(eta_21_os(i)+1)*(n_height/2);
    u_val=U_21(z);
    ddu_val=ddU_21(z);
    for k=0:N_21              
        
        T4_21(i,k+1)=(dedz_21^4)*d4T(eta_21_os(i),k);
        T2_21(i,k+1)=(dedz_21^2)*d2T(eta_21_os(i),k);
        T0_21(i,k+1)=T(eta_21_os(i),k);
        
        T0_U2_21(i,k+1)=T(eta_21_os(i),k)*ddu_val;
        T2_U0_21(i,k+1)=(dedz_21^2)*d2T(eta_21_os(i),k)*u_val;
        T0_U0_21(i,k+1)=T(eta_21_os(i),k)*u_val;                      
    end
end

%**************************************************************************
% derive squire matrices for the liquid

zeros_1_sq=zeros(N_1-1,N_1+1);
T2_1_sq=zeros_1_sq;
T0_1_sq=zeros_1_sq;
T0_U0_1_sq=zeros_1_sq;
cpl_1=zeros_1_sq;

for i=1:(N_1-1)
    z=(eta_1_sq(i)-1)*(1/2);
    du_val=dU_1(z);
    u_val=U_1(z);
    
    for k=0:N_1       
        
        T2_1_sq(i,k+1)=(dedz_1^2)*d2T(eta_1_sq(i),k);
        T0_1_sq(i,k+1)=T(eta_1_sq(i),k);
        
        T0_U0_1_sq(i,k+1)=T(eta_1_sq(i),k)*u_val;
        
        cpl_1(i,k+1)=-im*du_val*T(eta_1_sq(i),k);
          
    end
end

% derive squire matrices for the gas

zeros_21_sq=zeros(N_21-1,N_21+1);
T2_21_sq=zeros_21_sq;
T0_21_sq=zeros_21_sq;
T0_U0_21_sq=zeros_21_sq;
cpl_21=zeros_21_sq;

for i=1:(N_21-1)
    z=(eta_21_sq(i)+1)*(n_height/2);
    du_val=dU_21(z);
    u_val=U_21(z);
    
    for k=0:N_21       
        
        T2_21_sq(i,k+1)=(dedz_21^2)*d2T(eta_21_sq(i),k);
        T0_21_sq(i,k+1)=T(eta_21_sq(i),k);
        
        T0_U0_21_sq(i,k+1)=T(eta_21_sq(i),k)*u_val;
        
        cpl_21(i,k+1)=-im*du_val*T(eta_21_sq(i),k);
    end
end

% *************************************************************************
% *************************************************************************
%
% Construct matrix of energy norm.
%
% Liquid layer first
L_1=1;
L_21=n_height;
[T0,T1]=get_tees(N_1);
d0T_uniform_L=           (L_1/2)*T0;
d1T_uniform_L=(dedz_1^2)*(L_1/2)*T1;

% Gas layer now
[T0,T1]=get_tees(N_21);

d0T_uniform_G = (L_21/2)*T0;
d1T_uniform_G = (dedz_21^2)*(L_21/2)*T1;



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

nn1=N_1+1;
nn2=N_21+1;

n_oss=2*(nn1+nn2)+1;

B_mx=zeros(n_oss,n_oss);

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

alpha_param=1e-6;
beta_param=1;

k_val=sqrt(alpha_param*alpha_param+beta_param*beta_param);

B_mx(1:nn1,1:nn1)=(d1T_uniform_L+(k_val^2)*d0T_uniform_L);
B_mx(nn1+(1:nn2),nn1+(1:nn2))=r*(d1T_uniform_G+(k_val^2)*d0T_uniform_G);

B_mx(nn1+nn2+(1:nn1),nn1+nn2+(1:nn1))=d0T_uniform_L;
B_mx(2*nn1+nn2+(1:nn2),2*nn1+nn2+(1:nn2))=r*d0T_uniform_G;

% Typo in Yecko (2007), Equation 17.
B_mx(n_oss,n_oss)=(1/We)*k_val^4;

B_mx=B_mx/(2*k_val^2);

[~,L_aug,M_aug,~,~] =laminar_channel_solver_3deta(alpha_param,beta_param,...
    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);

[V,D]=eig(L_aug,M_aug);
lambda=diag(D);

for ii=1:length(lambda)
    if(abs(lambda(ii))>10000)
        lambda(ii)=-10000;
    end

    B=isnan(lambda(ii));
    if(B==1)
         lambda(ii)=-10000;
    end
end

[~,ix]=sort(-real(lambda));
Kc=20;
lc=real(lambda(ix(Kc)));

myD=zeros(length(lambda));
myD0=zeros(length(lambda));
% myD=diag(exp(lambda*t_final));

[~,myn]=size(L_aug);
id_m=eye(myn,myn);
Q=zeros(myn,myn);

for ii=1:myn
    for jj=1:myn
        ei=id_m(:,ii);
        xj=V(:,jj);
        Q(ii,jj)=dot(ei,xj);
    end
end
                
nrm=0*t_vec;

for i=1:length(t_vec)
    
        t_final=t_vec(i);
        
        %         for ii=1:length(lambda)
        %             if(real(lambda(ii))-lc>=0)
        %                 myD(ii,ii)=exp(lambda(ii)*t_final);
        %                 myD0(ii,ii)=1;
        %             end
        %         end
        
        % Etaug=Q*myD*(Q^(-1));
        
        dt=1e-4;
        maxtt=floor(t_final/dt);
  
        % Backward differentiation much better than trapezoidal rule
        Lplusaug=M_aug+0*(dt/2)*L_aug;
        Lminusaug=M_aug-(dt)*L_aug;
        Lminusaug_inv=Lminusaug^(-1);
        evolopaug=Lminusaug_inv*Lplusaug;

        Etaug=(evolopaug^maxtt);
       
        A_mx=ctranspose(Etaug)*B_mx*Etaug;
        A_mx=(A_mx+ctranspose(A_mx))/2;
        % B_mx0=myD0*B_mx*myD0;
        % B_mx0=(B_mx0+ctranspose(B_mx0))/2;
        
        [V,D]=eig(A_mx,B_mx);
        vec=diag(D);
        vec=real(vec);
        
        for ii=1:length(vec)
            B=isinf(vec(ii));
            if(B==1)
                vec(ii)=-1000;
            end
        end
        
        [my_max,ix]=max(vec);
        mysing0=V(:,ix);
        nrm(i)=my_max;
    
        display(i)
end


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


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 = ddU_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 = ddU_21(~)

global A2

d2u=2*A2;
    
end