%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Landau equation (singular kernel, Coulomb case) 
% dt f(v,t) = Q(f,f)(v,t)
% Approach CST2 based on conservative form (mass conservation) 
% Q(f,f)(v,t) = dv . QC(f,f)(v,t)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Space discretisation by Fourier spectral method 
% Time discretisation by explicit Runge-Kutta method 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Test = 3 (Figure 10 with final time T = 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all
clear all 
format long e 
clc
pause off

disp('Time integration by explicit Runge-Kutta method.')

PrecomputationYes = 0; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if PrecomputationYes == 0
    % Problem data 
    disp('See Carrillo et al. (Example 4, Figure 8)')
    Param.vAbsSquare = inline('v{1}.^2 + v{2}.^2 + v{3}.^2','v');
    Param.S = 10;
    Param.sigma = 0.3;
    Param.f0 ...
    = inline(['1/Param.S^2*exp(-Param.S/Param.sigma^2' ...
      '*(Param.vAbsSquare(v).^(1/2) - Param.sigma).^2)'],'v','Param');
    Param.CKernel = 1/(4*pi);
    Param.epsilon = 0;
    Param.phi ...
    = inline(['Param.CKernel*1./(Param.epsilon ' ...
      '+ Param.vAbsSquare(v)).^(3/2)'],'v','Param');
    % Dimension (velocity) 
    Param.d = 3;         
    % Truncated velocity domain 
    Param.a1 = - 10*ones(Param.d,1);
    Param.a2 = 10*ones(Param.d,1);
    for i = 1:Param.d
        Param.b{i}(1) = Param.a1(i);
        Param.b{i}(2) = Param.a2(i);
        Param.db(i,1) = Param.b{i}(2)-Param.b{i}(1);
        Param.bPlus(i,1) = Param.b{i}(1) + Param.b{i}(2);
    end
    % Velocity grid and neighbourhood of origin 
    M = 64;
    Width = 2;
    Param.M = M*ones(Param.d,1);
    Param.Width = Width*ones(Param.d,1);
    % Time grid  
    Param.t0 = 0;
    Param.T = Param.t0 + 1;
    Param.N = round(2000*(Param.T-Param.t0));            
    Param.dt = (Param.T - Param.t0)/Param.N;
    % Runge-Kutta method of order 4 
    Param.RKM = 4;
    % Plot 
    Param = SetBasicsPlot(Param);
    % Precomputation (Grid, Integrals, Kernel)
    disp('Precomputation')
    tic
    Param = Precomputation_Fourier(Param);
    Param = MyBasicIntegrals(Param);
    Param = Precomputation_Kernel_ApproachCST2_Simplified3d(Param);
    % Param = MyIndices(Param);
    TimePrecomputation = toc
    save D_d3_Precomputation.mat
    pause 
end 

if PrecomputationYes == 1
    load D_d3_Precomputation.mat
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time integration 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
n = 0;
f = Param.f0(Param.v,Param);
fs = Fourier_Real2Spectral(f,Param); 
Param = Quantities(n,f,fs,Param);
n 
Mass = Param.MassSpectral(n+1)
Momentum = Param.MomentumSpectral(n+1)
Energy = Param.EnergySpectral(n+1)
Entropy = Param.Entropy(n+1)
pause 

for n = 1:Param.N
    if n == 1 
        tic
    end 
    [f,fs] = ExplicitRKM(Param.RKM,f,fs,Param);
    Param = Quantities(n,f,fs,Param);            
    if n == 1 
        TimeComputationPerTimeStep = toc;
    end 
    if mod(n,Param.N/100) == 1 
        n 
        Mass = Param.MassSpectral(n+1)
        Momentum = Param.MomentumSpectral(n+1)
        Energy = Param.EnergySpectral(n+1)
        Entropy = Param.Entropy(n+1)
    end
    if n == Param.N/2
        Param.THalf = Param.Times(n+1);
        Param.fHalf = f;
    end    
end

save D_d3.mat
MyPlotPrint_Simplified3d(n,f,Param)        

TimePrecomputation
TimeComputationPerTimeStep

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Param = SetBasicsPlot(Param)

Param.FontSize = 15; 
Param.LineWidth = 1.5;
Param.ColorMyRed = [82.75,43.14,43.92]/100;
Param.ColorMyGrey = [85.88,85.88,85.88]/100;

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialisation (Fourier spectral method)  
% Grid (Trapezoidal quadrature nodes) 
% Fourier basis functions Fm(x) = C exp(i pi m/c (x+c))
% Associated eigenvalues for partial derivatives dx1,dx2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Param = Precomputation_Fourier(Param)

if max(mod(Param.M,2)) == 1 | max(Param.M <= 0) == 1
    error('Parameter M should comprise even positive integers.')
end

K = Param.M;
Param.dv = 1;
for i = 1:Param.d
    Param.vColumn{i} = linspace(Param.b{i}(1),Param.b{i}(2),K(i)+1)';
    Param.vColumn{i} = Param.vColumn{i}(1:end-1);
    w(i) = Param.vColumn{i}(2)-Param.vColumn{i}(1);
    mu{i} = 2*pi*1i/Param.db(i)*[-Param.M(i)/2:Param.M(i)/2-1]';
end
[Param.v{1},Param.v{2},Param.v{3}] ...
= ndgrid(Param.vColumn{1},Param.vColumn{2},Param.vColumn{3});
[Param.Mu{1},Param.Mu{2},Param.Mu{3}] = ndgrid(mu{1},mu{2},mu{3});

% Constants in spectral transformations and scalar product
Param.Constr2s = prod(sqrt(Param.db)./Param.M);
Param.Consts2r = prod(Param.M./sqrt(Param.db));
Param.Constscp = prod(w);

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fourier transform
% Inverse Fourier transform 
% !!! Real part of result !!! 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function fs = Fourier_Real2Spectral(f,Param)

fs = Param.Constr2s*fftshift(fftn(f));

end

function f = Fourier_Spectral2Real(fs,Param)
 
f = Param.Consts2r*real(ifftn(ifftshift(fs)));

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Basic integrals  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Param = MyBasicIntegrals(Param)

for i = 1:Param.d
    Param.I0{i} = sqrt(Param.db(i));
    mu = 2*pi*1i/Param.db(i)*[-Param.M(i)+1:Param.M(i)-1];
    I1_Aux{i} = sqrt(Param.db(i))./mu;
    I1_Aux{i}(Param.M(i)) ...
    = 1/2*(Param.b{i}(2)^2-Param.b{i}(1)^2)/sqrt(Param.db(i));
    I2_Aux{i} ...
    = ((Param.b{i}(2)^2-Param.b{i}(1)^2)*mu - 2*Param.db(i)) ...
      ./(sqrt(Param.db(i))*mu.^2);
    I2_Aux{i}(Param.M(i)) ...
    = 1/3*(Param.b{i}(2)^3-Param.b{i}(1)^3)/sqrt(Param.db(i));
    j = 0;
    Param.I1{i} = I1_Aux{i}(Param.M(i)+j)*diag(ones(Param.M(i)-j,1),j);
    Param.I2{i} = I2_Aux{i}(Param.M(i)+j)*diag(ones(Param.M(i)-j,1),j);
    for j = 1:Param.M(i)-1
        Param.I1{i} ...
        = Param.I1{i} ...
          + I1_Aux{i}(Param.M(i)+j)*diag(ones(Param.M(i)-j,1),j);
        Param.I1{i} ...
        = Param.I1{i} ...
          + I1_Aux{i}(Param.M(i)-j)*diag(ones(Param.M(i)-j,1),-j);
        Param.I2{i} ...
        = Param.I2{i} ...
          + I2_Aux{i}(Param.M(i)+j)*diag(ones(Param.M(i)-j,1),j);
        Param.I2{i} ...
        = Param.I2{i} ...
          + I2_Aux{i}(Param.M(i)-j)*diag(ones(Param.M(i)-j,1),-j);
    end
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Precomputations involving kernel  
% Evaluation of phi 
% Modification at isolated singularity (mean value)
% Cut-off large positive values of phi to define psi 
% Restriction of v and difference phi - psi to small domain
% In our paper: IWithoutFm_Psi{1,1} corresponds to I_{2,2} etc. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Param = Precomputation_Kernel_ApproachCST2_Simplified3d(Param)

% Evaluation of phi 
phi = Param.phi(Param.v,Param);

% Modification at isolated singularity 
[Ind1Inf,Ind2Inf,Ind3Inf] = find(phi == inf);
if length(Ind1Inf) > 0 
    phiAux = reshape(phi(Ind1Inf,:,:),Param.M(2),Param.M(3));
    [Ind2Inf,Ind3Inf] = find(phiAux == inf);
    disp('Modification at singularity')
    Indices = [Ind1Inf,Ind2Inf,Ind3Inf]
    phi(Ind1Inf,Ind2Inf,Ind3Inf) ...
    = 1/2*(phi(Ind1Inf-1,Ind2Inf,Ind3Inf) ...
      + phi(Ind1Inf+1,Ind2Inf,Ind3Inf));
end 

% Restriction of v to small neighbourhood of singularity 
for i = 1:Param.d
    MinM(i) = max(1,Param.M(i)/2 + 1 - Param.Width(i));
    MaxM(i) = min(Param.M(i),Param.M(i)/2 + 1 + Param.Width(i));
end 
for i = 1:Param.d 
    vRestricted{i} ...
    = Param.v{i}(MinM(1):MaxM(1),MinM(2):MaxM(2),MinM(3):MaxM(3));
end

% Cut-off to define psi 
psi = phi;
for m2 = MinM(2):MaxM(2)
    for m3 = MinM(3):MaxM(3)
        vFirst = Param.v{1}(MinM(1),m2,m3);
        psiFirst = psi(MinM(1),m2,m3);
        vLast = Param.v{1}(MaxM(1),m2,m3);
        psiLast = psi(MaxM(1),m2,m3);
        psi(MinM(1):MaxM(1),m2,m3) ...
        = psiFirst ...
          + (psiLast-psiFirst)/(vLast-vFirst) ...
          *(Param.v{1}(MinM(1):MaxM(1),m2,m3) - vFirst);
    end    
end

% Restriction of phi - psi to small neighbourhood of singularity
phiMinuspsiRestricted ...
= phi(MinM(1):MaxM(1),MinM(2):MaxM(2),MinM(3):MaxM(3)) ...
  - psi(MinM(1):MaxM(1),MinM(2):MaxM(2),MinM(3):MaxM(3));

% Auxiliary quantities for integrals involving psi 
psis = Fourier_Real2Spectral(psi,Param);
for i = 1:Param.d
    aux{i} ...
    = exp(-2*pi*1i/Param.db(i)*[-Param.M(i)/2:Param.M(i)/2-1]' ...
      *Param.b{i}(1));
end
[aux1,aux2,aux3] = ndgrid(aux{1},aux{2},aux{3});
MyGamma = aux1.*aux2.*aux3;
InnerSum = InnerSum3d('200',Param.I2{1},Param.I0{2},Param.I0{3},psis);
IWithoutFm_Psi{1,1} = MyGamma.*InnerSum;
InnerSum = InnerSum3d('110',Param.I1{1},Param.I1{2},Param.I0{3},psis);
IWithoutFm_Psi{1,2} = MyGamma.*InnerSum;
InnerSum = InnerSum3d('101',Param.I1{1},Param.I0{2},Param.I1{3},psis);
IWithoutFm_Psi{1,3} = MyGamma.*InnerSum;
InnerSum = InnerSum3d('020',Param.I0{1},Param.I2{2},Param.I0{3},psis);
IWithoutFm_Psi{2,2} = MyGamma.*InnerSum;
InnerSum = InnerSum3d('011',Param.I0{1},Param.I1{2},Param.I1{3},psis);
IWithoutFm_Psi{2,3} = MyGamma.*InnerSum;
InnerSum = InnerSum3d('002',Param.I0{1},Param.I0{2},Param.I2{3},psis);
IWithoutFm_Psi{3,3} = MyGamma.*InnerSum;

% Computation of integrals involving phi - psi by linear integral transform 
% u = v - w and quadrature 
for i = 1:Param.d
    E{i} = exp(-2*pi*1i*vRestricted{i}/Param.db(i)); 
end
dv(1) = vRestricted{1}(2,1,1) - vRestricted{1}(1,1,1);
dv(2) = vRestricted{2}(1,2,1) - vRestricted{2}(1,1,1);
dv(3) = vRestricted{3}(1,1,2) - vRestricted{3}(1,1,1);
dv123 = prod(dv);
for m1 = 1:Param.M(1)
    for m2 = 1:Param.M(2)
        for m3 = 1:Param.M(3)
            EAux ...
            = E{1}.^(m1-Param.M(1)/2-1) ...
              .*E{2}.^(m2-Param.M(2)/2-1) ...
              .*E{3}.^(m3-Param.M(3)/2-1);
            IWithoutFm_PhiMinusPsi{1,1}(m1,m2,m3) ...
            = dv123*sum(sum(sum(vRestricted{1}.^2 ...
              .*phiMinuspsiRestricted.*EAux)));
            IWithoutFm_PhiMinusPsi{1,2}(m1,m2,m3) ...
            = dv123*sum(sum(sum(vRestricted{1}.*vRestricted{2} ...
              .*phiMinuspsiRestricted.*EAux)));
            IWithoutFm_PhiMinusPsi{1,3}(m1,m2,m3) ...
            = dv123*sum(sum(sum(vRestricted{1}.*vRestricted{3} ...
              .*phiMinuspsiRestricted.*EAux)));
            IWithoutFm_PhiMinusPsi{2,2}(m1,m2,m3) ...
            = dv123*sum(sum(sum(vRestricted{2}.^2 ...
              .*phiMinuspsiRestricted.*EAux)));
            IWithoutFm_PhiMinusPsi{2,3}(m1,m2,m3) ...
            = dv123*sum(sum(sum(vRestricted{2}.*vRestricted{3} ...
              .*phiMinuspsiRestricted.*EAux)));
            IWithoutFm_PhiMinusPsi{3,3}(m1,m2,m3) ...
            = dv123*sum(sum(sum(vRestricted{3}.^2 ...
              .*phiMinuspsiRestricted.*EAux)));
        end
    end 
end 

% Sum of auxiliary quantities
Param.IWithoutFm{1,1} = IWithoutFm_Psi{1,1} + IWithoutFm_PhiMinusPsi{1,1};
Param.IWithoutFm{1,2} = IWithoutFm_Psi{1,2} + IWithoutFm_PhiMinusPsi{1,2};
Param.IWithoutFm{1,3} = IWithoutFm_Psi{1,3} + IWithoutFm_PhiMinusPsi{1,3};
Param.IWithoutFm{2,2} = IWithoutFm_Psi{2,2} + IWithoutFm_PhiMinusPsi{2,2};
Param.IWithoutFm{2,3} = IWithoutFm_Psi{2,3} + IWithoutFm_PhiMinusPsi{2,3};
Param.IWithoutFm{3,3} = IWithoutFm_Psi{3,3} + IWithoutFm_PhiMinusPsi{3,3};

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inner sums
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Value = InnerSum3d(Choice,I1,I2,I3,gs)

Value = zeros(size(gs));

if Choice == '100' | Choice == '200'
    Value = I2*I3*einsum(I1,gs,2,1);
end
if Choice == '010' | Choice == '020'
    Value = permute(gs,[2,1,3]);
    Value = I1*I3*einsum(I2,Value,2,1);
    Value = permute(Value,[2,1,3]);
end
if Choice == '001' | Choice == '002'
    Value = permute(gs,[3,2,1]);
    Value = I1*I2*einsum(I3,Value,2,1);
    Value = permute(Value,[3,2,1]);
end
if Choice == '110'
    Value = I3*einsum(I1,gs,2,1);
    Value = permute(Value,[2,1,3]);
    Value = einsum(I2,Value,2,1);
    Value = permute(Value,[2,1,3]);
end
if Choice == '101'
    Value = I2*einsum(I1,gs,2,1);
    Value = permute(Value,[3,2,1]);
    Value = einsum(I3,Value,2,1);
    Value = permute(Value,[3,2,1]);
end
if Choice == '011'
    Value = permute(gs,[2,1,3]);
    Value = I1*einsum(I2,Value,2,1);
    Value = permute(Value,[2,1,3]);
    Value = permute(Value,[3,2,1]);
    Value = einsum(I3,Value,2,1);
    Value = permute(Value,[3,2,1]);
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Landau operator 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Q,Qs,Param] = MyQ_Simplified3d(f,fs,Param)

