% Stochastic Lock Exchange test case
% (Navier-Stokes with Boussinesq Approximation)

LockExchangeSetup;

%Change the initialization of rho
if ~isfield(app,'var')
    app.var = 0.05;
end;
rho=[bcsDrho; bcsNrho; 9.81*0.5*tanh((XP(idsP)-a/2)*2/app.var)];

%Set default number of Monte-Carlo samples
if ~isfield(app, 'MC'), error('app.MC not specified'); end;
if ~isfield(app, 'S'), error('app.S not specified'); 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=zeros(length(rho),app.S(1));


%Initialize stochastic coefficients as normal random with variance as calculated
%in Stoch_Initial_Perturb
Yi=zeros(app.MC,app.S(1));
for i=1:app.S(1)
    Yi(:,i) = rand(app.MC, 1)*2-1; % Uniform distribution over [-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
for i = 1:app.S(1)
    tmp = cos(pi * (YP*2-1) * (max(app.S(1),5) - i));
    rhoi(pid,i) = tmp(idsP) .* ...
        (1 - sum([sum(abs(rhoi(pid,1:i)), 2), abs(rho(pid))], 2) / 9.81 * 2);
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;
    Yi(:,i) = Yi(:,i)*nrm;
end
Yi = sortrows(Yi,1);
