% Modified by Jing Lin from DORhoAdaptPostPlots.m in Feb. 2014.

% Saving and Plots for Generic Adaptive DO solvers

%% Load Data

%At first timestep do preliminaries
if first_pass
    EIGPLOT=[];
    TimePlot = [];
    MEAN_ENER=[];
    
    %Make limits for plotting
    XLIM = [min(XP(:)), max(XP(:))];
    YLIM = [0.6, max(YP(:))];
    %Select the colormap
    %colormap jet
    colormap(mk_colmap)
    %How many std dev's away to plot vorticity
    STDVORT = 3;
    STDRHO = 2;
    if length(app.S)>1 %Adaptation!
        Smin=app.S(1);
        Smax=app.S(2);
    else %No adaptation
        Smin=app.S;
        Smax=app.S;
    end
    first_pass=0;
    was_first = 1;
end

load(sprintf('%s/%s/%s', d0, app.savemat, DataNum));
%load (fullfile(app.savemat,[num2str(k,'%03d'),'.mat']));

u = u*(nu_x/h); v = v*(nu_x/h); rho = rho*app.DelS;
ui = ui*(nu_x/h); vi = vi*(nu_x/h); rhoi = rhoi*app.DelS;

u(1:Nbcs) = NaN; v(1:Nbcs) = NaN; rho(1:Nbcs) = NaN;
ui(1:Nbcs,:) = NaN; vi(1:Nbcs,:) = NaN; rhoi(1:Nbcs,:) = NaN;

pid = Nbcs+1:length(P);
uid = Nbcs+1:length(u);
vid = Nbcs+1:length(v);


%% Orthonormalizaztion