for i = 1:Param.d 
    % Spectral coefficients of first derivatives of f 
    dfs.dvfs{i} = Param.Mu{i}.*fs;
    % First derivatives of f 
    df.dvf{i} = Fourier_Spectral2Real(dfs.dvfs{i},Param);
end 

[Q,Qs,Param] = MyQ_ApproachCST2_Simplified3d(f,fs,df,dfs,Param);     

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computation of integrals involving psi by linear integral
% transform u = v - w and approximation of v - Omega by Omega
% ! Assumption that domain v - Omega = Omega permits use of IFFT for
% summation (I_m(v) = Fm(v) J_m) !
% Computation of integrals involving phi - psi by linear
% integral transform, approximation of v - Omega by Omega, and quadrature
% Computation of corresponding terms in Landau operator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Q,Qs,Param] = MyQ_ApproachCST2_Simplified3d(f,fs,df,dfs,Param)

for i = 1:Param.d
    MyI{i,i} = Fourier_Spectral2Real(fs.*Param.IWithoutFm{i,i},Param);
    for j = 1:Param.d
        if j == i
        else
            MyMuI{i,i}{j} ...
                = Fourier_Spectral2Real ...
                  (dfs.dvfs{j}.*Param.IWithoutFm{i,i},Param);
        end    
    end
    if i == 1
        MyI{i,i+1} ...
        = Fourier_Spectral2Real(fs.*Param.IWithoutFm{i,i+1},Param);
        for j = [1,2]
            MyMuI{i,i+1}{j} ...
            = Fourier_Spectral2Real ...
              (dfs.dvfs{j}.*Param.IWithoutFm{i,i+1},Param);
        end
        MyI{i,i+2} ...
        = Fourier_Spectral2Real(fs.*Param.IWithoutFm{i,i+2},Param);
        for j = [1,3]
            MyMuI{i,i+2}{j} ...
                = Fourier_Spectral2Real ...
                (dfs.dvfs{j}.*Param.IWithoutFm{i,i+2},Param);
        end    
    end
    if i == 2
        MyI{i,i+1} ...
        = Fourier_Spectral2Real(fs.*Param.IWithoutFm{i,i+1},Param);
        for j = [2,3]
            MyMuI{i,i+1}{j} ...
            = Fourier_Spectral2Real ...
              (dfs.dvfs{j}.*Param.IWithoutFm{i,i+1},Param);
        end    
    end    
