% 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
    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
     if isfield(app,'k')
        k = app.k;
    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.
% keyboard
%%
CYY=cov(YYt);
[VC,DC] = eig(CYY);  %Diagonal Covariance, and Diagonal Vectors
DC=diag(DC);
%Never flip
% if T>0.5
% the fliplr ensures that the first mode is the most important
if 0%k> 0;
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);
end

%%
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 = 4;
NPLOTSy = 4;
figure(1)
clf
% keyboard
%%
offu=size(ui,1);
offv=size(vi,1);
offP=size(Pi,1);
for i=1:app.S(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
    tmp = rhoi(offP*(i-1) + NodeP(2:app.Ny+1,2:app.Nx+1));
    pcolor(XP,YP,tmp),
    shading interp
 
    caxis([-std(tmp(:)),std(tmp(:))+eps]*STDRHO+mean(tmp(:)));
%     caxis([-0.5, 0.5] * app.drho);
    colorbar
    h=streamslice(XP,YP,tmp1,tmp2);
    set(h,'color','k');
    hold off
    axis equal
    xlim(XLIM);ylim(YLIM)
    axis off
    TT = title(sprintf('Mode %g',i));
    %formatplot([],gca,TT);view(2)
end
for i = 0:length(app.dref)-1
    %Deterministic run
    subplot(NPLOTSx,NPLOTSy,9 + i)
    det = load(sprintf('%s/Save/%s%d/%05d.mat',d0,app.ref,app.dref(i + 1),...
        k/app.PlotIntrvl),'rho','u','v');
    tmp = det.rho;
    pcolor(XP,YP,tmp(NodeP(2:app.Nx+1,2:app.Ny+1))),shading interp,
    hold on 
    tmp1 = det.u;
    tmp2 = det.v;

    tmp1 = 0.5*(tmp1(Nodeu(2:app.Ny+1,1:app.Nx  ))+tmp1(Nodeu(2:app.Ny+1,2:app.Nx+1)));
    tmp2 = 0.5*(tmp2(Nodev(2:app.Ny+1,2:app.Nx+1))+tmp2(Nodev(1:app.Ny,2:app.Nx+1)));
    h=streamslice(XP,YP,tmp1,tmp2);
    set(h,'color',[0.9 0.9 0.9]);
    axis equal
    xlim(XLIM);ylim(YLIM)
    axis off
    caxis([-0.5, 0.5]);
    TT = title(sprintf('Deterministic %g',i + 1));
%     keyboard 
%     if i == length(app.dref)-1
        colorbar
%     end
    %Realization
    
    %Find coefficients for actual projections
    YYT(i+1,:) =  (det.rho - rho)' * rhoi * dx * dy;
    
    subplot(NPLOTSx,NPLOTSy,5 + i)
    tmp = rhoi * YYt(i+1,:)' + rho;
%     tmp = rhoi * YYT(i+1,:)' + rho;
    pcolor(XP,YP,tmp(NodeP(2:app.Nx+1,2:app.Ny+1))),shading interp,
    hold on 
    tmp1 = ui*YYt(i+1,:)'+u;
    tmp2 = vi*YYt(i+1,:)'+v;

    tmp1 = 0.5*(tmp1(Nodeu(2:app.Ny+1,1:app.Nx  ))+tmp1(Nodeu(2:app.Ny+1,2:app.Nx+1)));
    tmp2 = 0.5*(tmp2(Nodev(2:app.Ny+1,2:app.Nx+1))+tmp2(Nodev(1:app.Ny,2:app.Nx+1)));
    h=streamslice(XP,YP,tmp1,tmp2);
    set(h,'color',[0.9 0.9 0.9]);
    axis equal
    xlim(XLIM);ylim(YLIM)
    axis off
    caxis([-0.5, 0.5]);
    TT = title(sprintf('Realization %g',i + 1));
    
%     if i == length(app.dref)-1
        colorbar
%     end
    
    %Difference
    subplot(NPLOTSx,NPLOTSy,13 + i)
    tmp = tmp - det.rho;
    tmp1 = ui*YYt(i+1,:)'+u - det.u;
    tmp2 = vi*YYt(i+1,:)'+v - det.v;
    nrm = sqrt((det.rho'*det.rho + det.u'*det.u + det.v'*det.v)*dx*dy);
    TT =  sqrt((tmp' * tmp + tmp1' * tmp1 + tmp2' * tmp2)*dx*dy) / nrm;
    tmp = tmp / nrm;
    pcolor(XP,YP,tmp(NodeP(2:app.Nx+1,2:app.Ny+1))),shading interp,
    
%     if i == length(app.dref) - 1
        colorbar%('SouthOutside')
%     end
%     caxis([-std(tmp(:)),std(tmp(:))+eps]*STDRHO*2+mean(tmp(:)));
    caxis([-1 1] * 0.01);
    hold on 
    

    tmp1 = 0.5*(tmp1(Nodeu(2:app.Ny+1,1:app.Nx  ))+tmp1(Nodeu(2:app.Ny+1,2:app.Nx+1)));
    tmp2 = 0.5*(tmp2(Nodev(2:app.Ny+1,2:app.Nx+1))+tmp2(Nodev(1:app.Ny,2:app.Nx+1)));
    h=streamslice(XP,YP,tmp1/nrm,tmp2/nrm);
    set(h,'color','k');
%     axis equal
    axis square
    xlim(XLIM);ylim(YLIM)
    axis off
    TT = title(sprintf('%0.1e', TT));  
end
% YYt(1:4,:)
% YYT
% YYt(1:4,:)-YYT
%%
    
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
pcolor(XP,YP,rho(NodeP(2:app.Ny+1,2:app.Nx+1))),
shading interp
caxis([-0.5, 0.5]);
colorbar

% LEVS = [0.49 0.4 0.3 0.2 0.1] * app.drho;
% contour(XP,YP,rho(NodeP(2:app.Ny+1,2:app.Nx+1)),LEVS,'k');
% contour(XP,YP,rho(NodeP(2:app.Ny+1,2:app.Nx+1)),-LEVS,'k--');
% contour(XP,YP,rho(NodeP(2:app.Ny+1,2:app.Nx+1)),[0 0],'w','LineWidth', 2);

h=streamslice(XP,YP,tmp1,tmp2);
set(h,'color',[0.9 0.9 0.9]);
% vort = sqrt(tmp1.^2+tmp2.^2);
% pcolor(XP,YP,vort),colorbar
% shading interp
% % caxis([-std(vort(:)),std(vort(:))]*STDVORT+mean(vort(:)));
% caxis([0,1]);
hold off
axis equal
xlim(XLIM);ylim(YLIM)
axis off
TT = title(sprintf('Time =%03.03g \nMean',T));

%%
GenericSaveScriptRho;
