function [rhoi, Yi] = Stoch_Initial_Perturb_rho(XP,YP,dx,dy,a1,b,b1,app,...
    Nbcs,NodeRho,idsP)
%function [rhoi, Yi] = Stoch_Initial_Perturb_rho (XP,YP,dx,dy,a1,b,b1,app,Nbcs,NodeP,idsP);
%
% INPUTS:
%       XP:     Matrix of x-coordinates corresponding to V-velocity CV
%              centers. YU(idsu) used as 'y' input to Fu. (see meshgrid)
%       YP:     Matrix of y-coordinates corresponding to V-velocity CV
%              centers. YV(idsv) used as 'y' input to Fv. (see meshgrid)
%       Nbcs:   Number of boundary conditions
%       NodeP:  Pressure node-numbering matrix (see NodePad.m)
%
% OUTPUTS:
%	rhoi:   Stochastic density modes
%
% Authors: Pat Haley and Pierre Lermusiaux

%%
% Find restricted region above plateau.

[jind, iind] = find ((XP<=a1) & (YP>=b1));

iselect = min(iind(:)):max(iind(:));
jselect = min(jind(:)):max(jind(:));

Xwk = XP(jselect,iselect);
Ywk = YP(jselect,iselect);

%-------------------------------------------------------------------------------
% Construct mollifier.
%-------------------------------------------------------------------------------

wgt = ones(size(Xwk)); % a "null" mollifier.

if ~isempty(app.MollifierWidth)
    wgt = Xwk./app.MollifierWidth;
    wgt = (1 - exp(-wgt.*wgt));
end;

wgt = wgt(:);

%-------------------------------------------------------------------------------
% Construct correlation matrix.
%-------------------------------------------------------------------------------

dxscale = a1/2;
dyscale = (b-b1)/2;

nxy = numel(Xwk);
CC = zeros(nxy,nxy);

for m = 1:nxy
    rx = (Xwk(m)-Xwk(:))./dxscale;
    ry = (Ywk(m)-Ywk(:))./dyscale;
    r2 = rx.*rx + ry.*ry;
    %  CC(:,m) = (1+r2./3).*exp(-r2);  % not a covariance
    CC(:,m) = exp(-r2);
end;

%-------------------------------------------------------------------------------
% Construct modes/coefficients from SVD of correlation matrix.
% (SVD seems less suseptible to negative values and EIG).
%-------------------------------------------------------------------------------

[Umat, Svals, Vmat] = svd(CC, 0);

[Ny, Nx] = size(XP);
rhoi = [zeros(Nbcs,app.S(1));zeros(length(idsP),app.S(1))];
pid = Nbcs + (1:length(idsP));

[Iselect,Jselect] = meshgrid (iselect,jselect);
plateau_index = Jselect + (Iselect-1).*Ny;

tmp = zeros(Ny,Nx);

for n = 1:app.S(1)
    tmp(plateau_index) = (app.sigma_nd(3)*sqrt(max(NodeRho(:))-Nbcs))*Umat(:,n);
    rhoi(pid,n) = tmp(idsP);
end;

%-------------------------------------------------------------------------------
% Create coefficients from SVD.
%-------------------------------------------------------------------------------

Svals = Svals(1:app.S(1),1:app.S(1));
Yi_Std = sqrt(diag(Svals))';
Yi_Std = Yi_Std*(app.varscale/sqrt(sum(Yi_Std.^2)));
Yi = repmat(Yi_Std, app.MC, 1).*normrnd(0,1,app.MC,app.S(1));