end
    
Q000{1} = - MyMuI{2,2}{1} - MyMuI{3,3}{1} + MyMuI{1,2}{2} + MyMuI{1,3}{3};
Q000{2} = - MyMuI{1,1}{2} - MyMuI{3,3}{2} + MyMuI{1,2}{1} + MyMuI{2,3}{3};
Q000{3} = - MyMuI{1,1}{3} - MyMuI{2,2}{3} + MyMuI{1,3}{1} + MyMuI{2,3}{2};
Q100{1} = MyI{2,2} + MyI{3,3};
Q010{2} = MyI{1,1} + MyI{3,3};
Q001{3} = MyI{1,1} + MyI{2,2};
Q010{1} = - MyI{1,2};
Q001{1} = - MyI{1,3};
Q100{2} = - MyI{1,2};
Q001{2} = - MyI{2,3};
Q100{3} = - MyI{1,3};
Q010{3} = - MyI{2,3};
for i = 1:Param.d 
    QC{i} = Q000{i}.*f + Q100{i}.*df.dvf{1} + Q010{i}.*df.dvf{2} ...
            + Q001{i}.*df.dvf{3};
end
for i = 1:Param.d
    Aux{i} = Fourier_Real2Spectral(QC{i},Param);
end
Qs = Param.Mu{1}.*Aux{1} + Param.Mu{2}.*Aux{2} + Param.Mu{3}.*Aux{3}; 
Q = Fourier_Spectral2Real(Qs,Param);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [f,fs] = ExplicitRKM(RKM,f,fs,Param)

