% 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 = 3;
NPLOTSy = 4;

if isfield(app, 'saveimg')
    figure(1)
    clf
    PLTNUM=1;

    subplot(NPLOTSx,NPLOTSy,1)
    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),
    % shading interp
    % caxis([min(-std(vort(:)),-eps),max(std(vort(:)),0)]*STDVORT+mean(vort(:))
    % );
    pcolor(XP, YP, P(NodeP(2:app.Ny+1,2:app.Nx+1)))
    shading interp
%     colorbar
    % caxis([min(-std(P(:)),-eps),max(std(P(:)),eps)]*STDVORT+mean(P(:)));
    caxis([-0.7 0.7]);

    %h=streamslice(XP,YP,tmp1,tmp2);
    %set(h,'color',[1 1 1]*0.3);
    % LEVU = linspace(0, max(P(pid)),20);
    % LEVD = linspace(min(P(pid)),0,20);
    % 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));
    %formatplot([],gca,TT);view(2)

    PLTNUM=PLTNUM+1;
    subplot(NPLOTSx,NPLOTSy,9:12)
    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','Mean','Location','SouthWest')
    % formatplot([],gca,TT);

    offu=size(ui,1);
    offv=size(vi,1);
    offP=size(Pi,1);
    for i=1:app.S(1)
        PLTNUM=PLTNUM+1;
        subplot(NPLOTSx,NPLOTSy,1+i)
        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;
    %     pcolor(XU(2:end,:),YV(:,2:end),vort),shading interp,
    %     caxis([-std(vort(:)),std(vort(:))]*STDVORT+mean(vort(:)));
        pcolor(XP, YP, Pi(NodeP(2:app.Ny+1,2:app.Nx+1) + offP*(i-1)))
        shading interp
        caxis([min(-std(Pi(:,i)),-eps),max(std(Pi(:,i)),eps)]*STDVORT+mean(Pi(:,i)));

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

    %     if i < 4
    %         LEVU = linspace(0, max(Pi(pid,i)),5);
    %         LEVD = linspace(min(Pi(pid,i)),0,5);
    %     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));
        %formatplot([],gca,TT);view(2)

        PLTNUM=PLTNUM+1;
        subplot(NPLOTSx,NPLOTSy,5 + i)
        %     [n x]=hist(YYt(:,i),linspace(min(YYt(:,i)),max(YYt(:,i)),100));
        [n x] = ksdensity(YYt(:,i));
        xnorm = (max(x)-min(x))/2; 
        plot(x/xnorm,n/max(n),'-','linewidth',1);
        [n x] = hist(YYt(:,i),min(app.MC/100,25));
    %     hold on 
    %     bar(x/xnorm,n/max(n),'g');
        xlim([-1,1]);
        title(sprintf('Marginal PDF\n Mode %d\n',i));
        hold off
        axis square
        %for j=1:app.S(1)
        %    coef(j) = InProd(ui(:,i),ui(:,j),dx,dy,uid)+InProd(vi(:,i),vi(:,j),dx,dy,vid);
        %end
        %fprintf('Mode %g Inner products:  %s\n',i,sprintf('%4.1g\t',coef));
    end
    subplot(NPLOTSx,NPLOTSy,5)
    title('Sample Scatter')
    plot3k([YYt;10 10 1; 10 10 -1],[],'.',3);
    % axis([min(YYt(:,1)) max(YYt(:,1)) min(YYt(:,2)) max(YYt(:,2))
    % min(YYt(:,3)) max(YYt(:,3)) ])
    axis([-2 2 -1.2 1.2])

    grid off
    % colorbar
    % xlabel('Y_1')
    % ylabel('Y_2')
    axis square
    view(2)
end
%%
GenericSaveScript;
