function [rhoi, Yi] = Stoch_Initial_Perturb_rho (XP,YP,dx,dy,a1,b,b1,app,Nbcs,NodeP,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
%       VAR:    Initial variances of stochastic modes (from eigenvalues of matrix)
%
% 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 = prod(size(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))];
rhoind = Nbcs+(1:length(idsP))-1;

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

tmp = zeros(Ny,Nx);

for n = 1:app.S(1)
    tmp(plateau_index) = Umat(:,n).*app.sigma_nd(3).*wgt;
    rhoi(rhoind,n) = tmp(idsP);
end;

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

Yi = zeros (app.MC,app.S(1));

for n=1:app.S(1)
    Yi(:,n) = sqrt(Svals(n,n)).*normrnd(0,1,app.MC,1);
end

Yi = Yi.*sqrt(app.varscale); % if M is realiz matrix, then P=M*M^T/n_r.

%keyboard;