dt = Param.dt;

% Explicit Runge-Kutta method (Order 1)
if RKM == 1
    [Q,Qs,Param] = MyQ_Simplified3d(f,fs,Param);
    f = f + dt*Q;
    fs = fs + dt*Qs;
end 

% Explicit Runge-Kutta method (Order 4)
if RKM == 4
    [QStage1,QsStage1,Param] = MyQ_Simplified3d(f,fs,Param);
    fAux = f + dt/2*QStage1;
    fsAux = fs + dt/2*QsStage1;
    [QStage2,QsStage2,Param] = MyQ_Simplified3d(fAux,fsAux,Param);
    fAux = f + dt/2*QStage2;
    fsAux = fs + dt/2*QsStage2;
    [QStage3,QsStage3,Param] = MyQ_Simplified3d(fAux,fsAux,Param);
    fAux = f + dt*QStage3;
    fsAux = fs + dt*QsStage3;
    [QStage4,QsStage4,Param] = MyQ_Simplified3d(fAux,fsAux,Param);
    f = f + dt/6*(QStage1 + 2*QStage2 + 2*QStage3 + QStage4);
    fs = fs + dt/6*(QsStage1 + 2*QsStage2 + 2*QsStage3 + QsStage4);
end 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mass, Momentum, Energy, Entropy  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Param = Quantities(n,f,fs,Param)

