clc
clear

% Initialize
rootdir = 'I:\GRL_DataOCF12000\';
imname = 'ocf12000_st_'; 
startim = 1500; 
overspace=1;
win = round(1152/1); % [px] window size for psd
w=60; % FIV windows
nu=1e-6;
U_b = 0.12; % (m/s) velocity reference scale
frsec = 0.25*0.083; % (s) | time between each frame.
mpx = 0.24/1152; % (m/px) | m to pixel scaleing ratio
mpy = 0.06/1152; % (m/px) | m to pixel scaleing ratio y axis
m2px = 1/mpx*frsec; % [(m/s)]*(px/m)*(s) | m/s to pixels convertion ratio
m2py = 1/mpy*frsec; %  [(m/s)]*(px/m)*(s) | m/s to pixels convertion ratio
k_DNS = 1/mpx;
f_DNS = 1/frsec;
windo = 5; % Window Filer | 1 rectangular | 2 Hamming | 3 Hanning | 4 Flat top | 5 Blackman
load([rootdir,'data\imagedata_w60.mat'])
win_FIV = round(length(imagedata{1,1}(1,:)))-1;
k_FIV = 1/(mpx*s{2,2});
FIV_order=5;
DNS_order=5;

% Load & Process Data
E11_DNS=0;
E22_DNS=0;
e0_DNS_E11=0;
e0_DNS_E22=0;
R11_DNS=0;
R22_DNS=0;
urms_DNS=0;
uprime_rms_DNS=0;
vrms_DNS=0;
vprime_rms_DNS=0;

