1

clc; close all

% [U, S, V] = svd(Asnap, 'econ');
% plot(diag(S)/sum(diag(S)), 'ro')
% % based on this plot, we set r = 10


%% DMD
X1 = Asnap(:, 1:end-1);
X2 = Asnap(:, 2:end);
r = 10;

[U, S, V] = svd(X1, 'econ');

% rank-r truncation
Ur = U(:, 1:r);
Sr = S(1:r, 1:r);
Vr = V(:, 1:r);

Atilde = Ur'*X2*Vr/Sr;

% eigenvalues and eigenvectors
[W, D] = eig(Atilde);

% DMD modes
Phi = X2*Vr/Sr*W;

% eigenvalues
lamdba = diag(D);

% frequencies
omega = log(lamdba)/(iconv*Deltat);

% initial condition
x1 = Asnap(:, 1);
b = Phi\x1;

%% Reconstruction
time_vector = (0:iconv*Deltat:nt*Deltat);
time_dynamics = zeros(r, size(Asnap, 2));
for iter = 1:size(Asnap, 2)
    time_dynamics(:, iter) = (b.*exp(omega*time_vector(iter)));
end
Asnap_dmd = Phi*time_dynamics;

% get vorticity vector at last time step
vorticity_dmd = real(Asnap_dmd(:, 100));

% reshape for plotting
vorticity_dmd = reshape(vorticity_dmd, ncx-1, ncy-1);
vorticity_dmd = vorticity_dmd';

% plot
fz=figure;
az=axes;
maxval=max(max(vorticity_dmd));
minval=min(min(vorticity_dmd));
vals=linspace(minval,maxval,100);
contour(X,Y,vorticity_dmd,vals);
set(az,'DataAspectRatio',[1 1 1],'XLim',[0 lx],'YLim',[0 ly]);
title(sprintf(' Vorticity (DMD) (%g...%g); Re_b=%g',...
                minval,maxval,Re));
colorbar;