Param.Times(n+1,1) = Param.t0 + n*Param.dt;

% Auxiliaries (Quadrature)
dv(1) = Param.v{1}(2,1,1) - Param.v{1}(1,1,1);
dv(2) = Param.v{2}(1,2,1) - Param.v{2}(1,1,1);
dv(3) = Param.v{3}(1,1,2) - Param.v{3}(1,1,1);
dv123 = prod(dv);    

% Auxiliaries (Spectral)
if n == 0
    for i = 1:Param.d 
        mu = 2*pi*1i/Param.db(i)*[-Param.M(i)/2:Param.M(i)/2-1].';
        Param.AuxI1{i} = sqrt(Param.db(i))./mu;
        Param.AuxI1{i}(Param.M(i)/2+1) ...
        = 1/2*(Param.b{i}(2)^2-Param.b{i}(1)^2)/sqrt(Param.db(i));       
        Param.AuxI2{i} ...    
        = ((Param.b{i}(2)^2-Param.b{i}(1)^2)*mu - 2*Param.db(i)) ...
          ./(sqrt(Param.db(i))*mu.^2);
        Param.AuxI2{i}(Param.M(i)/2+1) ...
        = 1/3*(Param.b{i}(2)^3-Param.b{i}(1)^3)/sqrt(Param.db(i));
    end 
