% Saving and Plots for Generic Adaptive DO solvers
clf
%At first timestep do preliminaries
if k==0
    close all
    EIGPLOT=[];
    MEAN_ENER=[];
    
    [a, maxids] = max(YYt);
    SIGNS1 = sign(diag(YYt(maxids,:)));
    %Make limits for plotting
    XLIM = [min(XP(:)),max(XP(:))];
    YLIM = [min(YP(:)),max(YP(:))];
    %Select the colormap
    %colormap jet
    colormap(mk_colmap)
    %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(YYt(maxids,:));
SIGNS1 = [SIGNS1, ones(app.S - length(SIGNS1), 1) ];
SIGNS1 = SIGNS1(1:length(SIGNS));
iidss = SIGNS(:)~=SIGNS1(:);
ui(uid, iidss)=-ui(uid, iidss);
vi(vid, iidss)=-vi(vid, iidss);
Pi(pid, iidss)=-Pi(pid, iidss);
YYt(:, iidss)=-YYt(:, iidss);
SIGNS1 = sign(YYt(maxids,:));

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=round(sqrt((app.S(1)+1)*2));
NPLOTSy=ceil((app.S(1)+1)*2/NPLOTSx);
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
% axis equal
xlim(XLIM);ylim(YLIM)
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');
TT(3)=title(sprintf('Time =%03.03g \nw=[%0.2g,%0.2g],u=[%0.2g,%0.2g]',T,...
    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)
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','o-'); 
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
    %     axis equal
    xlim(XLIM);ylim(YLIM)
    TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');
    TT(3)=title(sprintf('Mode =%g \nw=[%0.2g,%0.2g],u=[%0.2g,%0.2g]',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)
        [n x w] = ksdensity(YYt(:,i));
    [n x] = ksdensity(YYt(:,i),'width',w/2);
    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
    bar((x-min(x))/xnorm - max((x-min(x))/xnorm)/2,n/max(n),'g');
    xlim([-1,1]);
    title(sprintf('Mode %d Var = %g',i,EIGPLOT(i,end)));
    hold 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
%keyboard
GenericSaveScript;