% Analyze Wake Data
% Hannah Ross
% 8/2/2019

% Analyzes velocity data collected in the wake of a cross-flow turbine at
% two different blockage ratios.

clear

%% Load data
load('cleanData.mat')
fieldNames = fields(pos);

% Interpolate data to higher resolution
spacing = 0.1;
for i = 1:length(fieldNames)
    xLocs = unique(pos.(fieldNames{i})(:,1));
    for j = 1:length(xLocs)
        xInds = pos.(fieldNames{i}) == xLocs(j);
        posCells.(fieldNames{i}){j,1} = pos.(fieldNames{i})(xInds,:);
        velCells.(fieldNames{i}){j,1} = vel.(fieldNames{i})(xInds,:);
        dissCells.(fieldNames{i}){j,1} = diss.(fieldNames{i})(xInds,:);
        TKECells.(fieldNames{i}){j,1} = TKE.(fieldNames{i})(xInds,:);
        yMin = min(posCells.(fieldNames{i}){j,1}(:,2));
        yMax = max(posCells.(fieldNames{i}){j,1}(:,2));
        zMin = min(posCells.(fieldNames{i}){j,1}(:,3));
        zMax = max(posCells.(fieldNames{i}){j,1}(:,3));
        [posGrid.(fieldNames{i}){j,1}(:,:,2),posGrid.(fieldNames{i}){j,1}(:,:,3)] = meshgrid(yMin:spacing:yMax,zMin:spacing:zMax);
        posGrid.(fieldNames{i}){j,1}(:,:,1) = posGrid.(fieldNames{i}){j,1}(:,:,2)*0+posCells.(fieldNames{i}){j,1}(1,1);
        for k = 1:size(velCells.(fieldNames{i}){j,1},2)
            velGrid.(fieldNames{i}){j,1}(:,:,k) = griddata(posCells.(fieldNames{i}){j,1}(:,2),posCells.(fieldNames{i}){j,1}(:,3),velCells.(fieldNames{i}){j,1}(:,k),posGrid.(fieldNames{i}){j,1}(:,:,2),posGrid.(fieldNames{i}){j,1}(:,:,3),'v4');
        end
        dissGrid.(fieldNames{i}){j,1} = griddata(posCells.(fieldNames{i}){j,1}(:,2),posCells.(fieldNames{i}){j,1}(:,3),dissCells.(fieldNames{i}){j,1},posGrid.(fieldNames{i}){j,1}(:,:,2),posGrid.(fieldNames{i}){j,1}(:,:,3),'v4');
        dissGrid.(fieldNames{i}){j,1}(dissGrid.(fieldNames{i}){j,1} < 0) = 0;
        TKEGrid.(fieldNames{i}){j,1} = griddata(posCells.(fieldNames{i}){j,1}(:,2),posCells.(fieldNames{i}){j,1}(:,3),TKECells.(fieldNames{i}){j,1},posGrid.(fieldNames{i}){j,1}(:,:,2),posGrid.(fieldNames{i}){j,1}(:,:,3),'v4');
        clear xInds yMin yMax zMin zMax    
    end
    clear xLocs
end

%% Calculate wake and bypass flow positions
[ymask1,ymask2,zmask1,zmask2] = maskRegions(posGrid);

% Select TKE peaks along every row and column
for i = 1:length(fieldNames)
    for j = 1:length(TKEGrid.(fieldNames{i}))
        for k = 1:size(TKEGrid.(fieldNames{i}){j,1},1)
            [~,temp] = findpeaks(TKEGrid.(fieldNames{i}){j,1}(k,:).*ymask1.(fieldNames{i}){j,1}(k,:),'SortStr','descend');
            yLocTube.(fieldNames{i}){j,1}(k,1) = posGrid.(fieldNames{i}){j,1}(k,temp(1),2);
            clear temp
            [~,temp] = findpeaks(TKEGrid.(fieldNames{i}){j,1}(k,:).*ymask2.(fieldNames{i}){j,1}(k,:),'SortStr','descend');
            yLocTube.(fieldNames{i}){j,1}(k,2) = posGrid.(fieldNames{i}){j,1}(k,temp(1),2);
            clear temp
        end
        for m = 1:size(TKEGrid.(fieldNames{i}){j,1},2)
            [~,temp] = findpeaks(TKEGrid.(fieldNames{i}){j,1}(:,m).*zmask1.(fieldNames{i}){j,1}(:,m),'SortStr','descend');
            zLocTube.(fieldNames{i}){j,1}(1,m) = posGrid.(fieldNames{i}){j,1}(temp(1),m,3);
            clear temp
            [~,temp] = findpeaks(TKEGrid.(fieldNames{i}){j,1}(:,m).*zmask2.(fieldNames{i}){j,1}(:,m),'SortStr','descend');
            zLocTube.(fieldNames{i}){j,1}(2,m) = posGrid.(fieldNames{i}){j,1}(temp(1),m,3);
            clear temp
        end
    end
end

% Test each point to see if it is within the Y and Z bounds
for i = 1:length(fieldNames)
    for j = 1:length(posGrid.(fieldNames{i}))
        for k = 1:size(posGrid.(fieldNames{i}){j,1},1)
            for m = 1:size(posGrid.(fieldNames{i}){j,1},2)
                if posGrid.(fieldNames{i}){j,1}(k,m,2) > yLocTube.(fieldNames{i}){j,1}(k,1) && ... 
                    posGrid.(fieldNames{i}){j,1}(k,m,2) < yLocTube.(fieldNames{i}){j,1}(k,2) && ... 
                    posGrid.(fieldNames{i}){j,1}(k,m,3) > zLocTube.(fieldNames{i}){j,1}(1,m) && ... 
                    posGrid.(fieldNames{i}){j,1}(k,m,3) < zLocTube.(fieldNames{i}){j,1}(2,m)
                    indsCoreNaN.(fieldNames{i}){j,1}(k,m) = 1;
                else
                    indsCoreNaN.(fieldNames{i}){j,1}(k,m) = NaN;
                end
            end
        end
    end