end 

% Mass 
Param.MassQuadrature(n+1,1) = dv123*sum(sum(sum(f)));    
Param.MassSpectral(n+1,1) ...
= Param.I0{1}*Param.I0{2}*Param.I0{3} ...
  *fs(Param.M(1)/2+1,Param.M(2)/2+1,Param.M(3)/2+1);

% Momentum
Param.MomentumQuadrature(n+1,1) = dv123*sum(sum(sum(Param.v{1}.*f)));
Aux = fs(:,Param.M(2)/2+1,Param.M(3)/2+1);
Param.MomentumSpectral(n+1,1) ...
= real(Param.I0{2}*Param.I0{3}*sum(Param.AuxI1{1}.*Aux));    

% Energy 
Param.EnergyQuadrature1(n+1,1) = dv123*sum(sum(sum(Param.v{1}.^2.*f)));
Param.EnergyQuadrature2(n+1,1) = dv123*sum(sum(sum(Param.v{2}.^2.*f)));  
Param.EnergyQuadrature3(n+1,1) = dv123*sum(sum(sum(Param.v{3}.^2.*f)));  
Param.EnergyQuadrature(n+1,1) ...
= Param.EnergyQuadrature1(n+1,1) ...
  + Param.EnergyQuadrature2(n+1,1) + Param.EnergyQuadrature3(n+1,1);
Aux = fs(:,Param.M(2)/2+1,Param.M(3)/2+1);
Param.EnergySpectral1(n+1,1) ...
= real(Param.I0{2}*Param.I0{3}*sum(Param.AuxI2{1}.*Aux));
Aux = fs(Param.M(1)/2+1,:,Param.M(3)/2+1).';
Param.EnergySpectral2(n+1,1) ...
= real(Param.I0{1}*Param.I0{3}*sum(Param.AuxI2{2}.*Aux));
Aux = reshape(fs(Param.M(1)/2+1,Param.M(2)/2+1,:),Param.M(3),1);
Param.EnergySpectral3(n+1,1) ...
= real(Param.I0{1}*Param.I0{2}*sum(Param.AuxI2{3}.*Aux));
Param.EnergySpectral(n+1,1) ...
= Param.EnergySpectral1(n+1,1) ...
  + Param.EnergySpectral2(n+1,1) + Param.EnergySpectral3(n+1,1);

