clear all, clc, clf ,close all
d0=fileparts([pwd,filesep]);
d0=d0(1:strfind(d0,'trunk')+4);

%Set Re and MC
Re = 40;
app.MC = 100000;
app.DT_MUL = 1;
app.N = 15
path = ['Save/randombcs/BConlyCylinder_RE',num2str(Re),'_MC',num2str(app.MC/100),...
    '_NX',num2str(app.N),'_O'];
% path = ['Save/randombcs/Cylinder_RE',num2str(Re),'_MC',num2str(app.MC/100),...
%     '_NX',num2str(app.N),'_O'];
%load param file
load(sprintf('%s/%s/param.mat', d0, path));
%load final time
d = dir(sprintf('%s/%s/0*.mat', d0, path));
sol = load(sprintf('%s/%s/%s', d0, path, d(369).name));
[sol.YYt sol.ui sol.vi sol.Pi] = DOorthnorm(app, dx, dy, sol.YYt, sol.ui, uid, sol.vi, vid, sol.Pi, pid);
fprintf('%s/%s/%s\n', d0, path, d(end).name)
app.DT_MUL = 1;
%Sort the first columns of YYt
YYtsort = sortrows(sol.YYt, 1);

%% Plot the time evolution of the stochastic energy
% subplot(NPLOTSx,NPLOTSy,26:30)
figure
%subplot(2,1,1)
Tstart = (sol.k - size(sol.EIGPLOT,2) * app.PlotIntrvl) * app.dt;
if app.S(1)>0
    semilogy(linspace(Tstart,sol.T,size(sol.EIGPLOT,2)),sol.EIGPLOT','-','linewidth',1);
end
hold on;
semilogy(linspace(Tstart,sol.T,length(sol.MEAN_ENER)),sol.MEAN_ENER','--','linewidth',1);
hold off;
xlim([Tstart,max(sol.T,eps)]);
TT(1)=xlabel('Time');TT(2)=ylabel('Var(Y_i)');
TT(3)=title(sprintf('Variance'));
legend('Mode1','Mode2','Mode3','Mode4','Mean','Location','SouthWest')
ylim([10^-6, 100])
print('-dpng','-r400',sprintf('%s/%s/Fnl_yyt_evol', d0, path));



%% Plot the marginals vs RE
YYtsort = sortrows(sol.YYt, 1);
RE_locs = (YYtsort(:, 1) * sol.ui(2,1) + sol.u(2)) / app.nu;
%We have to create a grid to, well, grid the data, then do interpolation
%To know how large grid should be, we need to know the extreme values of
%all the pdfs
extremes = [min(YYt); max(YYt)];
%Whether or not to normalize the PDF at each reynolds number
normalize = true
mkY = true
logscale = false
plotdt = false
Stot = app.S
darkplot = false
yoff = 1e-2 / 5;

%Now divide the domain:
Nx = 350 + 1;
%Ny = 500;
Ny = round(2.5*sqrt(app.MC / Nx));

rex = linspace(min(RE_locs), max(RE_locs), Nx);
rey = zeros(Ny, app.S-1);
traceCyy = zeros(Nx, 1);
for j = 2:app.S
    rey(:, j-1) = linspace(extremes(1,j), extremes(2, j), Ny)';
end

%Now calculate rez
rez = cell(app.S-1, 1);
if mkY
    Y = cell(app.S-1, 1);
end
for j = 2:app.S
    rez{j-1} = zeros(Ny, Nx-1);
    if mkY
        Y{j-1} = zeros(Ny, Nx-1);
    end
end

if plotdt
    %YYtsort is assigned in the following script
   replot_getRHS;
   YYtsort = RHSY; 
   YYtsort = [sol.YYt(:, 1), YYtsort(:,2:end)];
   YYtsort = sortrows(YYtsort, 1);
end 
    

for i = 1:length(rex)-1
    fprintf('Calculating slice %d of %d\n', i, Nx-1)
    ids = find((RE_locs >= rex(i)) & (RE_locs < rex(i+1)));    
    YYtmean = mean(YYtsort(ids, :), 1);  
    for j = 2:app.S
        %[n x w] = ksdensity(YYtsort(ids, j));
        %[n x] = ksdensity(YYtsort(ids, j),'width',w/3);
        if logscale
            if mkY
                [n x] = hist(abs(YYtsort(ids, j)),Ny);
            else
                [n x] = hist(abs(YYtsort(ids, j)),sqrt(app.MC/Ny));
            end

        else
            if mkY
                [n x] = hist(YYtsort(ids, j),Ny);
            else
                [n x] = hist(YYtsort(ids, j),sqrt(app.MC/Ny));
            end
        end
        
        if max(x) > max (YYtsort(:, j))
            fprintf('Weird, max x is %g, while max YYt is %g\n', max(x), max(YYtsort(:, j)))
        end
        if normalize
            n = n / max(n);
        end

        if mkY
            %Y{j-1}(:, i) = linspace(min(x), max(x), Ny)';
            %rez{j-1}(:, i) = interp1(x, n, Y{j-1}(:, i));
            Y{j-1}(:, i) = x;
            rez{j-1}(:, i) = n;
        else
            rez{j-1}(:, i) = interp1(x, n, rey(:, j-1));            
        end                
        if darkplot
            rez{j-1}(:,i) = 0;
        end
    end
    traceCyy(i) = trace(cov(YYtsort(ids, :)));
end
%%    
figure
for Stot = 2:2:app.S
    clf
    Stot
    
    % Now we plot these suckers

    rex2 = (rex(1:end-1) + rex(2:end)) / 2;
    for j = 1:Stot-1
        if logscale
            subplot(Stot-1, 1, j, 'Color',[0 0 0],'YScale','log')
            %subplot(app.S-1, 1, j, 'Color',[0.25 0 0.25])
        else
            subplot(Stot-1, 1, j, 'Color',[0 0 0])
        end
        hold on
        colormap(mk_colmap)
        if ~mkY
            [X,yy] =meshgrid(rex2, rey(:, j));
        else
            [X] =meshgrid(rex2, rey(:, j));
            yy = Y{j};
        end
        yy(1, abs(yy(1, :) - yy(end,:)) < yoff) = yy(1, abs(yy(1, :) - yy(end,:)) < yoff) - yoff;
        
        pcolor(X, yy, rez{j})
        if plotdt
            plot([rex2(1), rex2(end)], [0,0], '--', 'Color',[0.25 0 0.25])
        end
        xlim([min(rex2), max(rex2)])
        ylim([min(yy(:)), max(yy(:))])
        hold off
        shading flat
    end
    if plotdt
        print('-dpng','-r400',sprintf('%s/%s/Fnl_dydt_v_Re_%d', d0, path, Stot));
    else
        print('-dpng','-r400',sprintf('%s/%s/Fnl_yyt_v_Re_%d', d0, path, Stot));
    end
end
%Sort the first columns of YYt
YYtsort = sortrows(sol.YYt, 1);
%%
close all
semilogy(rex, traceCyy, 'linewidth',2)
h=xlabel('Re');
%set(h,'FontSize',16)
xlim([min(rex2), max(rex2)])
% ylim([1e-5,1e-1])
if plotdt
    h=ylabel('Trace(C_{(dydt),y})')
    %set(h,'FontSize',16);
    print('-dpng','-r400',sprintf('%s/%s/Fnl_tracedCyydt_v_Re_%d', d0, path, Stot));
    print('-deps',sprintf('%s/%s/Fnl_tracedCyydt_v_Re_%d', d0, path, Stot));
else
    ylabel('Trace(C_{y,y})')
    print('-dpng','-r400',sprintf('%s/%s/Fnl_tracedCyy_v_Re_%d', d0, path, Stot));
    print('-deps',sprintf('%s/%s/Fnl_tracedCyy_v_Re_%d', d0, path, Stot));
    
end
%% Plot joint scatter plots
for k = 1:2
clf, close all,figure
colormap(mk_colmap)
YYt = sol.YYt;
TicksX = [];
TicksY = [];
Stot = k;
for j = 1:Stot
    i = j
    %h = subplot(Stot, Stot, Stot * (Stot-1) + i - Stot * (j-1));
    h = subplot(Stot, Stot, i + Stot * (j-1));
    [n x w] = ksdensity(YYt(:,i));
    [n x] = ksdensity(YYt(:,i),'width',w/3);
    %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
    %     xnorm = (max(x)-min(x))/2;
    %     bar((x-min(x))/xnorm - max((x-min(x))/xnorm)/2,n/max(n),'g');
    %xlim([-1,1]);
    plot(x, n,'-','linewidth',1, 'color',[0.25, 0, 0.25]);
    xlim([- std(x) * 2, std(x) * 2]);
    ylim([0, max(n(:))]);
    TicksX = [TicksX; [ceil(- std(x) * 2*100)/100,0, floor((std(x) * 2)*100)/100]];
    set(h, 'YTick', [0, floor(max(n(:))*100)/100],...
        'XTick', TicksX(i,:))            
    %hold off
    axis square
    axis off
end

for j = 1:Stot
    j
    for i = 1:Stot
        color = 1;
        %h = subplot(Stot, Stot, Stot * (Stot-1) + i - Stot * (j-1));
        h = subplot(Stot, Stot, i + Stot * (j-1));
        if i ~= j
            %plot3k(YYt(:,[i, j, j+1]),[],'.',3);
            plot3k(YYt(:,[i, j, color]),[],'.',3);
            axis([min(YYt(:,i)) max(YYt(:,i)) min(YYt(:,j)) max(YYt(:,j)) min(YYt(:,color)) max(YYt(:,color)) ])
            xlim(TicksX(i,[1,3]));
            ylim(TicksX(j,[1,3]));
            %set(h, 'YTick', [ceil(min(YYt(:, j))*100)/100, floor(max(YYt(:, j))*100)/100],...
            %    'XTick', TicksX(i,:));  
            set(h, 'YTick',TicksX(j,:) ,'XTick', TicksX(i,:));
            %grid off
            view(2)
            axis square
            axis off
        end
        
    end
end
print('-dpng','-r1000',sprintf('%s/%s/Fnl_yyts_%d', d0, path, Stot));        
end
%% Plot joint scatter plots colored by density
for k = 8:8
clf, close all,figure
colormap(mk_colmap)
YYt = sol.YYt;
TicksX = [];
TicksY = [];
Stot = k;
for j = 1:Stot
    i = j
    %h = subplot(Stot, Stot, Stot * (Stot-1) + i - Stot * (j-1));
    h = subplot(Stot, Stot, i + Stot * (j-1));
    [n x w] = ksdensity(YYt(:,i));
    [n x] = ksdensity(YYt(:,i),'width',w/3);
    %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
    %     xnorm = (max(x)-min(x))/2;
    %     bar((x-min(x))/xnorm - max((x-min(x))/xnorm)/2,n/max(n),'g');
    %xlim([-1,1]);
    plot(x, n,'-','linewidth',1, 'color',[0.25, 0, 0.25]);
    xlim([- std(x) * 2, std(x) * 2]);
    ylim([0, max(n(:))]);
    TicksX = [TicksX; [ceil(- std(x) * 2*100)/100,0, floor((std(x) * 2)*100)/100]];
    set(h, 'YTick', [0, floor(max(n(:))*100)/100],...
        'XTick', TicksX(i,:))            
    %hold off
    axis square
    axis off
end

for j = 1:Stot
    j
    for i = 1:Stot
        color = 1;
        %h = subplot(Stot, Stot, Stot * (Stot-1) + i - Stot * (j-1));
        h = subplot(Stot, Stot, i + Stot * (j-1));
        if i ~= j
            %plot3k(YYt(:,[i, j, j+1]),[],'.',3);
            x = YYt(:, i);
            y = YYt(:, j);
            xgrid = linspace(min(x), max(x), app.MC/100);
            ygrid = linspace(min(y), max(y), app.MC/100);
            color = hist2(x, y, xgrid, ygrid);
            color = log(interp2(xgrid', ygrid, color, x, y));
            
            plot3k([YYt(:,[i, j]), color],[],'.',3);
            
            axis([min(YYt(:,i)) max(YYt(:,i)) min(YYt(:,j)) max(YYt(:,j)) log(0) sqrt(max(color)) ])
            xlim(TicksX(i,[1,3]));
            ylim(TicksX(j,[1,3]));
            %set(h, 'YTick', [ceil(min(YYt(:, j))*100)/100, floor(max(YYt(:, j))*100)/100],...
            %    'XTick', TicksX(i,:));  
            set(h, 'YTick',TicksX(j,:) ,'XTick', TicksX(i,:));
            %grid off
            view(2)
            axis square
            axis off
        end
        
    end
end
print('-dpng','-r1000',sprintf('%s/%s/Fnl_yyts_densecol_%d', d0, path, Stot));        
end

%% Plot a bunch of realizations
Nbcs = length(u) - length(uid);
i = 1;
step = app.MC/100 * 0.8;
while i < length(YYt)
    %Make the realization
    u = sol.u + sol.ui * YYtsort(i, :)';
    v = sol.v + sol.vi * YYtsort(i, :)';
    P = sol.P + sol.Pi * YYtsort(i, :)';
    RE_loc = (YYtsort(i, 1) * sol.ui(2,1) + sol.u(2)) / app.nu;
    %Set Standard plotting stuff
    clf
    STREAMCOL = [0.25, 0, 0.25];
    %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 = 2;
    
    %Plot the field
    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
    tmp = u(1:Nbcs);
    u(1:Nbcs) = NaN;
    dudy = (u(Nodeu(2:end-2,2:end-1))-u(Nodeu(3:end-1,2:end-1)))/dy;
    u(1:Nbcs) = tmp;
    tmp = v(1:Nbcs);
    v(1:Nbcs) = NaN;
    dvdx = (v(Nodev(2:end-1,3:end-1))-v(Nodev(2:end-1,2:end-2)))/dx;
    v(1:Nbcs) = tmp;
    vort = dvdx-dudy;
     pcolor(XU(2:end,:),YV(:,2:end),vort),
     caxis([min(-std(vort(~isnan(vort))),-eps),max(std(vort(~isnan(vort))),0)]*STDVORT+mean(vort(~isnan(vort))));
     shading interp
    %tmp = P(1:Nbcs);
    %P(1:Nbcs) = NaN;
    %pcolor(XP,YP,P(NodeP(2:app.Ny +1, 2:app.Nx +1)));
    %shading interp
    %P(1:Nbcs) = tmp;
    %caxis([min(-std(P),-eps),max(std(P),0)]*STDVORT+mean(P));

    h=streamslice(XP,YP,tmp1,tmp2);
    set(h,'color',STREAMCOL);

    % LEVU = linspace(0, max(P(pid)),10);
    % LEVD = linspace(min(P(pid)),0,10);
    % contour(XP,YP,P(NodeP(2:app.Ny+1,2:app.Nx+1)),LEVU,'k');
    % contour(XP,YP,P(NodeP(2:app.Ny+1,2:app.Nx+1)),LEVD,'--k');

    hold off
    axis equal
    xlim(XLIM);ylim(YLIM)
    axis off
    TT = title(sprintf('Time =%03.03g \nRe = %3.1f',sol.T, round(RE_loc*10)/10));
    
    print('-dpng','-r200',sprintf('%s/%s/RE/%04d', d0, path, round(RE_loc*10)));
    if step + i > length(YYt)
        i = length(YYt)-1;
    else
        i = i + step;
    end
end

%% Now do it in groups of RE
close all
RE_locs = (YYtsort(:, 1) * sol.ui(2,1) + sol.u(2)) / app.nu;
i = 1;
RE0 = RE_locs(1);
dRE = 0.1;
step = 5;
while i < length(YYt)
    ids = find((RE_locs >= RE0) & (RE_locs < RE0 + dRE));
    k = 1;
    while (length(ids) < step * 2) && (k < 20)
        ids = find((RE_locs >= RE0) & (RE_locs < RE0 + dRE * 5 * k));
        k = k+1;
    end
    i = max(ids);
    YYtmean = mean(YYtsort(ids, :), 1);
    %Make the realization
    u = sol.u + sol.ui * YYtmean';
    v = sol.v + sol.vi * YYtmean';
    P = sol.P + sol.Pi * YYtmean';
    RE_loc = (YYtmean(1) * sol.ui(2,1) + sol.u(2)) / app.nu;
    %Set Standard plotting stuff
    clf
    STREAMCOL = [0.25, 0, 0.25];
    %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 = 2;
    
    %Plot the field
    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
    tmp = u(1:Nbcs);
    u(1:Nbcs) = NaN;
    dudy = (u(Nodeu(2:end-2,2:end-1))-u(Nodeu(3:end-1,2:end-1)))/dy;
    u(1:Nbcs) = tmp;
    tmp = v(1:Nbcs);
    v(1:Nbcs) = NaN;
    dvdx = (v(Nodev(2:end-1,3:end-1))-v(Nodev(2:end-1,2:end-2)))/dx;
    v(1:Nbcs) = tmp;
    vort = dvdx-dudy;
     pcolor(XU(2:end,:),YV(:,2:end),vort),
     caxis([min(-std(vort(~isnan(vort))),-eps),max(std(vort(~isnan(vort))),0)]*STDVORT+mean(vort(~isnan(vort))));
     shading interp
    %tmp = P(1:Nbcs);
    %P(1:Nbcs) = NaN;
    %pcolor(XP,YP,P(NodeP(2:app.Ny +1, 2:app.Nx +1)));
    %shading interp
    %P(1:Nbcs) = tmp;
    %caxis([min(-std(P),-eps),max(std(P),0)]*STDVORT+mean(P));

    h=streamslice(XP,YP,tmp1,tmp2);
    set(h,'color',STREAMCOL);

    % LEVU = linspace(0, max(P(pid)),10);
    % LEVD = linspace(min(P(pid)),0,10);
    % contour(XP,YP,P(NodeP(2:app.Ny+1,2:app.Nx+1)),LEVU,'k');
    % contour(XP,YP,P(NodeP(2:app.Ny+1,2:app.Nx+1)),LEVD,'--k');

    hold off
    axis equal
    xlim(XLIM);ylim(YLIM)
    axis off
    TT = title(sprintf('Time =%03.03g \nRe = %3.1f',sol.T, round(RE_loc*10)/10));
    
    print('-dpng','-r200',sprintf('%s/%s/RE_M/%04d', d0, path, round(RE_loc*10)));
    if i < length(YYt)
        RE0 = RE_locs(i+1);
    end
    if step + i > length(YYt)
        i = length(YYt)+1;
    else
        i = i + step;
    end
end

