% Plots for LidDrivenCavity with DO Modes.
%Save initially everything for plotting later
if k==0
    if isfield(app,'savemat')
        if isempty(app.savemat)
            d0=fileparts([pwd,filesep]);
            d0=d0(1:strfind(d0,'trunk')+4);
            save(sprintf('%s/Save/mat/param.mat',d0));
            fprintf('Saving %s mat-file in %s/Save/mat/parameters.mat',...
                SetupScript,d0);
        else
            if k < 10
                if ~exist(sprintf('%s/Save/%s/',d0,app.savemat),'dir')
                    mkdir(sprintf('%s/Save/%s/',d0,app.savemat))
                end
            end
            d0=fileparts([pwd,filesep]);
            d0=d0(1:strfind(d0,'trunk')+4);
            save(sprintf('%s/Save/%s/param.mat',d0,app.savemat));
            fprintf('Saving %s mat-file in %s/Save/%s/parameters.mat',...
                SetupScript,d0,app.savemat);
        end
    end
    colormap(mk_colmap)
    [a, maxids] = max(Yi);
    SIGNS1 = sign(diag(Yi(maxids,:)));
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(Yi);
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; viold=vi;
ui(uid,:)=fliplr(ui(uid,:)*(VC)); 
vi(vid,:)=fliplr(vi(vid,:)*(VC));
Pi(pid,:)=fliplr(Pi(pid,:)*(VC));
Yi=fliplr(Yi*(VC)); 
SIGNS = sign(diag(Yi(maxids,:))) .* SIGNS1;
SIGNS = [SIGNS, ones(length(SIGNS) - app.S, 1) ];
ui(uid,SIGNS<0)=-ui(uid,SIGNS<0);
vi(vid,SIGNS<0)=-vi(vid,SIGNS<0);
Pi(pid,SIGNS<0)=-Pi(pid,SIGNS<0);
Yi(:,SIGNS<0)=-Yi(:,SIGNS<0);
SIGNS1 = sign(diag(Yi(maxids,:)));

NPLOTSx=round(sqrt((app.S+1)*2));
NPLOTSy=ceil((app.S+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;
pcolor(XU(2:end,:),YV(:,2:end),dvdx-dudy),shading interp,caxis([-10,10])
h=streamslice(XP,YP,tmp1,tmp2);
set(h,'color','k');
hold off
xlim([0 1]);ylim([0 1])
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))));

PLTNUM=PLTNUM+1;
subplot(NPLOTSx,NPLOTSy,PLTNUM)
if k==0, EIGPLOT=[];end
EIGPLOT=[EIGPLOT, eig(CYY)];
if app.S>0
    semilogy(linspace(0,T,size(EIGPLOT,2)),EIGPLOT','.-');
end
xlim([0,max(T,eps)]);ylim([1e-3,1])
TT(1)=xlabel('Time(s)');TT(2)=ylabel('Var(Y_i)');
TT(3)=title(sprintf('Variance'));

offu=size(ui,1);
offv=size(vi,1);
offP=size(Pi,1);

for i=1:app.S
    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;
    pcolor(XU(2:end,:),YV(:,2:end),dvdx-dudy),shading interp, caxis([-30,30])
    h=streamslice(XP,YP,tmp1,tmp2);
    set(h,'color','k');
    hold off
    xlim([0 1]);ylim([0 1])
    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))));
        
    PLTNUM=PLTNUM+1;
    subplot(NPLOTSx,NPLOTSy,PLTNUM)
    [n x]=hist(Yi(:,i),linspace(min(Yi(:,i)),max(Yi(:,i)),100));
    plot(x,n/app.MC/(x(2)-x(1)),':.');
    xlim([-2,2]);
end


%Save image
if isfield(app,'saveimg')
    d0=fileparts([pwd,filesep]);
    d0=d0(1:strfind(d0,'trunk')+4);
    print('-djpeg85','-r200',...
        sprintf('%s/Save/img/%05d.jpg',d0,k/app.PlotIntrvl));
    fprintf('Saving %s image in %s/Save/img/%05d.jpg\n',...
             SetupScript,d0,k/app.PlotIntrvl);
end

%Save only essentials for plotting later
if isfield(app,'savemat')
    if isempty(app.savemat)
        d0=fileparts([pwd,filesep]);
        d0=d0(1:strfind(d0,'trunk')+4);
        save(sprintf('%s/Save/mat/%05d.mat',d0,k/app.PlotIntrvl),'P','u','v','T','rho');
        fprintf('Saving %s mat-file in %s/Save/mat/%05d.mat\n',...
             SetupScript,d0,k/app.PlotIntrvl);
    else
        if k < 10
            if ~exist(sprintf('%s/Save/%s/',d0,app.savemat),'dir')
                mkdir(sprintf('%s/Save/%s/',d0,app.savemat))
            end
        end
        d0=fileparts([pwd,filesep]);
        d0=d0(1:strfind(d0,'trunk')+4);
        save(sprintf('%s/Save/%s/%05d.mat',d0,app.savemat,k/app.PlotIntrvl),'P','u','v','T','rho');
        fprintf('Saving %s mat-file in %s/Save/%s/%05d.mat\n',...
             SetupScript,d0,app.savemat,k/app.PlotIntrvl);
    end
end

drawnow;
