% Cylinder Cross Flow setup file

%Define Forcing functions
Fu = @(time,p) 0;
Fv = @(time,p) 0;

%Set size of domain
a    = 16;          b    = 3;
app.Nx = a * app.N;
app.Ny = b * app.N ;
app.dt = a / app.Nx / 3;
app.PlotIntrvl = floor(1 / app.dt);

dx   = a./(app.Nx);    dy   = b./(app.Ny); 

%Set boundary conditions
bcsDrho = [0 -0.0 0 0];bcsNrho = [0];        
bcsDP= [0];         bcsNP = [0 0 0 0];        
bcsOP= [];
bcsDu= [0 1];  bcsNu= [0 0 0];        
bcsOu= [];
bcsDv= [0 0 0 0];  bcsNv=  [0];        

%Set number of timesteps
Nt = ceil(app.T/app.dt);

%% create geometry mask
Nbcs = length(bcsDP) + length(bcsNP);
bnd = [1 1 1 1];
Mask = zeros(app.Ny, app.Nx);

% specify center of the disk
xcdisk = 4 + 0.5;
ycdisk = b/2;
rdisk = 1/2;

xp = linspace(dx/2, a-dx/2, app.Nx  );yp = fliplr(linspace(dy/2, b-dy/2, app.Ny  ));
xu = linspace( dx   ,a-dx , app.Nx-1);yu = fliplr(linspace(dy/2, b-dy/2, app.Ny  ));
xv = linspace(dx/2, a-dx/2, app.Nx  );yv = fliplr(linspace(dy    , b-dy, app.Ny-1));
[XP YP] = meshgrid(xp,yp);[XU YU] = meshgrid(xu,yu);[XV YV] = meshgrid(xv,yv);

%Cylinder
% Mask((XP - xcdisk).^2 + (YP - ycdisk).^2 < rdisk.^2) = 1 ;
% Do a square instead
Mask(abs(XP - xcdisk) <= (rdisk) ...
    & abs(YP - ycdisk) <= (rdisk)) = 1 ;
Mask=(Mask+flipud(Mask));
Mask=Mask.*(max(1,Mask)).^-1;
clear xtemp ytemp 
% imagesc(Mask)
% keyboard
%create Node matrices
[NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask,bnd,1);
NodeRho=NodeP;
NodeP(find(NodeP==1))=-1;NodeP(find(NodeP==3))=1;NodeP(find(NodeP==-1))=3;
% Nodeu(find(Nodeu==3))=-1;Nodeu(find(Nodeu==5))=3;Nodeu(find(Nodeu==-1))=5;
Nodev(find(Nodev==3))=-1;Nodev(find(Nodev==5))=3;Nodev(find(Nodev==-1))=5;
NodeRho(find(NodeRho==3))=-1;NodeRho(find(NodeRho==5))=3;NodeRho(find(NodeRho==-1))=5;

%% Initialize fields
rhoinit=@(X,Y) 0*(X>5).*(X<6).*(Y>2);
pinit = @(X,Y) zeros(size(X));
uinit = @(X,Y) zeros(size(X));
vinit = @(X,Y) zeros(size(X));
%%
base = size(Nodeu,1)-2;
UScale = base./repmat(base-sum((Nodeu==1)),base,1);
UScale(:,1)=[];
rho=[bcsDrho' ;bcsNrho' ;rhoinit(XP(idsP), YP(idsP))];
P = [bcsDP';bcsNP';  pinit(XP(idsP), YP(idsP))];
% u = [bcsDu';bcsNu';  UScale(idsu)];
u = [bcsDu';bcsNu';  uinit(XU(idsu), YU(idsu))];
v = [bcsDv';bcsNv';  vinit(XV(idsv), YV(idsv))];
%Find u divergence:
% % Neu = [length(bcsDu)+1,Nbcs];
% % Nev = [length(bcsDv)+1,Nbcs];
% % VScale = [zeros(Nbcs,1);
% %     Div2D(u,v,Nodeu,Nodev,dx,dy,idsP,Neu,Nev)];
% % VScale = VScale(NodeP(2:app.Ny+1,2:app.Nx+1));
% % % figure,imagesc(VScale),colorbar
% % Vold = zeros(1, app.Nx);
% % for i = 1 : app.Ny - 1
% %     Vold = Vold + dy * VScale(i, :);
% %     VScale(i, :) = Vold;
% % end
% % VScale(app.Ny,:) = [];
% % v = [bcsDv';bcsNv';  VScale(idsv)];


% figure,imagesc(u(Nodeu)),colorbar
% figure,imagesc(v(Nodev)),colorbar
%Check
%  D2D = [zeros(Nbcs,1);Div2D(u,v,Nodeu,Nodev,dx,dy,idsP,Neu,Nev)];
%  figure,imagesc(D2D(NodeP)),colorbar;
%%
%Initialize Modes
ui=zeros(length(u),app.S(1));
vi=zeros(length(v),app.S(1));
Pi=zeros(length(P),app.S(1));
var = zeros(app.S(1),1);