E11_FIV=0;
E22_FIV=0;
e0_FIV_E11=0;
e0_FIV_E22=0;
R11_FIV=0;
R22_FIV=0;
urms_FIV=0;
uprime_rms_FIV=0;
vrms_FIV=0;
vprime_rms_FIV=0;
count=0;
if overspace==1
    for i = 0:7
        % Load in data
        load([rootdir,'mat\',imname,num2str(startim+i*600),'.mat']);
        u=u*U_b;
        v=v*U_b;
        dev = 13;
        for j = 1:dev-1
            count = count + 1;
            %~~~~~~~~~~~~~~~~~~~~~~~~~~
            %--------Initialize-------|
            %~~~~~~~~~~~~~~~~~~~~~~~~~~
            U(:,1) = u(round(end*j/dev),:)';
            V(:,1) = v(round(end*j/dev),:)';
%             U(:,1) = u(:,round(end*j/dev));
%             V(:,1) = v(:,round(end*j/dev));
            % FIV interpolated data {1=inpaint,2=linear,3=cubic}
            U_FIV(:,1) = imagedatainterpol{i*24+1,1}(round(end*j/dev),:)'./m2px;
            V_FIV(:,1) = imagedatainterpol{i*24+1,2}(round(end*j/dev),:)'./m2py;
            % RMS velocities
            urms_DNS_n=rms(U);
            uprime_rms_DNS_n=rms(U-mean(U));
            vrms_DNS_n=rms(V);
            vprime_rms_DNS_n=rms(V-mean(V));
            urms_FIV_n=rms(U_FIV);
            uprime_rms_FIV_n=rms(U_FIV-mean(U_FIV));
            vrms_FIV_n=rms(V_FIV);
            vprime_rms_FIV_n=rms(V_FIV-mean(V_FIV));
            %~~~~~~~~~~~~~~~~~~~~~~~~~~
            %---1-D Energy Spectrum---|
            %~~~~~~~~~~~~~~~~~~~~~~~~~~
            % DNS spectrum
            [E11_DNS_n,E22_DNS_n,K1_DNS,R11_DNS_n,R22_DNS_n] = EK_spectra(U,V,k_DNS,win,windo);
            % FIV spectrum
            [E11_FIV_n,E22_FIV_n,K1_FIV,R11_FIV_n,R22_FIV_n] = EK_spectra(U_FIV(:,1),V_FIV(:,1),k_FIV,win_FIV,windo);
            %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            %---1-D Instantanious Energy Dissipation---|
            %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            Udiff=gradient(U,mpx);
            index=log10(Udiff(:).^2)==-Inf;
            e0_DNS_E11_n=15*nu*10.^mean(log10(Udiff(index==0).^2));
            
            Vdiff=gradient(V,mpx);
            index=log10(Vdiff(:).^2)==-Inf;
            e0_DNS_E22_n=15/2*nu*10.^mean(log10(Vdiff(index==0).^2));
            
            Udiff=gradient(U_FIV,s{2,2}*mpx);
            index=log10(Udiff(:).^2)==-Inf;
            e0_FIV_E11_n=15*nu*10.^mean(log10(Udiff(index==0).^2));
        
            Vdiff=gradient(V_FIV,s{2,2}*mpx);
            index=log10(Vdiff(:).^2)==-Inf;
            e0_FIV_E22_n=15/2*nu*10.^mean(log10(Vdiff(index==0).^2));
            %~~~~~~~~~~~~~~~~~~~~~~~~~~
            %--------Data adding-------|
            %~~~~~~~~~~~~~~~~~~~~~~~~~~
            E11_DNS=E11_DNS+E11_DNS_n;
            E22_DNS=E22_DNS+E22_DNS_n;
            e0_DNS_E11=e0_DNS_E11+e0_DNS_E11_n;
            e0_DNS_E22=e0_DNS_E22+e0_DNS_E22_n;
            R11_DNS=R11_DNS+R11_DNS_n;
            R22_DNS=R22_DNS+R22_DNS_n;
            urms_DNS=urms_DNS+urms_DNS_n;
            uprime_rms_DNS=uprime_rms_DNS+uprime_rms_DNS_n;
            vrms_DNS=vrms_DNS+vrms_DNS_n;
            vprime_rms_DNS=vprime_rms_DNS+vprime_rms_DNS_n;
        
            E11_FIV=E11_FIV+E11_FIV_n;
            E22_FIV=E22_FIV+E22_FIV_n;
            e0_FIV_E11=e0_FIV_E11+e0_FIV_E11_n;
            e0_FIV_E22=e0_FIV_E22+e0_FIV_E22_n;
            R11_FIV=R11_FIV+R11_FIV_n;
            R22_FIV=R22_FIV+R22_FIV_n;
            urms_FIV=urms_FIV+urms_FIV_n;
            uprime_rms_FIV=uprime_rms_FIV+uprime_rms_FIV_n;
            vrms_FIV=vrms_FIV+vrms_FIV_n;
            vprime_rms_FIV=vprime_rms_FIV+vprime_rms_FIV_n;
        end
    end
else
    counti=0;
    win=180;
    for i = T(1):25:T(end)-25
        counti=counti+1;
        % Load in data
        % DNS
        load([rootdir,'mat\',imname,num2str(i),'.mat']);
        ydist=1:48:1152;
        xdist=1:73:1152;
        [y_FIV,x_FIV] = size(imagedatainterpol{1,1});
        countj=0;
        for j=xdist
            countj=countj+1;
            countk=0;
            for k=ydist
                countk=countk+1;
                Ublock_DNS(countj,countk,counti)=u(j,k)*U_b;
                Vblock_DNS(countj,countk,counti)=v(j,k)*U_b;
            end
        end
        % FIV interpolated data {1=inpaint,2=linear,3=cubic}
        Ublock_FIV(:,:,counti)=imagedatainterpol{counti,1}./m2px;
        Vblock_FIV(:,:,counti)=imagedatainterpol{counti,2}./m2py;
    end
    count=0;
    for j = 1:countj
        for k = 1:countk
            count = count + 1;
            %~~~~~~~~~~~~~~~~~~~~~~~~~~
            %--------Initialize-------|
            %~~~~~~~~~~~~~~~~~~~~~~~~~~
            U = Ublock_DNS(j,k,:);
            U = U(:);
            V = Vblock_DNS(j,k,:);
            V = V(:);
            U_FIV = Ublock_FIV(j,k,:);
            U_FIV = U_FIV(:);
            V_FIV = Vblock_FIV(j,k,:);
            V_FIV = V_FIV(:);
            % RMS velocities
            urms_DNS_n=rms(U);
            uprime_rms_DNS_n=rms(U-mean(U));
            vrms_DNS_n=rms(V);
            vprime_rms_DNS_n=rms(V-mean(V));
            urms_FIV_n=rms(U_FIV);
            uprime_rms_FIV_n=rms(U_FIV-mean(U_FIV));
            vrms_FIV_n=rms(V_FIV);
            vprime_rms_FIV_n=rms(V_FIV-mean(V_FIV));
            %~~~~~~~~~~~~~~~~~~~~~~~~~~
            %---1-D Energy Spectrum---|
            %~~~~~~~~~~~~~~~~~~~~~~~~~~
            win=length(U);
            % DNS spectrum
            [E11_DNS_n,E22_DNS_n,K1_DNS,R11_DNS_n,R22_DNS_n] = EK_spectra(U,V,f_DNS,win,windo);
            % FIV spectrum
            [E11_FIV_n,E22_FIV_n,K1_FIV,R11_FIV_n,R22_FIV_n] = EK_spectra(U_FIV(:,1),V_FIV(:,1),f_DNS,win,windo);
            %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            %---1-D Instantanious Energy Dissipation---|
            %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            Udiff=gradient(U,mpx);
            index=log10(Udiff(:).^2)==-Inf;
            e0_DNS_E11_n=15*nu*10.^mean(log10(Udiff(index==0).^2));
            
            Vdiff=gradient(V,mpx);
            index=log10(Vdiff(:).^2)==-Inf;
            e0_DNS_E22_n=15/2*nu*10.^mean(log10(Vdiff(index==0).^2));
            
            Udiff=gradient(U_FIV,s{2,2}*mpx);
            index=log10(Udiff(:).^2)==-Inf;
            e0_FIV_E11_n=15*nu*10.^mean(log10(Udiff(index==0).^2));
        
            Vdiff=gradient(V_FIV,s{2,2}*mpx);
            index=log10(Vdiff(:).^2)==-Inf;
            e0_FIV_E22_n=15/2*nu*10.^mean(log10(Vdiff(index==0).^2));
            %~~~~~~~~~~~~~~~~~~~~~~~~~~
            %--------Data adding-------|
            %~~~~~~~~~~~~~~~~~~~~~~~~~~
            E11_DNS=E11_DNS+E11_DNS_n;
            E22_DNS=E22_DNS+E22_DNS_n;
            e0_DNS_E11=e0_DNS_E11+e0_DNS_E11_n;
            e0_DNS_E22=e0_DNS_E22+e0_DNS_E22_n;
            R11_DNS=R11_DNS+R11_DNS_n;
            R22_DNS=R22_DNS+R22_DNS_n;
            urms_DNS=urms_DNS+urms_DNS_n;
            uprime_rms_DNS=uprime_rms_DNS+uprime_rms_DNS_n;
            vrms_DNS=vrms_DNS+vrms_DNS_n;
            vprime_rms_DNS=vprime_rms_DNS+vprime_rms_DNS_n;
        
            E11_FIV=E11_FIV+E11_FIV_n;
            E22_FIV=E22_FIV+E22_FIV_n;
            e0_FIV_E11=e0_FIV_E11+e0_FIV_E11_n;
            e0_FIV_E22=e0_FIV_E22+e0_FIV_E22_n;
            R11_FIV=R11_FIV+R11_FIV_n;
            R22_FIV=R22_FIV+R22_FIV_n;
            urms_FIV=urms_FIV+urms_FIV_n;
            uprime_rms_FIV=uprime_rms_FIV+uprime_rms_FIV_n;
            vrms_FIV=vrms_FIV+vrms_FIV_n;
            vprime_rms_FIV=vprime_rms_FIV+vprime_rms_FIV_n;
        end
    end
end
% Averaging
E11_DNS=E11_DNS/count;
E22_DNS=E22_DNS/count;
e0_DNS_E11=e0_DNS_E11/count;
e0_DNS_E22=e0_DNS_E22/count;
R11_DNS=R11_DNS/count;
R22_DNS=R22_DNS/count;
urms_DNS=urms_DNS/count;
uprime_rms_DNS=uprime_rms_DNS/count;
vrms_DNS=vrms_DNS/count;
vprime_rms_DNS=vprime_rms_DNS/count;

E11_FIV=(E11_FIV/count);
E22_FIV=(E22_FIV/count);
e0_FIV_E11=e0_FIV_E11/count;
e0_FIV_E22=e0_FIV_E22/count;
R11_FIV=R11_FIV/count;
R22_FIV=R22_FIV/count;
urms_FIV=urms_FIV/count;
uprime_rms_FIV=uprime_rms_FIV/count;
vrms_FIV=vrms_FIV/count;
vprime_rms_FIV=vprime_rms_FIV/count;

%3-D Energy Spectrum
K1_DNS=K1_DNS(2:end);
E11_DNS=E11_DNS(2:end);
E22_DNS=E22_DNS(2:end);

K1_FIV=K1_FIV(2:end);
E11_FIV=E11_FIV(2:end);
E22_FIV=E22_FIV(2:end);

K1_DNS_log=log10(K1_DNS);
ET_log=log10(1/2*E11_DNS+E22_DNS);
[pT,~,muT]=polyfit(K1_DNS_log,ET_log,DNS_order);
ET_poly=polyval(pT,K1_DNS_log,[],muT);

kT=polyder(pT);
dET=polyval(kT,K1_DNS_log,[],muT);
EK_DNS=-K1_DNS_log.*(10.^ET_poly)*log(10).*dET;

K1_FIV_log=log10(K1_FIV);
ET_log=log10(1/2*E11_FIV+E22_FIV);
[pT,~,muT]=polyfit(K1_FIV_log,ET_log,FIV_order);
ET_poly=polyval(pT,K1_FIV_log,[],muT);

kT=polyder(pT);
dET=polyval(kT,K1_FIV_log,[],muT);
EK_FIV=-K1_FIV_log.*(10.^ET_poly)*log(10).*dET;


dK_DNS=K1_DNS(2)-K1_DNS(1);
dK_FIV=K1_FIV(2)-K1_FIV(1);


%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%----integral length scale----|
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% DNS
fx_DNS=R11_DNS/R11_DNS(1); % f(x)
gx_DNS=R22_DNS/R22_DNS(1); % g(x)
if overspace==1
    L11_DNS=sum(fx_DNS)*mpx;
    L22_DNS=sum(gx_DNS)*mpx;
else
    L11_DNS=sum(fx_DNS)*frsec;
    L22_DNS=sum(gx_DNS)*frsec;
end

% FIV
fx_FIV=R11_FIV/R11_FIV(1); % f(x)
gx_FIV=R22_FIV/R22_FIV(1); % g(x)
if overspace==1
    L11_FIV=sum(fx_FIV)*s{2,2}*mpx;
    L22_FIV=sum(gx_FIV)*s{2,2}*mpx;
else
    L11_FIV=sum(fx_FIV)*frsec;
    L22_FIV=sum(gx_FIV)*frsec;
end
%~~~~~~~~~~~~~~~~~~~~~~~~~
%--------TKED rate-------|
%~~~~~~~~~~~~~~~~~~~~~~~~~
epsilon_DNS_EK=2*nu*sum(K1_DNS.^2.*EK_DNS)*dK_DNS;
epsilon_DNS_E11=15*nu*sum(K1_DNS.^2.*E11_DNS)*dK_DNS;
epsilon_DNS_E22=15/2*nu*sum(K1_DNS.^2.*E22_DNS)*dK_DNS;

epsilon_FIV_EK=2*nu*sum(K1_FIV.^2.*EK_FIV)*dK_FIV;
epsilon_FIV_E11=15*nu*sum(K1_FIV.^2.*E11_FIV)*dK_FIV;
epsilon_FIV_E22=15/2*nu*sum(K1_FIV.^2.*E22_FIV)*dK_FIV;
%~~~~~~~~~~~~~~~~~~~~~~~~~
%-------Lengthscales------|
%~~~~~~~~~~~~~~~~~~~~~~~~~
%------DNS------
eta_DNS_E11=(nu^3/epsilon_DNS_E11)^(1/4);
eta_DNS_E22=(nu^3/epsilon_DNS_E22)^(1/4);
eta_DNS_EK=(nu^3/epsilon_DNS_EK)^(1/4);

tau_DNS_E11=(nu/epsilon_DNS_E11)^(1/2);
tau_DNS_E22=(nu/epsilon_DNS_E22)^(1/2);
tau_DNS_EK=(nu/epsilon_DNS_EK)^(1/2);

upsilon_DNS_E11=(nu*epsilon_DNS_E11)^(1/4);
upsilon_DNS_E22=(nu*epsilon_DNS_E22)^(1/4);
upsilon_DNS_EK=(nu*epsilon_DNS_EK)^(1/4);

% turbulent & Taylor Reynolds number
TKE_DNS=sum(E11_DNS)*dK_DNS;
uprime_rms_DNS_E11 = sqrt(2/3*TKE_DNS);
L_DNS_E11=TKE_DNS^(3/2)/epsilon_DNS_E11;
Re_L_DNS_E11=TKE_DNS^(1/2)*L_DNS_E11/nu;
% Re_T_DNS_E11=uprime_rms_DNS_E11*L11_DNS/nu
L_f_DNS=sqrt(10)*eta_DNS_E11^(2/3)*L_DNS_E11^(1/3);
R_lambda_DNS_E11=uprime_rms_DNS_E11*L_f_DNS/nu;

f2p=gradient(gradient(fx_DNS,mpx),mpx);
L_f_DNS_2p=((-1/2)*f2p(1))^(-1/2)
uprime=sqrt((epsilon_DNS_E11*(L_f_DNS_2p/sqrt(2))^2)/(15*nu));
R_lambda_DNS_E11_2p=uprime*L_f_DNS_2p/nu

TKE_DNS=sum(E22_DNS)*dK_DNS;
uprime_rms_DNS;
uprime_rms_DNS_E22 = sqrt(2/3*TKE_DNS);
L_DNS_E22=TKE_DNS^(3/2)/epsilon_DNS_E22;
Re_L_DNS_E22=TKE_DNS^(1/2)*L_DNS_E22/nu;
% Re_T_DNS_E22=uprime_rms_DNS_E22*L22_DNS/nu
L_g_DNS=sqrt(10)*eta_DNS_E22^(2/3)*L_DNS_E22^(1/3);
R_lambda_DNS_E22=uprime_rms_DNS_E22*L_g_DNS/nu;

g2p=gradient(gradient(gx_DNS,mpx),mpx);
L_g_DNS_2p=((-1/2)*g2p(1))^(-1/2)
uprime=sqrt((epsilon_DNS_E22*L_g_DNS_2p^2)/(15*nu));
R_lambda_DNS_E22_2p=uprime*L_g_DNS_2p/nu

TKE_DNS=sum(EK_DNS)*dK_DNS;
uprime_rms_DNS;
uprime_rms_DNS_EK = sqrt(2/3*TKE_DNS);
L_DNS_EK=TKE_DNS^(3/2)/epsilon_DNS_EK;
Re_L_DNS_EK=TKE_DNS^(1/2)*L_DNS_EK/nu;
L_EK_DNS=sqrt(10)*eta_DNS_EK^(2/3)*L_DNS_EK^(1/3);
R_lambda_DNS_EK=uprime_rms_DNS_EK*L_EK_DNS/nu;


%------FIV------
eta_FIV_E11=(nu^3/epsilon_FIV_E11)^(1/4);
eta_FIV_E22=(nu^3/epsilon_FIV_E22)^(1/4);
eta_FIV_EK=(nu^3/epsilon_FIV_EK)^(1/4);

tau_FIV_E11=(nu/epsilon_FIV_E11)^(1/2);
tau_FIV_E22=(nu/epsilon_FIV_E22)^(1/2);
tau_FIV_EK=(nu/epsilon_FIV_EK)^(1/2);

upsilon_FIV_E11=(nu*epsilon_FIV_E11)^(1/4);
upsilon_FIV_E22=(nu*epsilon_FIV_E22)^(1/4);
upsilon_FIV_EK=(nu*epsilon_FIV_EK)^(1/4);

% turbulent & Taylor Reynolds number
TKE_FIV=sum(E11_FIV)*dK_FIV;
uprime_rms_FIV;
uprime_rms_FIV_E11 = sqrt(2/3*TKE_FIV);
L_FIV_E11=TKE_FIV^(3/2)/epsilon_FIV_E11;
Re_L_FIV_E11=TKE_FIV^(1/2)*L_FIV_E11/nu;
% Re_T_FIV_E11=uprime_rms_FIV_E11*L11_FIV/nu
L_f_FIV=sqrt(10)*eta_FIV_E11^(2/3)*L_FIV_E11^(1/3);
R_lambda_FIV_E11=uprime_rms_FIV_E11*L_f_FIV/nu;

f2p=gradient(gradient(fx_FIV,mpx),mpx);
L_f_FIV_2p=abs(((-1/2)*f2p(1))^(-1/2))
uprime=sqrt((epsilon_FIV_E11*(L_f_FIV_2p/sqrt(2))^2)/(15*nu));
R_lambda_FIV_E11_2p=uprime*L_f_FIV_2p/nu

TKE_FIV=sum(E22_FIV)*dK_FIV;
uprime_rms_FIV;
uprime_rms_FIV_E22 = sqrt(2/3*TKE_FIV);
L_FIV_E22=TKE_FIV^(3/2)/epsilon_FIV_E22;
Re_L_FIV_E22=TKE_FIV^(1/2)*L_FIV_E22/nu;
% Re_T_FIV_E22=uprime_rms_FIV_E22*L22_FIV/nu
L_g_FIV=sqrt(10)*eta_FIV_E22^(2/3)*L_FIV_E22^(1/3);
R_lambda_FIV_E22=uprime_rms_FIV_E22*L_g_FIV/nu;

g2p=gradient(gradient(gx_FIV,mpx),mpx);
L_g_FIV_2p=abs(((-1/2)*g2p(1))^(-1/2))
uprime=sqrt((epsilon_FIV_E22*L_g_FIV_2p^2)/(15*nu));
R_lambda_FIV_E22_2p=uprime*L_g_FIV_2p/nu

TKE_FIV=sum(EK_FIV)*dK_FIV;
uprime_rms_FIV;
uprime_rms_FIV_EK = sqrt(2/3*TKE_FIV);
L_FIV_EK=TKE_FIV^(3/2)/epsilon_FIV_EK;
Re_L_FIV_EK=TKE_FIV^(1/2)*L_FIV_EK/nu;
L_EK_FIV=sqrt(10)*eta_FIV_EK^(2/3)*L_FIV_EK^(1/3);
R_lambda_FIV_EK=uprime_rms_FIV_EK*L_EK_FIV/nu;


Eii={'E_{11}';'E_{22}';'E(K)'};
eta=[eta_DNS_E11;eta_DNS_E22;eta_DNS_EK];
tau=[tau_DNS_E11;tau_DNS_E22;tau_DNS_EK];
upsilon=[upsilon_DNS_E11;upsilon_DNS_E22;upsilon_DNS_EK];
lambda=[L_f_DNS;L_g_DNS;L_EK_DNS];
R_lambda=[R_lambda_DNS_E11;R_lambda_DNS_E22;R_lambda_DNS_EK];
L=[L_DNS_E11;L_DNS_E22;L_DNS_EK];
R_L=[Re_L_DNS_E11;Re_L_DNS_E22;Re_L_DNS_EK];
epsilon=[epsilon_DNS_E11;epsilon_DNS_E22;epsilon_DNS_EK];
e0_DNS_EK=1.931541074851625e-06;
e0=[e0_DNS_E11;e0_DNS_E22;e0_DNS_EK];
DNS_scales=table(eta,tau,upsilon,lambda,R_lambda,L,R_L,epsilon,e0,'RowName',Eii)

Eii={'E_{11}';'E_{22}';'E(K)'};
eta=[eta_FIV_E11;eta_FIV_E22;eta_FIV_EK];
tau=[tau_FIV_E11;tau_FIV_E22;tau_FIV_EK];
upsilon=[upsilon_FIV_E11;upsilon_FIV_E22;upsilon_FIV_EK];
lambda=[L_f_FIV;L_g_FIV;L_EK_FIV];
R_lambda=[R_lambda_FIV_E11;R_lambda_FIV_E22;R_lambda_FIV_EK];
L=[L_FIV_E11;L_FIV_E22;L_FIV_EK];
R_L=[Re_L_FIV_E11;Re_L_FIV_E22;Re_L_FIV_EK];
epsilon=[epsilon_FIV_E11;epsilon_FIV_E22;epsilon_FIV_EK];
e0_FIV_EK=0;
e0=[e0_FIV_E11;e0_FIV_E22;e0_FIV_EK];
FIV_scales=table(eta,tau,upsilon,lambda,R_lambda,epsilon,e0,'RowName',Eii)

% Model Spectrum
K1_DNS=K1_DNS(1:365);
E11_DNS=E11_DNS(1:365);
E22_DNS=E22_DNS(1:365);
EK_DNS=EK_DNS(1:365);
%~~~~~~~~~~~~~~~~~~~~~~~~~
%--------Ideal Spec-------|
%~~~~~~~~~~~~~~~~~~~~~~~~~
% TKED=1.931541074851625e-06;
K_eta_DNS_E11=K1_DNS*eta_DNS_E11;
K_eta_DNS_E22=K1_DNS*eta_DNS_E22;
K_eta_DNS_EK=K1_DNS*eta_DNS_EK;
K_eta_FIV_E11=K1_FIV*eta_FIV_E11;
K_eta_FIV_E22=K1_FIV*eta_FIV_E22;
K_eta_FIV_EK=K1_FIV*eta_FIV_EK;

%--------Kolomgorov Model---------
C11=18/55;
C22=24/55;
K0 = 1.8333;
E11_m_kol=C11*K0*epsilon_DNS_E11^(2/3)*K1_DNS.^(-5/3);
E22_m_kol=C22*K0*epsilon_DNS_E22^(2/3)*K1_DNS.^(-5/3);
EK_m_kol =K0*epsilon_DNS_EK^(2/3)*K1_DNS.^(-5/3);

%--------Exponential Model---------
B=8.4;
a=5.1;
zeta_DNS_E11=a*eta_DNS_E11*K1_DNS;
zeta_DNS_E22=a*eta_DNS_E22*K1_DNS;
E11_m_exp=-(1/2*B*exp(-zeta_DNS_E11).*(1-zeta_DNS_E11)+1/2*B*expint(zeta_DNS_E11).*(zeta_DNS_E11.^2-2));
E22_m_exp=-(1/4*B*exp(-zeta_DNS_E22).*(zeta_DNS_E22-1)-1/4*B*expint(zeta_DNS_E22).*(zeta_DNS_E22.^2+2));
EK_m_exp=B*exp(-a*K_eta_DNS_EK);
%-----------------------------------
%~~~~~~~Non-Dimential Spectra~~~~~~~
%-----------------------------------
xlimit=[K1_DNS(1),K1_DNS(end)];
ylimit=[1e-15 1e-5];
xlimit=xlimit*eta_DNS_E11;
ylimit=ylimit/(epsilon_DNS_E11*nu^5)^(1/4);
ND_Spec=figure;
ND_Spec.Position(3:4) = [1600 1350]; 
tiledlayout(6,3)
nexttile([2,1])
%----DNS----
loglog(K_eta_DNS_E11,abs(E11_DNS)/(epsilon_DNS_E11*nu^5)^(1/4),'Color','black',LineWidth=2)
hold on
%----FIV----
loglog(K_eta_FIV_E11,abs(E11_FIV)/(epsilon_FIV_E11*nu^5)^(1/4),'Color',[0.00,0.45,0.74],LineWidth=2)
%---Model----
loglog(K_eta_DNS_E11(1:60),E11_m_kol(1:60)/(epsilon_DNS_E11*nu^5)^(1/4),LineWidth=2,LineStyle="--")
loglog(K_eta_DNS_E11(21:end),E11_m_exp(21:end),LineWidth=2,LineStyle="-.")

xlabel('$K_1\eta$','Interpreter','latex','FontSize',24);
ylabel('$E_{11}/(\epsilon\nu^5)^{1/4}$','Interpreter','latex','FontSize',24);
xlim(xlimit)
ylim(ylimit)
lab1 = '$E_{11}(K_1)=\frac{18}{55}K_0\epsilon^{2/3}K_1^{-5/3}$';
lab2 = '$E_{11}(K_1)=\frac{1}{2}B\cdot e^{-\xi}[\xi-1]-\frac{1}{2}B\cdot E_{in}(\xi)[\xi^2-2]$';
legend('DNS','FIV',lab1,lab2,'FontSize',11,'Interpreter','latex','Location','southwest')
legend('boxoff')
text(1*10^(-2.2),1e-1,'A','FontSize',24)
hold off

nexttile([2,1])
%----DNS----
loglog(K_eta_DNS_E22,abs(E22_DNS)/(epsilon_DNS_E22*nu^5)^(1/4),'Color','black',LineWidth=2)
hold on
%----FIV----
loglog(K_eta_FIV_E22,abs(E22_FIV)/(epsilon_FIV_E22*nu^5)^(1/4),'Color',[0.00,0.45,0.74],LineWidth=2)
%---Model----
loglog(K_eta_DNS_E22(1:60),E22_m_kol(1:60)/(epsilon_DNS_E11*nu^5)^(1/4),LineWidth=2,LineStyle="--")
loglog(K_eta_DNS_E22(21:end),E22_m_exp(21:end),LineWidth=2,LineStyle="-.")

xlabel('$K_1\eta$','Interpreter','latex','FontSize',24);
ylabel('$E_{22}/(\epsilon\nu^5)^{1/4}$','Interpreter','latex','FontSize',24);
xlim(xlimit)
ylim(ylimit)
lab1 = '$E_{22}(K_1)=\frac{24}{55}K_0\epsilon^{2/3}K_1^{-5/3}$';
lab2 = '$E_{22}(K_1)=\frac{1}{4}B\cdot e^{-\xi}[1-\xi]+\frac{1}{4}B\cdot E_{in}(\xi)[\xi^2+2]$';
legend('DNS','FIV',lab1,lab2,'FontSize',11,'Interpreter','latex','Location','southwest')
legend('boxoff')
text(1*10^(-2.2),1e-1,'B','FontSize',24)
hold off

nexttile([2,1])
%----DNS----
loglog(K_eta_DNS_EK,abs(EK_DNS)/(epsilon_DNS_EK*nu^5)^(1/4),'Color','black',LineWidth=2)
hold on
%----FIV----
loglog(K_eta_FIV_EK,abs(EK_FIV)/(epsilon_FIV_EK*nu^5)^(1/4),'Color',[0.00,0.45,0.74],LineWidth=2)
%---Model----
loglog(K_eta_DNS_EK(1:60),EK_m_kol(1:60)/(epsilon_DNS_EK*nu^5)^(1/4),LineWidth=2,LineStyle="--")
loglog(K_eta_DNS_EK(21:end),EK_m_exp(21:end),LineWidth=2,LineStyle="-.")

xlabel('$K\eta$','Interpreter','latex','FontSize',24);
ylabel('$E(K)/(\epsilon\nu^5)^{1/4}$','Interpreter','latex','FontSize',24);
xlim(xlimit)
ylim(ylimit)
lab1 = '$E(K)=K_0\epsilon^{2/3}K_1^{-5/3}$';
lab2 = '$E(K)\sim B\cdot e^{-a\eta K}$';
legend('DNS','FIV',lab1,lab2,'FontSize',11,'Interpreter','latex','Location','southwest')
legend('boxoff')
text(1*10^(-2.2),1e-1,'C','FontSize',24)
hold off

%------------------------------------------------
%~~~~~~~Compensated inertial range spectra~~~~~~~
%------------------------------------------------
nexttile([2,1])
%----DNS----
semilogx(K_eta_DNS_E11,K1_DNS.^(5/3).*abs(E11_DNS)*(epsilon_DNS_E11)^(-2/3),'Color','black',LineWidth=2)
hold on
%----FIV----
semilogx(K_eta_FIV_E11,K1_FIV.^(5/3).*abs(E11_FIV)*(epsilon_FIV_E11)^(-2/3),'Color',[0.00,0.45,0.74],LineWidth=2)
semilogx([0.0815202,0.0815202],[0.1,1],'Color',[0.00,0.45,0.74],LineWidth=2,LineStyle=":")
%---Model----
semilogx(K_eta_DNS_E11(1:60),K1_DNS(1:60).^(5/3).*E11_m_kol(1:60)*(epsilon_DNS_E11)^(-2/3),'Color',[0.93,0.69,0.13],LineWidth=2,LineStyle="--")


xlabel('$K_1\eta$','Interpreter','latex','FontSize',24);
ylabel('$(K_1\eta)^{5/3}E_{11}/(\epsilon)^{2/3}$','Interpreter','latex','FontSize',24);
xlim(xlimit)
ylim([0 6])
lab2 = '$E_{11}(K_1)=\frac{18}{55}K_0\epsilon^{2/3}K_1^{-5/3}$';
lab1 = 'FIV resolution';
legend('DNS','FIV',lab1,lab2,'FontSize',11,'Interpreter','latex','Location','northwest')
legend('boxoff')
text(1*10^(-2.2),3,'D','FontSize',24)
hold off

nexttile([2,1])
%----DNS----
semilogx(K_eta_DNS_E22,K1_DNS.^(5/3).*abs(E22_DNS)*(epsilon_DNS_E22)^(-2/3),'Color','black',LineWidth=2)
hold on
%----FIV----
semilogx(K_eta_FIV_E22,K1_FIV.^(5/3).*abs(E22_FIV)*(epsilon_FIV_E22)^(-2/3),'Color',[0.00,0.45,0.74],LineWidth=2)
semilogx([0.125346,0.125346],[0.3,1.5],'Color',[0.00,0.45,0.74],LineWidth=2,LineStyle=":")
%---Model----
semilogx(K_eta_DNS_E22(1:60),K1_DNS(1:60).^(5/3).*E22_m_kol(1:60)*(epsilon_DNS_E22)^(-2/3),'Color',[0.93,0.69,0.13],LineWidth=2,LineStyle="--")


xlabel('$K_1\eta$','Interpreter','latex','FontSize',24);
ylabel('$(K_1\eta)^{5/3}E_{22}/(\epsilon)^{2/3}$','Interpreter','latex','FontSize',24);
xlim(xlimit)
ylim([0 6])
lab2 = '$E_{22}(K_1)=\frac{24}{55}K_0\epsilon^{2/3}K_1^{-5/3}$';
lab1 = 'FIV resolution';
legend('DNS','FIV',lab1,lab2,'FontSize',11,'Interpreter','latex','Location','northwest')
legend('boxoff')
text(1*10^(-2.2),3,'E','FontSize',24)
hold off

nexttile([2,1])
%----DNS----
semilogx(K_eta_DNS_EK,K1_DNS.^(5/3).*abs(EK_DNS)*(epsilon_DNS_EK)^(-2/3),'Color','black',LineWidth=2)
hold on
%----FIV----
semilogx(K_eta_FIV_EK,K1_FIV.^(5/3).*abs(EK_FIV)*(epsilon_FIV_EK)^(-2/3),'Color',[0.00,0.45,0.74],LineWidth=2)

K_win_eta=0.130129;
semilogx([K_win_eta,K_win_eta],[2,5],'Color',[0.00,0.45,0.74],LineWidth=2,LineStyle=":")
%---Model----
semilogx(K_eta_DNS_EK(1:60),K1_DNS(1:60).^(5/3).*EK_m_kol(1:60)*(epsilon_DNS_EK)^(-2/3),'Color',[0.93,0.69,0.13],LineWidth=2,LineStyle="--")

xlabel('$K\eta$','Interpreter','latex','FontSize',24);
ylabel('$(K\eta)^{5/3}E(K)/(\epsilon)^{2/3}$','Interpreter','latex','FontSize',24);
xlim(xlimit)
ylim([0 12])
lab2 = '$E(K)=K_0\epsilon^{2/3}K_1^{-5/3}$';
lab1 = 'FIV resolution';
legend('DNS','FIV',lab1,lab2,'FontSize',11,'Interpreter','latex','Location','northwest')
legend('boxoff')
text(1*10^(-2.2),6,'F','FontSize',24)
hold off

%---------------------------------------------------
%~~~~~~~~~~~~~~~~~~~K2 spectra~~~~~~~~~~~~~~~~~~~~~~
%---------------------------------------------------
nexttile([2,1])
%----DNS----
semilogx(K_eta_DNS_E11,K1_DNS.^(2).*abs(E11_DNS)*(epsilon_DNS_E11)^(-2/3),'Color','black',LineWidth=2)
hold on
%----FIV----
semilogx(K_eta_FIV_E11,K1_FIV.^(2).*abs(E11_FIV)*(epsilon_FIV_E11)^(-2/3),'Color',[0.00,0.45,0.74],LineWidth=2)
semilogx([0.0815202,0.0815202],[1,4],'Color',[0.00,0.45,0.74],LineWidth=2,LineStyle=":")
%---Model----
semilogx(K_eta_DNS_E11(1:60),K1_DNS(1:60).^(2).*E11_m_kol(1:60)*(epsilon_DNS_E11)^(-2/3),'Color',[0.93,0.69,0.13],LineWidth=2,LineStyle="--")


xlabel('$K_1\eta$','Interpreter','latex','FontSize',12);
ylabel('$(K_1\eta)^{2}E_{11}/(\epsilon)^{2/3}$','Interpreter','latex','FontSize',12);
xlim(xlimit)
ylim([0 18])
lab2 = '$E_{11}(K_1)=\frac{18}{55}K_0\epsilon^{2/3}K_1^{-5/3}$';
lab1 = 'FIV resolution';
legend('DNS','FIV',lab1,lab2,'FontSize',9,'Interpreter','latex','Location','northwest')
legend('boxoff')
text(1*10^(-2.2),9,'G','FontSize',12)
hold off

nexttile([2,1])
%----DNS----
semilogx(K_eta_DNS_E22,K1_DNS.^(2).*abs(E22_DNS)*(epsilon_DNS_E22)^(-2/3),'Color','black',LineWidth=2)
hold on
%----FIV----
semilogx(K_eta_FIV_E22,K1_FIV.^(2).*abs(E22_FIV)*(epsilon_FIV_E22)^(-2/3),'Color',[0.00,0.45,0.74],LineWidth=2)
semilogx([0.125346,0.125346],[2,6],'Color',[0.00,0.45,0.74],LineWidth=2,LineStyle=":")
%---Model----
semilogx(K_eta_DNS_E22(1:60),K1_DNS(1:60).^(2).*E22_m_kol(1:60)*(epsilon_DNS_E22)^(-2/3),'Color',[0.93,0.69,0.13],LineWidth=2,LineStyle="--")


xlabel('$K_1\eta$','Interpreter','latex','FontSize',12);
ylabel('$(K_1\eta)^{2}E_{22}/(\epsilon)^{2/3}$','Interpreter','latex','FontSize',12);
xlim(xlimit)
ylim([0 18])
lab2 = '$E_{22}(K_1)=\frac{24}{55}K_0\epsilon^{2/3}K_1^{-5/3}$';
lab1 = 'FIV resolution';
legend('DNS','FIV',lab1,lab2,'FontSize',9,'Interpreter','latex','Location','northwest')
legend('boxoff')
text(1*10^(-2.2),9,'H','FontSize',12)
hold off

nexttile([2,1])
%----DNS----
semilogx(K_eta_DNS_EK,K1_DNS.^(2).*abs(EK_DNS)*(epsilon_DNS_EK)^(-2/3),'Color','black',LineWidth=2)
hold on
%----FIV----
semilogx(K_eta_FIV_EK,K1_FIV.^(2).*abs(EK_FIV)*(epsilon_FIV_EK)^(-2/3),'Color',[0.00,0.45,0.74],LineWidth=2)

K_win_eta=0.130129;
semilogx([K_win_eta,K_win_eta],[16,22],'Color',[0.00,0.45,0.74],LineWidth=2,LineStyle=":")
%---Model----
semilogx(K_eta_DNS_EK(1:60),K1_DNS(1:60).^(2).*EK_m_kol(1:60)*(epsilon_DNS_EK)^(-2/3),'Color',[0.93,0.69,0.13],LineWidth=2,LineStyle="--")

xlabel('$K\eta$','Interpreter','latex','FontSize',12);
ylabel('$(K\eta)^{2}E(K)/(\epsilon)^{2/3}$','Interpreter','latex','FontSize',12);
xlim(xlimit)
ylim([0 40])
lab2 = '$E(K)=K_0\epsilon^{2/3}K_1^{-5/3}$';
lab1 = 'FIV resolution';
legend('DNS','FIV',lab1,lab2,'FontSize',9,'Interpreter','latex','Location','northwest')
legend('boxoff')
text(1*10^(-2.2),20,'I','FontSize',12)
hold off