% Saving and Plots for Generic Adaptive DO solvers

%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;
    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
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));
rhoi(pid,:)=fliplr(rhoi(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);
rhoi(pid,SIGNS<0)=-rhoi(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 = 5;
NPLOTSy = 1;
figure(1)
clf
PLTNUM=1;

subplot(NPLOTSx,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),
shading interp
caxis([min(-std(vort(:)),-eps),max(std(vort(:)),0)]*STDVORT+mean(vort(:)));
h=streamslice(XP,YP,tmp1,tmp2);
set(h,'color','k');
hold off
xlim(XLIM);ylim(YLIM)
axis equal
axis off
% TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');
% TT(3)=title(sprintf('Time =%03.03g \nw=[%0.2f,%0.2f],u=[%0.2f,%0.2f]',T+eps,...
%     min(dvdx(:)-dudy(:)),max(dvdx(:)-dudy(:)),...
%     min(sqrt(tmp1(:).^2+tmp2(:).^2)),max(sqrt(tmp1(:).^2+tmp2(:).^2))));
%formatplot([],gca,TT);view(2)
PLTNUM=PLTNUM+1;
subplot(NPLOTSx,NPLOTSy,PLTNUM);
pcolor(XP,YP,rho(NodeP(2:app.Ny+1,2:app.Nx+1))),shading interp,
xlim(XLIM);ylim(YLIM)
% TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');
% TT(3)=title(sprintf('Mean \nrho=[%0.2f,%0.2f]',...
%     min(rho(pid)),max(rho(pid))));
caxis([-std(rho(:)),std(rho(:))]*STDRHO+mean(rho(:)));
axis equal
axis off
% PLTNUM=PLTNUM+1;
% subplot(NPLOTSx,NPLOTSy,PLTNUM)
% if app.S(1)>0
%     semilogy(linspace(0,T,size(EIGPLOT,2)),EIGPLOT','.:');
% end
% hold on; semilogy(linspace(0,T,length(MEAN_ENER)),MEAN_ENER','--k'); hold off;
% xlim([0,max(T,eps)]);
% TT(1)=xlabel('Time(s)');TT(2)=ylabel('Var(Y_i)');
% TT(3)=title(sprintf('Variance'));
% 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,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;
    pcolor(XU(2:end,:),YV(:,2:end),vort),shading interp,
    caxis([-std(vort(:)),std(vort(:))]*STDVORT+mean(vort(:)));
    h=streamslice(XP,YP,tmp1,tmp2);
    set(h,'color','k');
    hold off
    xlim(XLIM);ylim(YLIM)
    axis equal
    axis off
    
%     TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');
%     TT(3)=title(sprintf('Mode =%g \nw=[%0.2f,%0.2f],u=[%0.2f,%0.2f]',i,...
%         min(dvdx(:)-dudy(:)),max(dvdx(:)-dudy(:)),...
%         min(sqrt(tmp1(:).^2+tmp2(:).^2)),max(sqrt(tmp1(:).^2+tmp2(:).^2))));
    %formatplot([],gca,TT);view(2)
    PLTNUM=PLTNUM+1;
    subplot(NPLOTSx,NPLOTSy,PLTNUM);
    tmp = rhoi(:,i);
    pcolor(XP,YP,tmp(NodeP(2:end-1,2:end-1))),shading interp,
%     caxis([]);
    xlim(XLIM);ylim(YLIM)
    axis equal
    axis off
%     TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');
%     TT(3)=title(sprintf('Mode =%g \nrho=[%0.2f,%0.2f]',i,...
%         min(rhoi(pid,i)),max(rhoi(:,i))));
    caxis([-std(rhoi(:,i)),std(rhoi(:,i))+eps]*STDRHO+mean(rhoi(:,i)));
    
    PLTNUM=PLTNUM+1;
    subplot(NPLOTSx,NPLOTSy,PLTNUM)
    %     [n x]=hist(YYt(:,i),linspace(min(YYt(:,i)),max(YYt(:,i)),100));
    [n x w] = ksdensity(YYt(:,i));
    [n x w] = ksdensity(YYt(:,i),'width',w/3);
    xnorm = (max(x)-min(x))/2; 
    plot((x-min(x))/xnorm - 1,n/max(n),'-');
    [n x] = hist(YYt(:,i),min(app.MC/100,50));
    hold on 
    bar((x-min(x))/xnorm - max((x-min(x))/xnorm)/2,n/max(n),'g');
    xlim([-1,1]);
    hold off
%     title(sprintf('Mode %d Var = %0.2f',i,EIGPLOT(i,end)));
    axis square
    axis off
    %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
% fprintf('\n');
% while PLTNUM<NPLOTSx*NPLOTSy
%     if k == 0
%        [dum ndex] = sort(-abs(YYt(:,1)));
%        ndex = ndex(1:NPLOTSx*NPLOTSy);
%     end
%     PLTNUM = PLTNUM+1;
%     subplot(NPLOTSx,NPLOTSy,PLTNUM);
%     tmp = rhoi*YYt(ndex(PLTNUM),:)'+rho;
%     pcolor(XP,YP,tmp(NodeP(2:app.Nx+1,2:app.Ny+1))),shading interp,
%     hold on 
%     uI = ui*YYt(ndex(PLTNUM),:)'+u;
%     vI = vi*YYt(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)));
%     h=streamslice(XP,YP,tmp1,tmp2);
%     set(h,'color','k');
%     %     axis equal
%     xlim(XLIM);ylim(YLIM)
%     TT(3)=title(sprintf('Realiz YYt=%0.2f \nrho=[%0.2f,%0.2f]',sum(abs(YYt(ndex(PLTNUM),:))),...
%         min(tmp),max(tmp)));
%     caxis([-std(tmp),std(tmp)+eps]*STDRHO+mean(tmp));
% end

GenericSaveScriptRho;