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

%Set Case number
CASE = 2;

%Find the directory for this case number
path = dir([d0, '*Case', num2str(CASE)]);

%Error checking
if isempty(path)
    fprintf('This case number %d does not exist.\n', CASE)
elseif length(path) >1
    fprintf('Multiple directories with the same case number %d exist.\n', CASE)
else
    d0 = [d0, path(1).name, '/'];
end
    
%Load the params file
load([d0, 'param.mat']);

%Find all the casefiles for this time
casefiles = dir([d0, '0*.mat']);

for i=1:length(casefiles)
    %%
    DOsol = load([d0, casefiles(i).name]);
    Nbcs = length(DOsol.u) - length(uid);

    %In this loop we make some of the DO realizations, let's make 100
    Nrealiz = 200;
    
    %Let's record the Reynold's number for the realization we calculate
    RE = zeros(Nrealiz, 1);
    
    %Let's also record the recirculation zone length
    l_r = zeros(Nrealiz, 1);
    
    realid = floor(linspace(1, app.MC, Nrealiz));
    for j=1:Nrealiz
        %Let's create the j'th realization
        u = DOsol.u + DOsol.ui * DOsol.YYt(realid(j), :)';
        v = DOsol.v + DOsol.vi * DOsol.YYt(realid(j), :)';
        P = DOsol.P + DOsol.Pi * DOsol.YYt(realid(j), :)';
        %And calculate the Reynold's number for this realization
        RE(j) = u(2,1) / app.nu;
        
        %Transform to X, Y coordinate space
        u(1:Nbcs) = NaN; %This is for plotting
        u = u(Nodeu(2:app.Ny+1, 2:app.Nx));
        v(1:Nbcs) = NaN; %This is for plotting
        v = v(Nodev(2:app.Ny, 2:app.Nx+1));
        P(1:Nbcs) = NaN; %This is for plotting
        P = P(NodeP(2:app.Ny+1,2:app.Nx+1));

        %Lets plot the u-velocity
        %pcolor(XU, YU, u), colorbar, shading interp
        %pcolor(XV, YV, v), colorbar, shading interp
        %pcolor(XP, YP, P), colorbar, shading interp

        %Find a particular x, y value range in u (We could have done this
        %outside the loop -- it would be faster
        %This find the center-line u-velocity
        ids = (XU > 5) & abs(YU - 1.5) < (dy * 0.9);

        %Next, let's find out where u changes sign, i.e, the length of the 
        %recirculation zone
        u_c = u(ids);
        x_c = XU(ids);

        id_r = find(sign(u_c(2:end)) + sign(u_c(1:end-1)) == 0, 1) + 1;
        if id_r == 1
            id_r = 2;
        end
        %Now we know the zero velocity lies somewhere between 
        %u_centerline(id_r) and u_centerline(ids_r-1), so, simple linear
        %interpolation:
        l_r(j) = ( u_c(id_r-1) * x_c(id_r) - u_c(id_r) * x_c(id_r-1) )  ...
                / (u_c(id_r-1) - u_c(id_r) );
    end
    %Finally, let's plot the recirculation zone length over RE for this
    %time step
    plot(RE, l_r, '.')
    ylim([6, 8])
    xlim([37, 48])
    title(sprintf('Recirculation zone length vs. Reynolds number\n at time %0.2g', DOsol.T))
    drawnow
    
end


