%% Execute and manual analyse script: growth
close all
clear all

%% solve
% executePbeCase('init')
executePbeCase('seededGrowthUsingDensity');

%% load
% loadPbeResult('init')
[result pbeDefinition] = loadPbeResult('seededGrowthUsingDensity');
result

%% analyse
iTime = 20
% result(iTime).time
% figure,plot3(result(iTime).hcMatrix(1,:), ...
%           result(iTime).hcMatrix(2,:), ...
%           result(iTime).numberVector .* result(iTime).probabilityVector, 'rd')

figure, hold on

hc1max = 0;hc2max = 0;
for iTime = 1:10:(length(result))
    
%     plot(result(iTime).hcMatrix(1,1), ...
%          result(iTime).hcMatrix(2,1), ...
%           'm.', 'MarkerSize', 20)
    switch mod((iTime-1)/10, 3)
        case 0, markerString = 'kx';
        case 1, markerString = 'bx';
        case 2, markerString = 'rx';
        case 3, markerString = 'mx';
    end
    result(iTime)
    plot3(result(iTime).hcMatrix(1,:), ...
          result(iTime).hcMatrix(2,:), ...
          result(iTime).numberVector .* result(iTime).probabilityVector, ...
          markerString)

      
    hc1max = max(hc1max, result(iTime).hcMatrix(1,:));
    hc1max = max(hc1max);
    hc2max = max(hc2max, result(iTime).hcMatrix(2,:));
    hc2max = max(hc2max);
end

xlim([0 hc1max]*1.05)
ylim([0 hc2max]*1.05)
zlim([0 5e21])
% set(gca,'zScale', 'log')
view([13, 16]) 

%% plot volume over time
figure,hold on
[X, Y] = extractDataFromPbeResult('integralProperty', 'volume', result, pbeDefinition);
plot(X, Y(1,:), 'k-'), plot(X, Y(1,:)-Y(2,:), 'r.'), plot(X, Y(1,:)+Y(2,:), 'r.')
title('total volume');
% error extra y-axis
axisHandle = axes;
plot(X, Y(2,:)./Y(1,:), 'k-', 'parent', axisHandle)
set(axisHandle, 'color', 'none', 'yaxisLocation', 'right', ...
    'xTick', [])

%% plot mean width over time
figure,hold on
[X, Y] = extractDataFromPbeResult('integralProperty', 'meanWidth', result, pbeDefinition);
plot(X, Y(1,:), 'k-'), plot(X, Y(1,:)-Y(2,:), 'r.'), plot(X, Y(1,:)+Y(2,:), 'r.')
title('mean width');
% error extra y-axis
axisHandle = axes;
plot(X, Y(2,:)./Y(1,:), 'k-', 'parent', axisHandle)
set(axisHandle, 'color', 'none', 'yaxisLocation', 'right', ...
    'xTick', [])

%% plot shape coordinates
figure,hold on
[X1, Y1] = extractDataFromPbeResult('genericintegralproperty', @(H, up)(H(1,:)), result, pbeDefinition);
plot(X1, Y1(1,:), 'k-'), plot(X1, Y1(1,:)-Y1(2,:), 'r.'), plot(X, Y1(1,:)+Y1(2,:), 'r.')
[X2, Y2] = extractDataFromPbeResult('genericintegralproperty', @(H, up)(H(2,:)), result, pbeDefinition);
plot(X2, Y2(1,:), 'r-'), plot(X2, Y2(1,:)-Y2(2,:), 'k.'), plot(X2, Y2(1,:)+Y2(2,:), 'k.')
% error extra y-axis
axisHandle = axes;hold on
plot(X1, Y1(2,:)./Y1(1,:), 'r-', 'parent', axisHandle)
plot(X2, Y2(2,:)./Y2(1,:), 'k-', 'parent', axisHandle)
set(axisHandle, 'color', 'none', 'yaxisLocation', 'right', ...
    'xTick', [])

title('shape coordinates');

%% volume distribution
% volume distribution density, represented by volume bins
% binLower = 1e-15; binUpper = 1.5e-13; binSpacing = 2e-15;
binLower = 1e-5; binUpper = 8e-5; binSpacing = 4e-6;
boundDefinitionMatrix = [binLower:binSpacing:(binUpper - binSpacing);...
        (binLower + binSpacing):binSpacing:binUpper];
[T, Y] = extractDataFromPbeResult('propertyDistribution', ...
    {'number', 'meanWidth', boundDefinitionMatrix}, result, pbeDefinition);

pivotPositionVector = mean(boundDefinitionMatrix, 1);

nPivot = length(pivotPositionVector);
figure,hold on
for iTime = 1:10:(length(result))
    
%     plot(result(iTime).hcMatrix(1,1), ...
%          result(iTime).hcMatrix(2,1), ...
%           'm.', 'MarkerSize', 20)
    switch mod((iTime-1)/10, 3)
        case 0, markerString = 'kx-';
        case 1, markerString = 'bx-';
        case 2, markerString = 'rx-';
        case 3, markerString = 'mx-';
    end
    plot(pivotPositionVector, Y(1:nPivot, iTime), markerString)
    % error
    plot(pivotPositionVector, ...
         Y(1:nPivot, iTime) + Y(nPivot+(1:nPivot), iTime), ...
         [markerString(1) '.'])
    plot(pivotPositionVector, ...
         Y(1:nPivot, iTime) - Y(nPivot+(1:nPivot), iTime), ...
         [markerString(1) '.'])
end
