% Naviers Stokes Flow Lock Exchange test case

LockExchangeSetup_numdet;

%Set default number of Monte-Carlo samples
if ~isfield(app,'MC')
    app.MC = 40000;
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));

%Initialize stochastic coefficients as normal random with variance as calculated
%in Stoch_Initial_Perturb
YYt=zeros(app.MC,app.S(1));
randn('state',0);
for i=1:app.S(1)
%     tmp = normrnd(0, (sqrt(eps) * exp(-i)),app.MC/2,1);
%     YYt(:,i) = [tmp;-tmp];    
    YYt(:,i) = 1e-3*sqrt(eps)*exp(0)*repmat([app.dref/100/2 - app.drho/2]',app.MC/(app.S+1),1);
end
%Need these for Orthonormalization
uid = Nbcs+1:length(u);
vid = Nbcs+1:length(v);
pid = Nbcs+1:length(rho);

% Now initialize density
rhoi(pid,1) = rho(pid) / app.drho * 2;
OO = ones(app.MC/4,1);
YYt(:,1) = repmat([app.dref/100/2 - app.drho/2]',app.MC/(app.S+1),1);

% for i = 2 : app.S(1)
%     if i == 2 
%         tmp = cos(XP / 2 * pi) .* cos(YP / 2 * pi);
%     else
%         tmp = sin(YP * pi * (i - 2)).*cos(XP / 2 * pi);
%     end
%     tmp = sin(YP * pi * app.Ny / i / 2) .* sin(XP * pi * app.Nx / i / 2);
%     rhoi(pid,i) = tmp(idsP);    
% end
%
if app.S(1) > 1
    tmp = sin(pi * YP * (2 - 1));
    tmp = tmp(idsP) .* (app.drho / 2 - abs(rho(pid))) .* sign(rho(pid));
    tmp(tmp<0) = 0;
    rhoi(pid,2) = tmp .* sign(YP(idsP));
end
if app.S(1) > 2
    tmp = sin(pi * YP * (2 - 1));
    tmp = tmp(idsP) .* (app.drho / 2 - abs(rho(pid))) .* sign(rho(pid));
    tmp(tmp>0) = 0;
    rhoi(pid,3) = tmp .* sign(YP(idsP));
    rhoi(1:Nbcs,:) = 0;
end

%%
%%
% keyboard
% OO = ones(app.MC/4,1);
% sstart = 1;
% send = ceil(app.Ny / 3);
% for i = 1 : app.S(1)
%     tmp = zeros(size(XP));
%     tmp(sstart : min(send, app.Ny), 1 : app.Nx) = 1;
%     tmp = rho(NodeP(2 : app.Ny + 1, 2 : app.Nx + 1)) .* tmp / app.drho * 2;    
%     rhoi(pid,i) = tmp(idsP);    
%     sstart = send + 1;
%     send = send + ceil(app.Ny / 3);
%     YYt(:,i) = repmat([app.dref/100/2 - app.drho/2]',app.MC/4,1);
% end
% keyboard
%
% for i=2:app.S(1)
%     tmp = normrnd(0, (sqrt(eps) * exp(-i)),app.MC/2,1);
%     YYt(:,i) = [tmp;-tmp];    
% end
% %Make sure dummy variances are uncorrelated
% for i=2:app.S(1)
%     %Commented below in order to check with Themis' code
%     for j=1:app.S(1)
%         if i~=j
%             coef = YYt(:,i)' * YYt(:,j) / (app.MC-1);
%             YYt(:,i) = YYt(:,i) - coef * YYt(:,j) / std(YYt(:,j)).^2;
%         end
%     end
% end

%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;
%     YYt(:,i) = YYt(:,i) * nrm;
% end
for i = 1:1
    M = (ui(uid, :)' * ui(uid, :) + vi(vid, :)' * vi(vid, :) ...
        + rhoi(pid, :)' * rhoi(pid, :)) * dx * dy;
    [VC DC] = eig(M);
    YYt = fliplr(YYt * VC * sqrt(DC));
    DC = diag (1 ./ diag(sqrt(DC)));
    ui(uid, :) = fliplr(ui(uid, :) * VC * DC);
    vi(vid, :) = fliplr(vi(vid, :) * VC * DC);
    rhoi(pid, :) = fliplr(rhoi(pid, :) * VC * DC);
    Pi(pid, :) = fliplr(Pi(pid, :) * VC * DC);
%     keyboard
end
CYY=cov(YYt);
[VC,DC] = eig(CYY) ; %Diagonal Covariance, and Diagonal Vectors
ui(uid,:) = fliplr(ui(uid, :) * VC);
vi(vid,:) = fliplr(vi(vid, :) * VC);
Pi(pid,:) = fliplr(Pi(pid, :) * VC);
rhoi(pid,:) = fliplr(rhoi(pid, :) * VC);
YYt=fliplr(YYt * VC);

rhoi(1:Nbcs,:) = 0;