% 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 [alpha,beta,norm_mx0,singvec0]=run_script_spatiotemporal_rapid(t_final)

[Re,N_1]=fix_all_parameters();
im=sqrt(-1);

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

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

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

dedz_1 = 2;

%**************************************************************************
% 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)/2;
    u_val=Re*z*(1-z);
    ddu_val=-2*Re;
    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 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)/2;
    du_val=Re*(1-2*z);
    u_val=Re*z*(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

% %**************************************************************************
% Construct matrix of energy norm.

[T0,T1]=get_tees(N_1);
d0T_uniform_L=           (1/2)*T0;
d1T_uniform_L=(dedz_1^2)*(1/2)*T1;

n_oss=2*(N_1+1);
nn=N_1+1;

B_mx=zeros(n_oss,n_oss);

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

alpha=1e-4:0.2:4;
beta=0:0.2:8;

% alpha=2;
% beta=0;

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

norm_mx0=zeros(length(alpha),length(beta));

for i=1:length(alpha)
    for j=1:length(beta)
        
        alpha_param=alpha(i);
        beta_param=beta(j);

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

        B_mx(1:nn,1:nn)=d1T_uniform_L+(k_val^2)*d0T_uniform_L;
        B_mx(nn+(1:nn),nn+(1:nn))=d0T_uniform_L;
        
        B_mx=(1/2)*(k_val^(-2))*B_mx;
        
        [~,L_aug,M_aug,~,~] =OSS_solver(alpha_param,beta_param,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);
        
        dt=1e-5;
        maxtt=floor(t_final/dt);
  
        Lplusaug=M_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;
        
        [V,D]=eig(A_mx,B_mx);
        vec=diag(D);
        vec=real(vec);
        [~,ix]=max(vec);
        
        singvec0=V(:,ix);
        
        norm_mx0(i,j)=max(vec);

    end
    
    display(strcat('i= ',num2str(i),'max-i-alpha=',num2str(length(alpha))))
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

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

