% 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
    %For saving
    d0=fileparts([pwd,filesep]);
	d0=d0(1:strfind(d0,'trunk')+4);
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:
%This is annoying, but for the reference solution, the sign is -0.00000000
%for the middle one, so we have to adjust for it. (This is because two
%modes were switched!) --> Have to adjust for this in the general case
%later
if k == 0 
    SIGNS = [1 -1 1];
else
    SIGNS = sign([sum(uiold.*ui)]);
end
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)];

Pc = P;
Pic = Pi;
uic = ui;
vic = vi;
uc = u;
vc = v;
rhoc = rho;
YYtc =  YYt;
MEAN_ENERc = MEAN_ENER;
EIGPLOTc = EIGPLOT;
%load the reference solution
RSOL = load(sprintf('%s/Save/%s/%05d.mat',d0,app.ref,k/app.PlotIntrvl),...
    'P','Pi','rho','u','ui','v','vi','YYt','MEAN_ENER','EIGPLOT');

%Interpolate unto the coarser mesh and do differences
Nxe = size(REF.XP,2);
Nye = size(REF.XP,1);
INTERPTYPE = 'cubic';
uiold = interp2(REF.XP, REF.YP, RSOL.P(REF.NodeP(2:Nye +1, 2:Nxe+1)), XP, YP, INTERPTYPE); 
P = Pc - mean(Pc(idsP)) - [bcsDP'; bcsNP'; uiold(idsP) - mean(uiold(idsP))];
%P = Pc - [bcsDP'; bcsNP'; uiold(idsP)];

uiold = interp2(REF.XP, REF.YP, RSOL.rho(REF.NodeP(2:Nye +1, 2:Nxe+1)), XP, YP, INTERPTYPE); 
rho = rhoc - [bcsDrho'; bcsNrho'; uiold(idsP)];

uiold = interp2(REF.XU, REF.YU, RSOL.u(REF.Nodeu(2:Nye+1, 2:Nxe)), XU, YU, INTERPTYPE); 
u = uc - [bcsDu'; bcsNu'; uiold(idsu)];

uiold = interp2(REF.XV, REF.YV, RSOL.v(REF.Nodev(2:Nye, 2:Nxe+1)), XV, YV, INTERPTYPE); 
v = vc - [bcsDv'; bcsNv'; uiold(idsv)];

tmp = zeros(size(Pic));
for i = 1:app.S(1)
    uiold = RSOL.Pi(:,i);
    uiold = interp2(REF.XP, REF.YP, uiold(REF.NodeP(2:Nye +1, 2:Nxe+1)), XP, YP, INTERPTYPE); 
    tmp(:,i) = abs(Pic(:,i)) - abs([0*bcsDP'; 0*bcsNP'; uiold(idsP)]);
end    
Pi = tmp;
tmp = zeros(size(uic));
for i = 1:app.S(1)
    uiold = RSOL.ui(:,i);
    uiold = interp2(REF.XU, REF.YU, uiold(REF.Nodeu(2:Nye +1, 2:Nxe)), XU, YU, INTERPTYPE); 
    tmp(:,i) = abs(uic(:,i)) - abs([0*bcsDu'; 0*bcsNu';  uiold(idsu)]);
end 
ui = tmp;
tmp = zeros(size(vic));
for i = 1:app.S(1)   
    uiold = RSOL.vi(:,i);
    uiold = interp2(REF.XV, REF.YV, uiold(REF.Nodev(2:Nye , 2:Nxe+1)), XV, YV, INTERPTYPE); 
    tmp(:,i) = abs(vic(:,i)) - abs([0*bcsDv'; 0*bcsNv'; uiold(idsv)]);
end
vi = tmp;

YYt = (abs(YYtc) - abs(RSOL.YYt));
MEAN_ENER = abs(abs(MEAN_ENERc) - abs(RSOL.MEAN_ENER));
EIGPLOT = abs(abs(EIGPLOTc) - abs(RSOL.EIGPLOT));

%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),
vort = sqrt(tmp1.^2+tmp2.^2);
pcolor(XP,YP,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,
    vort = sqrt(tmp1.^2+tmp2.^2);
    pcolor(XP,YP,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]=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),':.');
    [n x] = hist(YYt(:,i),app.MC/100);
    hold on 
    bar(x/xnorm,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

%% Calculate error norms and save to file

PL2(k/app.PlotIntrvl+1) = P'*P*dx*dy;
PiL2(:,k/app.PlotIntrvl+1) = diag(Pi'*Pi*dx*dy);
uL2(k/app.PlotIntrvl+1) = u'*u*dx*dy;
uiL2(:,k/app.PlotIntrvl+1) = diag(ui'*ui*dx*dy);
vL2(k/app.PlotIntrvl+1) = v'*v*dx*dy;
viL2(:,k/app.PlotIntrvl+1) = diag(vi'*vi*dx*dy);
rhoL2(k/app.PlotIntrvl+1) = rho'*rho*dx*dy;
MEAN_ENERL2 = MEAN_ENER.^2;
EIGPLOTL2 = EIGPLOT.^2;
YYtL2(:,k/app.PlotIntrvl+1) = diag(YYt'*YYt / app.MC);

PLinf(k/app.PlotIntrvl+1) = max(abs(P));
PiLinf(:,k/app.PlotIntrvl+1) = max(abs(Pi),[],1)';
uLinf(k/app.PlotIntrvl+1) = max(abs(u));
uiLinf(:,k/app.PlotIntrvl+1) = max(abs(ui),[],1)';
vLinf(k/app.PlotIntrvl+1) = max(abs(v));
viLinf(:,k/app.PlotIntrvl+1) = max(abs(vi),[],1)';
rhoLinf(k/app.PlotIntrvl+1) = max(abs(rho));
MEAN_ENERLinf = (abs(MEAN_ENER));
EIGPLOTLinf = abs(EIGPLOT);
YYtLinf(:,k/app.PlotIntrvl+1) = max(abs(YYt),[],1)';

if k == 0
    if ~exist(sprintf('%s/Save/%s/',d0,app.savemat),'dir')
        mkdir(sprintf('%s/Save/%s/',d0,app.savemat));
    end
end

 save(sprintf('%s/Save/%s/RE%d_NX%03d_dt%03d.mat',d0,app.savemat,1/app.nu,app.Nx,1/app.dt),...
     'PL2','PiL2','uL2','uiL2','vL2','viL2','rhoL2','MEAN_ENERL2','EIGPLOTL2','YYtL2',...
     'PLinf','PiLinf','uLinf','uiLinf','vLinf','viLinf','rhoLinf','MEAN_ENERLinf','EIGPLOTLinf','YYtLinf');
fprintf('%05d: Saving %s/Save/%s/RE%d_NX%03d_dt%03d.mat\n',k,d0,app.savemat,1/app.nu,app.Nx,1/app.dt); 
%Restore variable so that calculatation may be continued
P = Pc;
Pi = Pic;
ui = uic;
vi = vic;
u = uc;
v = vc;
rho = rhoc;
YYt =  YYtc;
MEAN_ENER = MEAN_ENERc;
EIGPLOT = EIGPLOTc;

% figure(2),
% clf
% MEAN_ENER = MEAN_ENER(1:end-1);
% EIGPLOT = EIGPLOT(:,1:end-1);
% DOAdaptGenericPlots;
