close all
clear all
clc

%% Example script to undertake the POD decomposition of a flow field

%% Make directories and load data

FlowFieldDir = 'FlowFieldDirectory/';
OutputDir = 'OutputDirectory/';
mkdir(OutputDir);

load([FlowField Dir 'filename.mat']); % Flow field with X,Y,U,V

%% Preallocate matrices

% Get matrix sizes, remove NaNs
sx = size(U,1); sy = size(U,2); st = size(U,3);
U(isnan(U)) = 0; V(isnan(V)) = 0;

% Data matrix for decomposition
A = zeros(sx*sy*2,st);
for i = 1:st

    u_vect = reshape(U(:,:,i),[],1);
    v_vect = reshape(V(:,:,i),[],1);
    A(:,i) = [u_vect;v_vect]; % Stack velocity components

end

%% Perform SVD

[UU,S,VV] = svd(A,'econ');

%% Reshape matrices

u_mat = UU(1:sx*sy,:); v_mat = UU((sx*sy+1):end,:);

u_mode = reshape(u_mat,sx,sy,st);
v_mode = reshape(v_mat,sx,sy,st);

%% Cumulative energy and 90% index

s1 = diag(S); cumenergy = cumsum(s1.^2)/sum(s1.^2);

[val,ind90] = min(abs(cumenergy-0.9));
if val < 0.9
    ind90 = ind90+1;
end

save([OutputDir 'SVD_MeanInclusion.mat'],...
        'A','cumenergy','DataName','ExpName','ind90','S','s1','st',...
        'sx','sy','st','time','u_mode','v_mode','X','Y','UU','VV','-v7.3');

%%