function usePrivate = analyseScript
%% initialize
clc,close all
load localData/exp255_500rpm_f_resultStruct.mat

%% comments
% air bubbles: can be well identified by relative additive roundness (see
%   figure 20) without filtering

%% evaluate data
% figure(1),hold on,title('lambda')
figure(2),hold on,title('width')
figure(3),hold on,title('roughness')

figure(4),hold on,title('relative additive sphericity')
figure(5),hold on,title('roughness by width')

figure(20),hold on,title('phase: lambda vs. mean width')
% figure(21),hold on,title('phase: roughness vs. mean width')
figure(22),hold on,title('phase: lambdaH vs. lambda M')
figure(23),hold on,title('phase: roundness vs. brightFactorX')

figure(30),hold on,title('distribution: mean width');
figure(31),hold on,title('distribution: roundness');

nSampleVector = nan(1, length(timeResult));
timeVector = nan(1, length(timeResult));

for iTime = 1:length(timeResult)
    % data vectors
    thisTime = timeGrid(iTime);
    timeVector(iTime) = thisTime;
    thisColor = iTime/length(timeResult) * [1 0 0];
    thisLambdaVector = timeResult(iTime).lambda;
    thisWidthVector = timeResult(iTime).meanWidth;
    thisRoughnessVector = timeResult(iTime).roughness;
    thisConvexityVector = timeResult(iTime).convexity;
    thisSphericityVector = timeResult(iTime).sphericity;
%     thisSphericityVector = min([timeResult(iTime).sphericityH; timeResult(iTime).sphericityM]);
    thisSphericity2Vector = timeResult(iTime).sphericity2;
    thisLambdaHVector = timeResult(iTime).lambdaH;
    thisLambdaMVector = timeResult(iTime).lambdaM;
    thisRoundnessVector = (2*thisLambdaVector ./ thisWidthVector);
    thisBrightHVector = timeResult(iTime).brightFactorH;
    thisBrightMVector = timeResult(iTime).brightFactorM;
                      
    % check for relative roundness < 1 (sanity)
    if any(thisRoundnessVector > 1.01)
        warning('input data is not sane: 2[lambda] > mean width');
    end
          
    filterVector = true(size(thisLambdaVector));
    % filter for roundness (bubbles)
%     filterVector = filterVector & abs(thisRoundnessVector - 1) > 0.01;
    % filter roughness (goodness of fit)
%     filterVector = filterVector & (thisRoughnessVector < 15);
%     filterVector = filterVector & (thisConvexityVector > 0.95);
%     filterVector = filterVector & (thisSphericityVector > 0.1);
%     filterVector = filterVector & (thisSphericity2Vector < 0.96);
    filterVector = filterVector & (thisRoughnessVector./thisWidthVector < 0.1);
    % filter for mean width (jitter in very-large particles)
%     filterVector = filterVector & (thisWidthVector < ...
%         mean(thisWidthVector) + 3*std(thisWidthVector));
    filterVector = filterVector & (thisWidthVector < 2000);
    % filter for mean width (minimum size)
    filterVector = filterVector & (thisWidthVector > 100);
    filterVector = filterVector & (thisBrightHVector > 0.015);
    filterVector = filterVector & (thisBrightMVector > 0.015);
    
    % apply filter
    thisLambdaVector = thisLambdaVector(filterVector);
    thisWidthVector = thisWidthVector(filterVector);
    thisRoughnessVector = thisRoughnessVector(filterVector);
    thisRoundnessVector = thisRoundnessVector(filterVector);
       
    thisConvexityVector =thisConvexityVector(filterVector);
    thisSphericityVector = thisSphericityVector(filterVector);
    thisSphericity2Vector = thisSphericity2Vector(filterVector);
    thisLambdaHVector = thisLambdaHVector(filterVector);
    thisLambdaMVector = thisLambdaMVector(filterVector);
    thisBrightHVector = thisBrightHVector(filterVector);
    thisBrightMVector = thisBrightMVector(filterVector);
    
    nSample = length(thisWidthVector);
    nSampleVector(iTime) = nSample;
%     figure(1)
%     plotTimeVector(thisTime, thisLambdaVector)
%     
    figure(2)
    plotTimeVector(thisTime, thisWidthVector)
%     
    figure(3)
    plotTimeVector(thisTime, thisRoughnessVector)

    figure(4)
    plotTimeVector(thisTime, 2*thisLambdaVector ./ thisWidthVector)
    
    figure(5)
    plotTimeVector(thisTime, thisRoughnessVector./thisWidthVector)

    figure(20)
    plot(thisWidthVector, thisLambdaVector, 'r.', 'color', thisColor)

%     figure(21)
%     plot(thisWidthVector, thisRoughnessVector, 'r.', 'color', thisColor)   
    
    figure(22)
    plot(thisLambdaHVector, thisLambdaMVector, 'r.', 'color', thisColor)
    
    figure(23)
    plot(thisRoundnessVector, thisBrightHVector, 'r.', 'color', thisColor)
    plot(thisRoundnessVector, thisBrightMVector, 'r.', 'color', thisColor)
    
    figure(30)
    thisWidthVector = sort(thisWidthVector);
    plotDistribution(thisWidthVector, 20, thisColor, ...
        [0 100 150 200 250 300 400 500 600 700])
%     plot(thisWidthVector, (1:nSample)./nSample, 'r.-', 'color', thisColor)
    
    figure(31)
    thisRoundnessVector = sort(thisRoundnessVector);
    plotDistribution(thisRoundnessVector, 20, thisColor, ...
        [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1]);
%     plot(thisRoundnessVector, (1:nSample)./nSample, 'r.-', 'color', thisColor)


%     figure(40+iTime)
%     plotDistribution2D(thisWidthVector, thisRoundnessVector);
    
end

nSampleVector
timeVector