% Naviers Stokes Flow Lock Exchange test case

LockExchangeSetup_numdet;

%Set default number of Monte-Carlo samples
if ~isfield(app,'MC')
    app.MC = 1000;
end
if ~isfield(app,'S') 
    app.S(1)=1;
end
if app.S(1)>3
    display('Only works for S=[1,2,3])')
    app.S(1)=3;
end

%%Initialize modes
ui = zeros(length(u),app.S(1));
vi = zeros(length(v),app.S(1));
Pi = zeros(length(P),app.S(1));
rhoi = ones(length(rho),app.S(1));
%%
k = app.k;%app.PlotIntrvl * 8;
%%
rho(:) = 0;
u(:) = 0;
v(:) = 0;
for i = 1:length(app.dref)
    det{i} = load(sprintf('%s/Save/%s%d/%05d.mat',d0,app.ref,app.dref(i),...
        k/app.PlotIntrvl),'rho','u','v');
    rho = rho + det{i}.rho;
    u = u + det{i}.u;
    v = v + det{i}.v;
end
rho = rho / (length(app.dref));
u = u / (length(app.dref));
v = v / (length(app.dref));
for i = 1:min(length(app.dref), app.S(1))
    rhoi(:,i) = det{i}.rho - rho;
    ui(:,i) = det{i}.u - u;
    vi(:,i) = det{i}.v - v;
end
%Orhonormalize
TMP = (rhoi'*rhoi + ui'*ui + vi'*vi)*dx*dy;
[VT DT] = eig(TMP);
rhoi = fliplr(rhoi * VT);
ui = fliplr(ui * VT);
vi = fliplr(vi * VT);

TMP = sqrt(1./diag((rhoi'*rhoi + ui'*ui + vi'*vi)*dx*dy));
rhoi = rhoi * diag(TMP);
ui = ui * diag(TMP);
vi = vi * diag(TMP);

%check
% i = 2;
% LL = (rhoi'*(det{i}.rho - rho) + ui'*(det{i}.u - u) + vi'*(det{i}.v - v))*dx*dy
% tmp = rhoi * LL + rho - det{i}.rho;
% pcolor(tmp(NodeP)),colorbar
% figure
% tmp = ui * LL + u - det{i}.u;
% pcolor(tmp(Nodeu)),colorbar
% figure
% tmp = vi * LL + v - det{i}.v;
% pcolor(tmp(Nodev)),colorbar

%%
uid = Nbcs+1:length(u);
vid = Nbcs+1:length(v);
pid = Nbcs+1:length(rho);
%Ortho-Normalization of basis fields to avoid round-off errors
for i=1:app.S(1)
    %Commented below in order to check with Themis' code
    for j=1:app.S(1)
        if i~=j
            coef = InProd(ui(:,i),ui(:,j),dx,dy,uid) + ...
                InProd(vi(:,i),vi(:,j),dx,dy,vid) + ...
                InProd(rhoi(:,i),rhoi(:,j),dx,dy,pid);
            ui(:,i) = ui(:,i) - coef * ui(:,j);
            vi(:,i) = vi(:,i) - coef * vi(:,j);
            rhoi(:,i) = rhoi(:,i) - coef * rhoi(:,j);
        end
    end
    nrm=sqrt(InProd(ui(:,i),ui(:,i),dx,dy,uid) + ...
        InProd(vi(:,i),vi(:,i),dx,dy,vid) + ...
        InProd(rhoi(:,i),rhoi(:,i),dx,dy,pid));
    ui(:,i) = ui(:,i) / nrm;
    vi(:,i) = vi(:,i) / nrm;
    rhoi(:,i) = rhoi(:,i) / nrm;
end
rhoi(1:Nbcs,:) = 0;

%% Now find the corresponding coefficients
for i = 1:length(app.dref)
    YYt(i,:) = ((det{i}.rho - rho)' * rhoi + (det{i}.u - u)' * ui + ...
        (det{i}.v - v)' * vi) * dx * dy;
end
%%
% keyboard