% Saving and Plots for Generic Adaptive DO solvers
clf
%At first timestep do preliminaries
if k==0
    close all
    EIGPLOT=[];
    MEAN_ENER=[];
    
    %Make limits for plotting
    XLIM = [min(XP(:)),max(XP(:))];
    YLIM = [min(YP(:)),max(YP(:))];
    %Select the colormap
    colormap jet
    %How many std dev's away to plot vorticity
    STDVORT = 3;
    
    if length(app.S)>1 %Adaptation!
        Smin=app.S(1);
        Smax=app.S(2);
    else %No adaptation
        Smin=app.S;
        Smax=app.S;
    end
end

% 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(YYt);
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));
YYt=fliplr(YYt*(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);
YYt(:,SIGNS<0)=-YYt(:,SIGNS<0);

%%
EIGPLOT=[ EIGPLOT, [flipud(DC);NaN(Smax-length(DC),1)] ];
MEAN_ENER=[ MEAN_ENER, InProd(u(:),u(:),dx,dy,uid)+InProd(v(:),v(:),dx,dy,vid)];

%Plotting

NPLOTSx = 6;
NPLOTSy = 5;
figure(1)
clf
subplot(NPLOTSx,NPLOTSy,1:4)
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
tmp = u(1:Nbcs);
u(1:Nbcs) = NaN;
dudy = (u(Nodeu(2:end-2,2:end-1))-u(Nodeu(3:end-1,2:end-1)))/dy;
u(1:Nbcs) = tmp;
tmp = v(1:Nbcs);
v(1:Nbcs) = NaN;
dvdx = (v(Nodev(2:end-1,3:end-1))-v(Nodev(2:end-1,2:end-2)))/dx;
v(1:Nbcs) = tmp;
vort = dvdx-dudy;
% pcolor(XU(2:end,:),YV(:,2:end),vort),
% caxis([min(-std(vort(~isnan(vort))),-eps),max(std(vort(~isnan(vort))),0)]*STDVORT+mean(vort(~isnan(vort))));
% shading interp
tmp = P(1:Nbcs);
P(1:Nbcs) = NaN;
pcolor(XP,YP,P(NodeP(2:app.Ny +1, 2:app.Nx +1)));
shading interp
P(1:Nbcs) = tmp;
caxis([min(-std(P),-eps),max(std(P),0)]*STDVORT+mean(P));

h=streamslice(XP,YP,tmp1,tmp2);
set(h,'color','k');

% LEVU = linspace(0, max(P(pid)),10);
% LEVD = linspace(min(P(pid)),0,10);
% contour(XP,YP,P(NodeP(2:app.Ny+1,2:app.Nx+1)),LEVU,'k');
% contour(XP,YP,P(NodeP(2:app.Ny+1,2:app.Nx+1)),LEVD,'--k');

hold off
axis equal
xlim(XLIM);ylim(YLIM)
axis off
TT = title(sprintf('Time =%03.03g \nMean',T));

% subplot(NPLOTSx,NPLOTSy,26:30)
figure(2)
subplot(2,1,1)
if app.S(1)>0
    semilogy(linspace(0,T,size(EIGPLOT,2)),EIGPLOT','-','linewidth',1);
end
hold on;
semilogy(linspace(0,T,length(MEAN_ENER)),MEAN_ENER','--','linewidth',1);
hold off;
xlim([0,max(T,eps)]);
TT(1)=xlabel('Time');TT(2)=ylabel('Var(Y_i)');
TT(3)=title(sprintf('Variance'));
% legend('Mode1','Mode2','Mode3','Mode4','Mean','Location','SouthWest')
figure(1)
offu=size(ui,1);
offv=size(vi,1);
offP=size(Pi,1);
for i=1: min(5,app.S(1))
    subplot(NPLOTSx,NPLOTSy, 5 * i + 1: 5 * (i + 1) - 1)
    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
    tmp = ui(1:Nbcs,i);
    ui(1:Nbcs,i) = NaN;
    dudy = (ui((i-1)*offu+Nodeu(2:end-2,2:end-1))-ui((i-1)*offu+Nodeu(3:end-1,2:end-1)))/dy;
    ui(1:Nbcs,i) = tmp;
    tmp = vi(1:Nbcs,i);
    vi(1:Nbcs,i) = NaN;
    dvdx = (vi((i-1)*offv+Nodev(2:end-1,3:end-1))-vi((i-1)*offv+Nodev(2:end-1,2:end-2)))/dx;
    vi(1:Nbcs,i) = tmp;
%     vort = dvdx-dudy;
%     pcolor(XU(2:end,:),YV(:,2:end),vort),shading interp,
%     caxis([-std(vort(~isnan(vort))),std(vort(~isnan(vort)))]*STDVORT+mean(vort(~isnan(vort))));
    
        tmp = Pi(1:Nbcs,i);
        Pi(1:Nbcs,i) = NaN;
        pcolor(XP,YP,Pi(NodeP(2:app.Ny +1, 2:app.Nx +1)+offP*(i-1)));
        shading interp
        Pi(1:Nbcs,i) = tmp;
        caxis([min(-std(Pi(:,i)),-eps),max(std(Pi(:,i)),0)]*STDVORT+mean(Pi(:,i)));
    
        h=streamslice(XP,YP,tmp1,tmp2);
        set(h,'color','k');
    
%     if i < 3
%         LEVU = linspace(0, max(Pi(pid,i)),10);
%         LEVD = linspace(min(Pi(pid,i)),0,10);
%     else
%         LEVU = linspace(0, max(Pi(pid,i)),5);
%         LEVD = linspace(min(Pi(pid,i)),0,5);
%     end
%     contour(XP,YP,Pi(NodeP(2:app.Ny+1,2:app.Nx+1)+offP*(i-1)),LEVU,'k');
%     contour(XP,YP,Pi(NodeP(2:app.Ny+1,2:app.Nx+1)+offP*(i-1)),LEVD,'--k');
    
    hold off
    axis equal
    xlim(XLIM);ylim(YLIM)
    axis off
    TT = title(sprintf('\nMode %g',i));
    
    subplot(NPLOTSx,NPLOTSy, 5 * (i + 1))
    [n x w] = ksdensity(YYt(:,i));
    [n x] = ksdensity(YYt(:,i),'width',w/3);
    xnorm = (max(x)-min(x))/2;
    plot((x-min(x))/xnorm - 1,n/max(n),'-','linewidth',1);
    [n x] = hist(YYt(:,i),min(max(app.MC/100,sqrt(app.MC)),25));
    hold on
%     xnorm = (max(x)-min(x))/2;
%     bar((x-min(x))/xnorm - max((x-min(x))/xnorm)/2,n/max(n),'g');
    xlim([-1,1]);
%     title(sprintf('Marginal PDF\n Mode %d\n',i));
    hold off
%     axis off
    axis square
end
subplot(NPLOTSx,NPLOTSy,5)
% % % % title('Sample Scatter')
% % % plot3k(YYt(:,1:3),[],'.',3);
% % % axis([min(YYt(:,1)) max(YYt(:,1)) min(YYt(:,2)) max(YYt(:,2)) min(YYt(:,3)) max(YYt(:,3)) ])
% % % grid off
% % % % axis off
% % axis square
% % % view(2)
tmp = P(1:Nbcs);
P(1:Nbcs) = NaN;
pcolor(XP,YP,P(NodeP(2:app.Ny +1, 2:app.Nx +1)));
shading interp
P(1:Nbcs) = tmp;
caxis([min(-std(P),-eps),max(std(P),0)]*STDVORT+mean(P));
axis off
colorbar
%%
GenericSaveScript;
