function [YYt ui vi Pi rhoi] = DOorthnorm(app, dx, dy, YYt, ui, uid, vi, vid, Pi, pid, rhoi)
% function [YYt u v P rho a] = DOorthnorm(YYt, ui, uid, vi, vid, Pi, pid, rhoi)
% This function orthonormalizes the DO modes while preserving the original
% realizations (upto a correction for the stochastic energy) and the
% stochastic energy.
%
% INPUTS:
%   dx:         Size of x discretization
%   dy:         Size of y discretization
%   YYt(MC,S):  Stochastic coefficients
%   ui(:,S):    Mode shapes for u-velocity
%   uid:        Vector of ids for interior u-velocities
%   vi(:,S):    Mode shapes for v-velocity
%   vid:        Vector of ids for interior v-velocities
%   Pi(:,S):    Mode shapes for Pressure
%   Pid:        Vector of ids for interior Pressures
%   rhoi(Optional): Mode shapes for density
%
% OUTPUTS:
%   YYt(MC,S):  Stochastic coefficients
%   ui(:,S):    Mode shapes for u-velocity
%   vi(:,S):    Mode shapes for v-velocity
%   Pi(:,S):    Mode shapes for Pressure
%   rhoi(Optional): Mode shapes for density
%
% Written by: Matt Ueckermann
% Last Modified: Tapovan Lolla, on 4/25/2014.

if ~isfield(app,'IP'),app.IP=[1,1,1];end
if ~isfield(app,'Sb'), app.Sb = 0;end

%% Old orthonormalization procedure
if isfield(app,'GSorth')
    disp('Gram Schmidt Orthonormalization has been removed from the FV Framework, proceeding with regular orthonormalization instead ...')
end

% First rotate the stochastic coefficients so that we have uncorrelated samples
Sintid = app.Sb + 1:app.S; % IDs of modes to be orthonormalized.
CYY = cov(YYt(:, Sintid)); % Covariance matrix of coefficients.
[VC, tmp] = eig(CYY) ;     % Diagonalizing the covariance matrix.

ui(:, Sintid) = fliplr(ui(:, Sintid) * VC); % Need to multiply modes by Eigenvector matrix to preserve realizations. YYt*[ui',vi',rhoi'].
vi(:, Sintid) = fliplr(vi(:, Sintid) * VC);
Pi(:, Sintid) = fliplr(Pi(:, Sintid) * VC);
if nargin == 11
    rhoi(:,Sintid) = fliplr(rhoi(:, Sintid) * VC);
else
    rhoi = 0;
end
YYt(:,Sintid) = fliplr(YYt(:,Sintid) * VC); % Multiply the coefficients by VC to make them uncorrelated.

% Next orthonormalize the modes (ensure their inner product is I). 
% Warning: THIS IP IS NOT A SIMPLE MATRIX TRANSPOSE MATRIX MULTIPLICATION.
% Calculate current inner product. Note that this operation uses only inner elements of modes, excluding boundaries.

if nargin == 11
    M = (app.IP(1)* (ui(uid, :)' * ui(uid, :)) + app.IP(2)* (vi(vid, :)' * vi(vid, :)) ...
        + app.IP(3) * (rhoi(pid, :)' * rhoi(pid, :)))*(dx*dy);
else
    M = (ui(uid, :)' * ui(uid, :) + vi(vid, :)' * vi(vid, :))*(dx*dy);
end

% Calculate EigenValues and Vectors of M, to multiply modes matrix with them.

[VC DC] = eig(M);
YYt = fliplr(YYt * VC * sqrt(DC)); % Update coefficients to ensure realizations are preserved.
DC = diag (1./ diag(sqrt(DC)));    % Inverse of eigenvalue matrix needed for the next two lines.
ui = fliplr(ui * VC * DC);         % This multiplication by VC and DC ensures that the inner product of the new modes is I.
vi = fliplr(vi * VC * DC);
if nargin == 11
    rhoi = fliplr(rhoi * VC * DC);
end
Pi = fliplr(Pi * VC * DC);

% Now rotate all the modes such that the boundaries are orthogonal once again.
% WARNING: The boundaries have a different Inner Product.
% It is aimed to only make boundaries orthogonal to other boundaries, not the interior.
% The Inner Product for boundaries uses the usual dot product norm.

Nbcs = length(ui) - length(uid);
if nargin == 11
    [VC,DC] = eig(ui(1:Nbcs,:)'*ui(1:Nbcs,:) + vi(1:Nbcs,:)'*vi(1:Nbcs,:) + ...
        rhoi(1:Nbcs,:)'*rhoi(1:Nbcs,:));
    rhoi = fliplr(rhoi*(VC));
else
    [VC,DC] = eig(ui(1:Nbcs,:)'*ui(1:Nbcs,:) + vi(1:Nbcs,:)'*vi(1:Nbcs,:));
end
ui = fliplr(ui*(VC)); % Note that the entire modes matrix is affected by this operation, not just the boundaries.
vi = fliplr(vi*(VC)); % It may be verified that this operation preserves orthonormality of inner elements of the modes.
Pi = fliplr(Pi*(VC));
YYt = fliplr(YYt*(VC)); % Modify the coefficients matrix also, to preserve realizations.

% Rotate the stochastic coefficients to achieve uncorrelated behavior.

CYY = cov(YYt(:, Sintid));

[VC,DC] = eig(CYY) ; %Diagonal Covariance, and Diagonal Vectors

ui(:, Sintid) = fliplr(ui(:, Sintid) * VC); % This operation preserves orthonormality of the modes.
vi(:, Sintid) = fliplr(vi(:, Sintid) * VC);
Pi(:, Sintid) = fliplr(Pi(:, Sintid) * VC);
if nargin == 11
    rhoi(:, Sintid) = fliplr(rhoi(:, Sintid) * VC);
end

% In this final step, we allow for mismatch in the realizations, to ensure the stochastic energy is preserved.
% Also correct for, or make sure that, stochastic energy is preseved.
YYt(:, Sintid) = fliplr(YYt(:,Sintid)*VC*sqrt(sum(tmp(:))/sum(DC(:))));