end

% Delete errant regions
indsCoreNaN.MongoBMSC{2,1}(73:78,750) = NaN;

% Convert to logical
for i = 1:length(fieldNames)
    for j = 1:length(indsCoreNaN.(fieldNames{i}))
        for k = 1:size(indsCoreNaN.(fieldNames{i}){j,1},1)
            for m = 1:size(indsCoreNaN.(fieldNames{i}){j,1},2)
                indsCore.(fieldNames{i}){j,1}(k,m) = nansum([indsCoreNaN.(fieldNames{i}){j,1}(k,m) 0]);
            end
        end
        indsBypass.(fieldNames{i}){j,1} = ~indsCore.(fieldNames{i}){j,1};
    end
end

% Convert bypass inds to NaN
for i = 1:length(fieldNames)
    for j = 1:length(indsBypass.(fieldNames{i}))
        for k = 1:size(indsBypass.(fieldNames{i}){j,1},1)
            for m = 1:size(indsBypass.(fieldNames{i}){j,1},2)
                if indsBypass.(fieldNames{i}){j,1}(k,m) == 0
                    indsBypassNaN.(fieldNames{i}){j,1}(k,m) = NaN;
                elseif indsBypass.(fieldNames{i}){j,1}(k,m) == 1
                    indsBypassNaN.(fieldNames{i}){j,1}(k,m) = double(indsBypass.(fieldNames{i}){j,1}(k,m));
                end
            end
        end
    end
end

% Detect boundaries of continuous regions
for i = 1:length(fieldNames)
    for j = 1:length(indsCore.(fieldNames{i}))
        bounds.(fieldNames{i})(j,1) = bwboundaries(indsCore.(fieldNames{i}){j,1});
    end
end

% Convert from pixels to cm
for i = 1:length(fieldNames)
    for j = 1:length(posGrid.(fieldNames{i}))
        xBounds.(fieldNames{i}){j,1} = [min(min(posGrid.(fieldNames{i}){j,1}(:,:,2))) max(max(posGrid.(fieldNames{i}){j,1}(:,:,2)))];
        xConv.(fieldNames{i}){j,1} = (xBounds.(fieldNames{i}){j,1}(2) - xBounds.(fieldNames{i}){j,1}(1))/size(posGrid.(fieldNames{i}){j,1},2);
        yBounds.(fieldNames{i}){j,1} = [min(min(posGrid.(fieldNames{i}){j,1}(:,:,3))) max(max(posGrid.(fieldNames{i}){j,1}(:,:,3)))];
        yConv.(fieldNames{i}){j,1} = (yBounds.(fieldNames{i}){j,1}(2) - yBounds.(fieldNames{i}){j,1}(1))/size(posGrid.(fieldNames{i}){j,1},1);
        aConv.(fieldNames{i}){j,1} = ((xBounds.(fieldNames{i}){j,1}(2) - xBounds.(fieldNames{i}){j,1}(1))*(yBounds.(fieldNames{i}){j,1}(2) - yBounds.(fieldNames{i}){j,1}(1)))/(size(posGrid.(fieldNames{i}){j,1},2)*size(posGrid.(fieldNames{i}){j,1},1));
        bounds.(fieldNames{i}){j,1}(:,1) = bounds.(fieldNames{i}){j,1}(:,1)*yConv.(fieldNames{i}){j,1} + yBounds.(fieldNames{i}){j,1}(1);
        bounds.(fieldNames{i}){j,1}(:,2) = bounds.(fieldNames{i}){j,1}(:,2)*xConv.(fieldNames{i}){j,1} + xBounds.(fieldNames{i}){j,1}(1);
    end
end        

%% Calculate mean core and bypass flow rates
for i = 1:length(fieldNames)
    for j = 1:length(velGrid.(fieldNames{i}))
        velMean.(fieldNames{i})(j,1) = mean(mean(velGrid.(fieldNames{i}){j,1}(:,:,1)));
        coreVelMean.(fieldNames{i})(j,1) = nanmean(nanmean(velGrid.(fieldNames{i}){j,1}(:,:,1).*indsCoreNaN.(fieldNames{i}){j,1}));
        bypassVelMean.(fieldNames{i})(j,1) = nanmean(nanmean(velGrid.(fieldNames{i}){j,1}(:,:,1).*indsBypassNaN.(fieldNames{i}){j,1}));
    end
end

%% Calculate mean TKE
for i = 1:length(fieldNames)
    for j = 1:length(TKEGrid.(fieldNames{i}))
        TKEMean.(fieldNames{i})(j,1) = mean(mean(TKEGrid.(fieldNames{i}){j,1}));
    end
end

%% Calculate wake area
for i = 1:length(fieldNames)
    for j = 1:length(indsCore.(fieldNames{i}))
        b = vision.BlobAnalysis;
        [area] = b(logical(indsCore.(fieldNames{i}){j,1}));
        wakeArea.(fieldNames{i})(j,1) = area*aConv.(fieldNames{i}){j,1};
        clear b area
    end
end

%% Calculate mean dissipation rate
for i = 1:length(fieldNames)
    for j = 1:length(dissGrid.(fieldNames{i}))
        dissMean.(fieldNames{i})(j,1) = mean(mean(dissGrid.(fieldNames{i}){j,1}));
    end
end