function [ui, vi, VAR] = Stoch_Initial_Perturb(XU,YU,XV,YV,Nbcs,NbcsD,NodeP,NUMMODEBASIS)
%function [ui, vi, VAR] = 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
%       VAR:    Initial variances of stochastic modes (from eigenvalues of matrix)
%
% Authors: Themis Sapsis and Matt Ueckermann

% 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;
dy = YV(1)-YV(2);
miny = min(YV(:))-dy;
maxy = max(YV(:))+dy;
b=maxy-miny;

Ncoarse = 30*30;
imax_isp=max(min(ceil(Ncoarse^(Nx/(Nx+Ny)))-2,60),15);  % For the coarse grid used to create initial perturbations
jmax_isp=max(min(ceil(Ncoarse^(Ny/(Nx+Ny)))-2,60),15);   % 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);
DDy_isp = b/(jmax_isp+1);
[X_isp_c,Y_isp_c] = meshgrid(minx:DDx_isp:maxx,miny:DDy_isp:maxy);
[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]);
%First run the setup script, to setup the geometry

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% wgt is a mollifier function to have consistent boundary conditions %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wgt=ones(Ny,Nx); %Ones in interior
wgt(NodeP<=Nbcs)=0; %Zeros on DIRICHLET (only) boundaries
wgt(NodeP<=NbcsD)=1; %Zeros on DIRICHLET (only) boundaries
%Now interpolate onto coarse grid
wgt = interp2(XP,YP,wgt,X_isp_c,Y_isp_c);
for k=1:0 %Change to higher value if you want smoother transition
    wgs=wgt;
for j=2:size(wgt,1)-1
    for i=2:size(wgt,2)-1
    wgt(j,i)=wgs(j,i)*(wgs(j-1,i)+wgs(j-1,i+1)+wgs(j-1,i-1)+wgs(j,i-1)+wgs(j,i+1)+wgs(j+1,i-1)+wgs(j+1,i)+wgs(j+1,i+1))/8;
    end
end
end

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

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

b=5;


for m=1:(imax_isp+2)*(jmax_isp+2)
    
    for n=1:(imax_isp+2)*(jmax_isp+2)
        
        j_m=mod(m-1,jmax_isp+2);
        i_m=(m-1-j_m)/(jmax_isp+2);
        
        j_n=mod(n-1,jmax_isp+2);
        i_n=(n-1-j_n)/(jmax_isp+2);
         
        xm=(i_m)*DDx_isp;
        ym=(j_m)*DDy_isp;
        
        xn=(i_n)*DDx_isp;
        yn=(j_n)*DDy_isp;
        
        r=norm([xm-xn,ym-yn]);
        
        CC(m,n)=(1+b*r+b^2*r^2/3)*exp(-b*r)*10^-6;   % Definition of Correlation function
        
        CC(m,n)=wgt(j_n+1,i_n+1)*wgt(j_m+1,i_m+1)*CC(m,n)*DDx_isp*DDy_isp;
        
%         CC(m,n)=cos(r/b)*exp(-r^2/b^2);
        
        
    end
end
display('Correlation matrix created')
[smf_cor,D]=eig(CC);
display('Eigen functions found')
for m=1:(imax_isp+2)*(jmax_isp+2)
D1(m)=D(m,m);
end
clear D
figure
% for k=1:12
%     
%     subplot(3,4,k); pcolor(X_Mod_gr,Y_Mod_gr,reshape(smf_cor(:,m-k+1),imax_cor-1,jmax_cor-1)'); shading interp; hold on;
%     contour(X_Mod_gr,Y_Mod_gr,reshape(smf_cor(:,m-k+1),imax_cor-1,jmax_cor-1)','k')
% end

for k=1:NUMMODEBASIS
    VAR(k)=sqrt(D1(m-k+1));
    if k<=12
    subplot(3,4,k); pcolor(X_isp_c,Y_isp_c,reshape(smf_cor(:,m-k+1),jmax_isp+2,imax_isp+2)); shading interp; hold on;
        
    end
%     contour(XGrid_c,YGrid_c,reshape(smf_cor(:,m-k+1),jmax+2,imax+2)','k')
% 
    [FX,FY] = gradient(reshape(smf_cor(:,m-k+1),jmax_isp+2,imax_isp+2),DDx_isp,DDy_isp);
%     
%     
    U_T=FY;
    V_T=-FX;
     
    if k<=12  
    quiver(X_isp_c,Y_isp_c,U_T,V_T,2);hold off;
    end

    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)];

VAR=VAR/VAR(1);

%clear CC smf_cor D D1;