% Entropy
g = f(find(f > 0));
Param.Entropy(n+1,1) = dv123*sum(sum(sum(g.*log(g))));

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function MyPlotPrint_Simplified3d(n,f,Param)

% Auxiliary function 
Param.AsColumn = inline('reshape(f,prod(size(f)),1)');

% Mass 
close all 
FactorMin = 1 - 1e-10;
FactorMax = 1 + 1e-10;
Times = Param.Times;
Values = Param.MassSpectral;
MinValues = FactorMin*min(Values);
MaxValues = FactorMax*max(Values);
close all 
plot(Times,Values, ...
'LineWidth',Param.LineWidth, ...
'Color',Param.ColorMyRed)   
grid on
xlabel('t')
title({['Landau equation (Test problem D, d = 3)']; ...
       ['Mass (RKM = ',num2str(Param.RKM), ...
        ', dt = ',num2str(Param.dt,'%10.0e'),')']})   
axis([Param.t0,Param.T,MinValues,MaxValues])
yticks([MinValues,Values(1),MaxValues])
ytickformat('%10.10e')
set(gcf,'Color',Param.ColorMyGrey)
set(gca,'Color',Param.ColorMyGrey)
set(gca,'FontSize',Param.FontSize)
CheckMassConservation = norm(Param.AsColumn(Values-Values(1)))
print -f1 -r600 -depsc Plot_ProblemD_3d_Mass

% Momentum
close all 
FactorMin = 1;
FactorMax = 1;
Times = Param.Times;
Values = Param.MomentumSpectral;
MinValues = min(Values);
MaxValues = max(Values);
close all 
plot(Times,Values, ...
'LineWidth',Param.LineWidth, ...
'Color',Param.ColorMyRed)   
grid on
xlabel('t')
title({['Landau equation (Test problem D, d = 3)']; ...
       ['Momentum (RKM = ',num2str(Param.RKM), ...
        ', dt = ',num2str(Param.dt,'%10.0e'),')']})   
axis([Param.t0,Param.T,MinValues,MaxValues])
yticks([MinValues,MaxValues])
set(gcf,'Color',Param.ColorMyGrey)
set(gca,'Color',Param.ColorMyGrey)
set(gca,'FontSize',Param.FontSize)
CheckMomentumConservation = norm(Param.AsColumn(Values-Values(1)))
print -f1 -r600 -depsc Plot_ProblemD_3d_Momentum

% Energy
close all 
FactorMin = 0.9;
FactorMax = 1.1;
Times = Param.Times;
Values = Param.EnergySpectral;
MinValues = min(Values);
MaxValues = max(Values);
close all 
plot(Times,Values, ...
'LineWidth',Param.LineWidth, ...
'Color',Param.ColorMyRed)   
grid on
xlabel('t')
title({['Landau equation (Test problem D, d = 3)']; ...
       ['Energy (RKM = ',num2str(Param.RKM), ...
        ', dt = ',num2str(Param.dt,'%10.0e'),')']})  
axis([Param.t0,Param.T,MinValues,MaxValues])
yticks([MinValues,MaxValues])
ytickformat('%10.11e')
set(gcf,'Color',Param.ColorMyGrey)
set(gca,'Color',Param.ColorMyGrey)
set(gca,'FontSize',Param.FontSize)
CheckEnergyConservation = norm(Param.AsColumn(Values-Values(1)))
print -f1 -r600 -depsc Plot_ProblemD_3d_Energy

% Entropy
close all 
FactorMin = 0.9;
FactorMax = 1.1;
Times = Param.Times;
Values = Param.Entropy;
MinValues = min(Values);
MaxValues = max(Values);
close all 
plot(Times,Values, ...
'LineWidth',Param.LineWidth, ...
'Color',Param.ColorMyRed)   
grid on
xlabel('t')
title({['Landau equation (Test problem D, d = 3)']; ...
       ['Entropy (RKM = ',num2str(Param.RKM), ...
        ', dt = ',num2str(Param.dt,'%10.0e'),')']})   
axis([Param.t0,Param.T,MinValues,MaxValues])
yticks([MinValues,MaxValues])
set(gcf,'Color',Param.ColorMyGrey)
set(gca,'Color',Param.ColorMyGrey)
set(gca,'FontSize',Param.FontSize)
print -f1 -r600 -depsc Plot_ProblemD_3d_Entropy

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

