function [ui, vi, variance] = Stoch_Initial_Perturb(XU,YU,XV,YV,Nbcs,NbcsD,NodeP,NUMMODEBASIS)
%function [ui, vi, variance] = Stoch_Initial_Perturb(XU,YU,XV,YV,Nbcs,NbcsD,NodeP,NUMMODEBASIS)
%
% INPUTS:
%       XU:     Matrix of x-coordinates corresponding to U-velocity CV
%              centers. XU(idsu) used as 'x' input to Fu. (see meshgrid)
%       XV:     Matrix of x-coordinates corresponding to V-velocity CV
%              centers. YU(idsu) used as 'y' input to Fu. (see meshgrid)
%       YU:     Matrix of y-coordinates corresponding to U-velocity CV
%              centers. XV(idsv) used as 'x' input to Fv. (see meshgrid)
%       YV:     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
%       NbcsD: Number of Dirichlet boundary conditions
%       NodeP:  Pressure node-numbering matrix (see NodePad.m)
%       NUMMODEBASIS: Size of stochastic dimension
%
% OUTPUTS:
%	ui:     Stochastic u-velocity modes
%	vi: 	Stochastic v-velocity modes
%       variance:    Initial varianceiances of stochastic modes (from eigenvalues of matrix)
%
% Authors: Themis Sapsis, Matt Ueckermann
% Last Modified: Tapovan Lolla, 13 Feb. 2014.

% The input coordinate matrices should be for interior points only!
if all(size(XU,2)==size(NodeP,2))
    error('The input coordinate matrices should be for interior points only!')
end

%
[Ny Nx] = size(NodeP);
dx = XU(1,2)-XU(1);
minx = min(XU(:))-dx;
maxx = max(XU(:))+dx;
a = maxx-minx;	% Width of domain.

dy = YV(1)-YV(2);
miny = min(YV(:))-dy;
maxy = max(YV(:))+dy;
b = maxy-miny;	% Height of domain.

Ncoarse = 30*30;  	
imax_isp=max(min(ceil(Ncoarse^(Nx/(Nx+Ny)))-2,60),15); % No. of interior coarse grid points in x-direction. (Limit between 15-60)

% For the coarse grid used to create initial perturbations
jmax_isp=max(min(ceil(Ncoarse^(Ny/(Nx+Ny)))-2,60),15); % No. of interior coarse grid points in x-direction. (Limit between 15-60)

% In general set it such that the product is less than 5*10^3 for a reasonable comp. time
DDx_isp = a/(imax_isp+1);	% x-spacing between coarse grid.
DDy_isp = b/(jmax_isp+1);	% y-spacing between coarse grid.

[X_isp_c,Y_isp_c] = meshgrid(minx:DDx_isp:maxx,miny:DDy_isp:maxy);
[Ny_c Nx_c] = size(X_isp_c);

[XP YP] = meshgrid([minx,linspace(minx+dx/2,maxx-dx/2,Nx-2),maxx],...
    [miny,linspace(miny+dy/2,maxy-dy/2,Ny-2),maxy]);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% wgt is a mollifier function to have consistent boundary conditions %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wgt = ones(Ny,Nx);   % Ones in interior
wgt(NodeP<=Nbcs)=0;  % Zeros on all boundaries

% Now interpolate onto coarse grid
wgt = interp2(XP,YP,wgt,X_isp_c,Y_isp_c);

kmax = 5;  % Averaging passes. 
           % Change to higher value if you want smoother transition.
for k = 1:kmax
   % Create shift matrices.
   wgt_int = wgt(2:end-1,2:end-1);	 % Interior.

   wgt_right = wgt(2:end-1,3:end);	 % Right.
   wgt_left = wgt(2:end-1,1:end-2);      % Left.
   wgt_top = wgt(1:end-2,2:end-1);	 % Top.
   wgt_bottom = wgt(3:end,2:end-1);      % Bottom.

   wgt_rt = wgt(1:end-2,3:end);	        % Right Top.
   wgt_lt = wgt(1:end-2,1:end-2);	% Left Top.
   wgt_rb = wgt(3:end,3:end);	        % Right Bottom.
   wgt_lb = wgt(3:end,1:end-2);	        % Left Bottom.
   
   wgt(2:end-1,2:end-1) = wgt_int.*(wgt_right + wgt_left + wgt_top + wgt_bottom + wgt_rt + wgt_lt + wgt_rb + wgt_lb)/8;
end
clear wgt_int wgt_right wgt_left wgt_top wgt_bottom wgt_rt wgt_lt wgt_rb wgt_lb

pcolor(X_isp_c,Y_isp_c,wgt); shading flat;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Correlation Function Calculation %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

b = 1;	% Length Scale.

[X1_isp_c,X2_isp_c] = meshgrid(X_isp_c(:),X_isp_c(:));	% Meshgrid is used as a trick to compute the matrix of differences of x coordinates.
[Y1_isp_c,Y2_isp_c] = meshgrid(Y_isp_c(:),Y_isp_c(:));	% Meshgrid is used as a trick to compute the matrix of differences of y coordinates.
[wgt_1,wgt_2] = meshgrid(wgt(:),wgt(:));		% Meshgrid wgt for tensor product computation. 
r = sqrt((X1_isp_c-X2_isp_c).^2 + (Y1_isp_c-Y2_isp_c).^2);	% Compute matrix of distances.
C = (1 + b*r + b^2*(r.^2)/3).*exp(-b*r)*10^-6;		% Compute Correlation Matrix.
C = C.*(wgt_1.*wgt_2)*DDx_isp*DDy_isp;

[smf_cor,D1] = eig(C);
D1 = flipud(diag(D1));		% Flipping to ensure high energy modes are at the beginning.
smf_cor = fliplr(smf_cor);	

for k=1:NUMMODEBASIS
    variance(k)= D1(k);	% variance of modes.

    if k<=12
        subplot(3,4,k);
        pcolor(X_isp_c,Y_isp_c,reshape(smf_cor(:,k),jmax_isp+2,imax_isp+2));
        shading interp;
        hold on;
    end

    [FX,FY] = gradient(reshape(smf_cor(:,k),Ny_c,Nx_c),DDx_isp,DDy_isp);
    U_T = FY;
    V_T = -FX;    
    
    UU_COR(:,:,k)=interp2(X_isp_c,Y_isp_c,U_T,XU,YU,'linear',0);
    VV_COR(:,:,k)=interp2(X_isp_c,Y_isp_c,V_T,XV,YV,'linear',0);
end

[Ny Nx] = size(XU);
ui = [zeros(Nbcs,NUMMODEBASIS);reshape(UU_COR,(Ny)*(Nx),NUMMODEBASIS)];
[Ny Nx] = size(XV);
vi = [zeros(Nbcs,NUMMODEBASIS);reshape(VV_COR,(Ny)*(Nx),NUMMODEBASIS)];

variance = variance/variance(1);	% Normalize varianceiances.