%%
%Create fall-off function around cylinder:
Mask = Nodeu(2 : app.Ny + 1, 2 : app.Nx);
mids = (Mask == 1);
Mask(:) = 1;
Mask(mids) = 0;
Maskv = Nodev(2 : app.Ny, 2 : app.Nx + 1);
midv = (Maskv == 1);
Maskv(:) = 1;
Maskv(midv) = 0;
%Now blend the mask:
for i = 1 : ceil(app.Ny) *2
   Mask(2:end-1,2:end-1) = 0.2 * (Mask(2:end-1,2:end-1) + ...
       Mask(1:end-2, 2:end-1) + Mask(3:end, 2:end-1) + ...
       Mask(2:end-1,1:end-2) + Mask(2:end-1,3:end));
   Maskv(2:end-1,2:end-1) = 0.2 * (Maskv(2:end-1,2:end-1) + ...
       Maskv(1:end-2, 2:end-1) + Maskv(3:end, 2:end-1) + ...
       Maskv(2:end-1,1:end-2) + Maskv(2:end-1,3:end));
   Maskv(midv) = 0;
   Mask(mids) = 0;
end

Mask = Mask(idsu);
Maskv = Maskv(idsv);
%%
uinit = @(x,y,M,N) - 1 / b * sin(pi .* x / a) .* sin(pi .* M .* x / a)...
    .* (cos(pi .* y / b) .* sin(pi .* N .* y / b) + N .* sin(pi .* y / b) .* cos(pi .* N .* y / b));
vinit = @(x,y,M,N)  1 / a * sin(pi .* y / b) .* sin(pi .* N .* y / b)...
    .* (cos(pi .* x / a) .* sin(pi .* M .* x / a) + M .* sin(pi .* x / a) .* cos(pi .* M .* x / a));

i=1;
for ii = 1 : app.S(1)
    for iiy = 1 : ceil(ii/2)
        if (i <= app.S(1))
            ix = ii - iiy + 1;
            iy = iiy;
            nrm = sqrt(((ix^2 + 1))/16*(3/2)^(iy==1) +...
                ((iy^2 + 1))/16*(3/2)^(ix==1));
            ui(:,i) = [zeros(Nbcs,1);
                uinit(XU(idsu), YU(idsu), ix, iy) / nrm .* Mask];
            vi(:,i) = [zeros(Nbcs,1);
                vinit(XV(idsv), YV(idsv), ix, iy) / nrm .* Maskv];
            var(i) = app.var * exp(-(ix+iy) + 2);
            i = i+1;     
%             [i-1 ix iy ii + 1]
        end
        if (i <= app.S(1)) && (ix ~= iy)
            %Note ix and iy interchanged
            iy = ii - iiy + 1;
            ix = iiy;
            nrm = sqrt(((ix^2 + 1))/16*(3/2)^(iy==1) +...
                ((iy^2 + 1))/16*(3/2)^(ix==1));
            ui(:,i) = [zeros(Nbcs,1);
                uinit(XU(idsu), YU(idsu), ix, iy) / nrm .* Mask];
            vi(:,i) = [zeros(Nbcs,1);
                vinit(XV(idsv), YV(idsv), ix, iy) / nrm .* Maskv];
            var(i) = app.var * exp(-(ix+iy) + 2);
            i = i+1;   
%             [i-1 ix iy ii + 1]
        end
    end
end

% CylinderPlots
%Initialize stochastic coefficients as normal random with variance as calculated
%in Stoch_Initial_Perturb
YYt=zeros(app.MC,app.S(1));
randn('state', 0);
for i=1:app.S
    tmp = (normrnd(0,sqrt(var(i)),app.MC/2,1));
    tmp = [-tmp;tmp];
    tmp = tmp * sqrt(var(i)) / std(tmp);
    YYt(:,i) = tmp;
end

%Need these for Orthonormalization
uid = Nbcs+1:length(u);
vid = Nbcs+1:length(v);
%%
%Ortho-Normalization of basis fields to avoid round-off errors
%Also decorrelate modes

for i=1:app.S(1)
     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);
         ui(:,i) = ui(:,i) - coef * ui(:,j);
         vi(:,i) = vi(:,i) - coef * vi(:,j);
         coef = YYt(:,i)' * YYt(:,j) / (app.MC-1);
         YYt(:,i) = YYt(:,i) - coef * YYt(:,j) / std(YYt(:,j)).^2;
         end
     end
    nrm = sqrt(InProd(ui(:,i),ui(:,i),dx,dy,uid)+InProd(vi(:,i),vi(:,i),dx,dy,vid));
    ui(:,i) = ui(:,i)/nrm;
    vi(:,i) = vi(:,i)/nrm;
     
    YYt(:,i) = YYt(:,i) * sqrt(var(i)) / std(YYt(:,i));
end
% for i = 1:2
%     M = (ui(:, :)' * ui(:, :) + vi(:, :)' * vi(:, :)) * dx * dy;
%     [VC DC] = eig(M);
%     YYt = fliplr(YYt * VC * sqrt(DC));
%     DC = diag (1 ./ diag(sqrt(DC)));
%     ui(:, :) = fliplr(ui(:, :) * VC * DC);
%     vi(:, :) = fliplr(vi(:, :) * VC * DC);
%     Pi(:, :) = fliplr(Pi(:, :) * VC * DC);
%     %     keyboard
% end
for k = 1:3
    for i=1:app.S(1)
        for j=1:app.S(1)
            if i~=j
                coef = YYt(:,i)' * YYt(:,j) / (app.MC-1);
                YYt(:,i) = YYt(:,i) - coef * YYt(:,j) / std(YYt(:,j)).^2;
            end
        end
        YYt(:,i) = YYt(:,i) * sqrt(var(i)) / std(YYt(:,i));
    end
end