% % First Step is to diagonalize the modes (for Plotting, should have been done
% % in DOAdapt, but in general we need to make sure.
% CYY=cov(Yi);
% clear DC
% [VC,DC] = eig(CYY);  %Diagonal Covariance, and Diagonal Vectors
% DC=diag(DC);
% 
% % the fliplr ensures that the first mode is the most important
% uiold=ui;
% ui(uid,:)=fliplr(ui(uid,:)*(VC));
% vi(vid,:)=fliplr(vi(vid,:)*(VC));
% Pi(pid,:)=fliplr(Pi(pid,:)*(VC));
% rhoi(pid,:)=fliplr(rhoi(pid,:)*(VC));
% Yi=fliplr(Yi*(VC));
% 
% %Keep the signs the same for plotting purposes:
% SIGNS = sign(sum(uiold.*ui));
% ui(uid,SIGNS<0)=-ui(uid,SIGNS<0);
% vi(vid,SIGNS<0)=-vi(vid,SIGNS<0);
% Pi(pid,SIGNS<0)=-Pi(pid,SIGNS<0);
% rhoi(pid,SIGNS<0)=-rhoi(pid,SIGNS<0);
% Yi(:,SIGNS<0)=-Yi(:,SIGNS<0);
% 
EIGPLOT=[EIGPLOT, var(Yi)'];
TimePlot = [TimePlot, k];
% MEAN_ENER=[ MEAN_ENER, InProd(u(:),u(:),dx,dy,uid)+InProd(v(:),v(:),dx,dy,vid)];


%% Plotting

clf;
Fig_ID = 1;
set(figure(Fig_ID), 'Position', [0, 100, 1200, 1400]);

NPLOTSx=(app.S(1)+1)*3;
NPLOTSy=3;

% Mean Velocity Plot
PLTNUM=1;
subplot(RowPerFig, NPLOTSy, PLTNUM)
tmp1 = 0.5*(u(Nodeu(2:app.Ny+1,1:app.Nx))+u(Nodeu(2:app.Ny+1,2:app.Nx+1)));
%tmp2 = 0.5*(v(Nodev(2:app.Ny+1,2:app.Nx+1))+v(Nodev(1:app.Ny,2:app.Nx+1)));
hold on
%dudy = (u(Nodeu(2:end-2,2:end-1))-u(Nodeu(3:end-1,2:end-1)))/dy;
%dvdx = (v(Nodev(2:end-1,3:end-1))-v(Nodev(2:end-1,2:end-2)))/dx;
%vort = dvdx-dudy;
%pcolor(XU(2:end,:),YV(:,2:end),vort);
pcolor(XU,YU,u(Nodeu(2:end-1,2:end-1)));
shading interp
%vort = vort(~isnan(vort));
%caxis([min(-std(vort(:)),-eps),max(std(vort(:)),0)]*STDVORT+mean(vort(:)));
%Fig_Hd=streamslice(XP,YP,tmp1,tmp2);
%set(Fig_Hd,'color','k');
hold off
xlim(XLIM); ylim(YLIM);
TT(3)=title(sprintf('Time =%03.03gh, Mean u, u=[%0.1e,%0.1e]',k/60, min(tmp1(pid)),max(tmp1(pid))));
% TT(3)=title(sprintf('Time =%03.03g \nw=[%0.2f,%0.2f],u=[%0.2f,%0.2f]',T+eps,...
%     min(dvdx(pid)-dudy(pid)),max(dvdx(pid)-dudy(pid)),...
%     min(sqrt(tmp1(pid).^2+tmp2(pid).^2)),max(sqrt(tmp1(pid).^2+tmp2(pid).^2))));

% ridgex = [0 33.2667 33.2667 7.414 0].*(h/1000);
% ridgey = [0 0 4.25 5.37 5.37].*(h/1000);
% ridgeColor = [105 105 105]/255;
% patch(ridgex,ridgey,ridgeColor);

% Mean Rho Plot
PLTNUM = PLTNUM+1;
subplot(RowPerFig, NPLOTSy, PLTNUM);
pcolor(XP,YP,rho(NodeP(2:app.Ny+1,2:app.Nx+1))),shading interp,
xlim(XLIM); ylim(YLIM)
TT(3)=title(sprintf('Mean rho, rho=[%0.1e,%0.1e]',...
    min(rho(pid)),max(rho(pid))));
%caxis([-std(rho(pid)),std(rho(pid))]*STDRHO+mean(rho(pid)));

% Mode Variance Evolution Plot
PLTNUM = PLTNUM+1;
subplot(RowPerFig, NPLOTSy, PLTNUM)
if app.S(1)>0
    semilogy(TimePlot, EIGPLOT', '.:');
end;
% hold on; semilogy(linspace(0,T,length(MEAN_ENER)),MEAN_ENER','--k'); hold off;
%xlim([0,max(T,eps)]);
TT(3)=title(sprintf('Variance'));

% Plot modes
offu=size(ui,1);
offv=size(vi,1);
offP=size(rhoi,1);
for i = 1 : app.S(1)
    % Velocity mode plot
    PLTNUM = PLTNUM+1;
    subplot(RowPerFig, NPLOTSy, PLTNUM)
    tmp1 = 0.5*(ui((i-1)*offu+Nodeu(2:app.Ny+1,1:app.Nx))+ui((i-1)*offu+Nodeu(2:app.Ny+1,2:app.Nx+1)));
    %tmp2 = 0.5*(vi((i-1)*offv+Nodev(2:app.Ny+1,2:app.Nx+1))+vi((i-1)*offv+Nodev(1:app.Ny,2:app.Nx+1)));
    hold on
    %dudy = (ui((i-1)*offu+Nodeu(2:end-2,2:end-1))-ui((i-1)*offu+Nodeu(3:end-1,2:end-1)))/dy;
    %dvdx = (vi((i-1)*offv+Nodev(2:end-1,3:end-1))-vi((i-1)*offv+Nodev(2:end-1,2:end-2)))/dx;
    %vort = dvdx-dudy;
    %tmp2 = ui(:,i);
    pcolor(XP,YP,tmp1);
    shading interp
    %vort = vort(~isnan(vort));
    %caxis([min(-std(vort(:)),-eps),max(std(vort(:)),0)]*STDVORT+mean(vort(:)));
    %Fig_Hd=streamslice(XP,YP,tmp1,tmp2);
    %set(Fig_Hd,'color','k');
    hold off
    xlim(XLIM); ylim(YLIM)
    TT(3)=title(sprintf('Mode %d, u=[%0.1e,%0.1e]',i,min(ui(uid,i)),max(ui(uid,i))));
    
    % Rho mode plot
    PLTNUM = PLTNUM+1;
    subplot(RowPerFig, NPLOTSy, PLTNUM);
    tmp = rhoi(:,i);
    pcolor(XP,YP,tmp(NodeP(2:end-1,2:end-1))),shading interp,
    xlim(XLIM);ylim(YLIM)
    TT(3)=title(sprintf('Mode %d, rho=[%0.1e,%0.1e]',i,...
        min(rhoi(pid,i)),max(rhoi(pid,i))));
    %caxis([-std(rhoi(pid,i)),std(rhoi(pid,i))+eps]*STDRHO+mean(rhoi(pid,i)));
    
    % Coeffcient distribution plot
    PLTNUM = PLTNUM+1;
    subplot(RowPerFig, NPLOTSy, PLTNUM)
    [~, ~, w] = ksdensity(Yi(:,i));
    [n, x] = ksdensity(Yi(:,i),'width',w/2);
    plot(x, n,'-','linewidth',1);
    [n2, x2] = hist(Yi(:,i),min(max(app.MC/100,sqrt(app.MC)),25));
    hold on
    bar(x2,n2/max(n2)*max(n),'g');
    title(sprintf('Mode %d, Var = %0.1e',i,EIGPLOT(i,end)));
    hold off
    
    if (mod(PLTNUM/NPLOTSy, RowPerFig)==0)
        drawnow;
        set(figure(Fig_ID), 'PaperPosition', [0, 0, 12, 14]);
        if (Fig_ID == 1)
            figname = sprintf('%s/%s/Mean_Mode1to%d_%s',...
                d0, app.saveimg, RowPerFig-1, DataNum);
        else
            figname = sprintf('%s/%s/Mode%dto%d_%s', d0, app.saveimg, ...
                (Fig_ID-1)*RowPerFig, Fig_ID*RowPerFig - 1, DataNum);
        end;
        print('-dpng', '-r400', figname);
        Fig_ID = Fig_ID + 1;
        set(figure(Fig_ID), 'Position', [0, 100, 1200, 1400]);
        PLTNUM = 0;
        %colormap jet
        colormap(mk_colmap)
    end;
end

% Plot realizations to fill blank positions left.
flag = false;
while (PLTNUM > 0)&&(PLTNUM < RowPerFig*NPLOTSy)
    flag = true;
    if was_first
        [~, ndex] = sort(-abs(Yi(:,1)));
        ndex = ndex(1:RowPerFig*NPLOTSy);
        was_first = 0;
    end;
    PLTNUM = PLTNUM+1;
    subplot(RowPerFig, NPLOTSy, PLTNUM);
    tmp = rhoi*Yi(ndex(PLTNUM),:)'+rho;
    pcolor(XP,YP,tmp(NodeP(2:app.Ny+1,2:app.Nx+1))),shading interp,
    hold on
    uI = ui*Yi(ndex(PLTNUM),:)'+u;
    vI = vi*Yi(ndex(PLTNUM),:)'+v;
    
    %tmp1 = 0.5*(uI(Nodeu(2:app.Ny+1,1:app.Nx  ))+uI(Nodeu(2:app.Ny+1,2:app.Nx+1)));
    %tmp2 = 0.5*(vI(Nodev(2:app.Ny+1,2:app.Nx+1))+vI(Nodev(1:app.Ny,2:app.Nx+1)));
    %Fig_Hd = streamslice(XP,YP,tmp1,tmp2);
    %set(Fig_Hd,'color','k');
    
    xlim(XLIM);ylim(YLIM)
    tmp = tmp(pid);
    TT(3)=title(sprintf('Realiz Yi=%0.1e, rho=[%0.1e,%0.1e]',sum(abs(Yi(ndex(PLTNUM),:))),...
        min(tmp),max(tmp)));
    %caxis([-std(tmp),std(tmp)+eps]*STDRHO+mean(tmp));
end;

if flag
    drawnow;
    set(figure(Fig_ID), 'PaperPosition', [0, 0, 12, 14]);
    figname = sprintf('%s/%s/Mode%dto%d_%s', d0, app.saveimg, ...
        (Fig_ID-1)*RowPerFig, app.S(1), DataNum);
    print('-dpng', '-r400', figname);
end;


