% caseName = 'growthDensitySampleUniform'; % 0.16
% caseName = 'growthDensitySampleNormal'; % 7%
% caseName = 'growthDensitySampleNormalHighRes'; % 2.1%
% caseName = 'growthDensitySampleMetro'; % vol < 4.5%; width < 1.4% (both fgo down)
% caseName = 'growthDensitySampleMetroWidth'; % vol < 2.4; width < 1.4% (mw goes up; vol equilibrates at 1.3%)
% caseName = 'growthDensitySampleMetroVolume'; % vol < 4.5%; width < 4.5% (both go up)
% caseName = 'growthDensitySampleMetroMix'; % vol < 1.9% ; width < 1.9%
caseName = 'growthSampleMetroMixHighRes';

% 1) to obtain the accuracy of mix by normal 1360 sample points would be
%    required compared to 100
% 1.1) Assume you require an accuracy of 0.1%. With normal, you require
%    490000 sample points with metroMix 36100
% 1.2) Over 10 times faster with almost no extra costs (metro and
%    integration is neglible)

%% execute solver
executeCase('init');
tic
executeCase(caseName);
disp('--------HERE-----------')
toc
disp('^^^^^')

%% load result
[result pbeDefinition] = loadCase(caseName)

%% 2D distribution
figure, hold on

hc1max = 0;hc2max = 0;
for iTime = 1:5:(length(result))
    
    % switch through some marker colors
    switch mod((iTime-1)/5, 4)
        case 0, markerString = 'kx';
        case 1, markerString = 'bx';
        case 2, markerString = 'rx';
        case 3, markerString = 'mx';
    end

    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

% add validity cone
% hrep.A = [pbeDefinition.dec.validityCone.A; ...
%           1*eye(2)];
% hrep.h = [pbeDefinition.dec.validityCone.h;...
%           hc1max; hc2max];
% hrep.A = [hrep.A zeros(length(hrep.h), 1); ...
%           0 0 1; 0 0 -1];
% hrep.h = [hrep.h;0; 1];
% viewPolytope(hrep);
% set(gca, 'PlotBoxAspectRatioMode', 'auto', ...
%          'DataAspectRatioMode', 'auto')
         
% axis scaling
xlim([0 hc1max]*1.05)
ylim([0 hc2max]*1.05)
% zlim([0 9e21])
% zlim([0 2e25])
zlim([0 2e23])
% set(gca,'zScale', 'log')
view([13, 16]) 
title('particle density distribution')
zlabel('particle density in #/m^2')
xlabel('hC1 in m');ylabel('hC2 in m')


% %% plot volume over time
% [X, Y] = extractDataFromPbeResult('integralProperty', 'volume', result, pbeDefinition);
% 
% figure,hold on
% plot(X, Y(1,:), 'k-'), plot(X, Y(1,:)-Y(2,:), 'r.'), plot(X, Y(1,:)+Y(2,:), 'r.')
% title('total volume');
% legend({'total volume (left axis)', 'lower error (std)', 'upper error (std)'},...
%     'Location', 'NorthWest');
% xlabel('time in s')
% ylabel('total volume in m^3')
% % error extra y-axis
% axisHandle = axes;
% plot(X, Y(2,:)./Y(1,:), 'k--', 'parent', axisHandle)
% set(axisHandle, 'color', 'none', 'yaxisLocation', 'right', ...
%     'xTick', [])
% legend('coef. of variation (right axis)', 'Location', 'NorthEast')
% ylabel('coefficient of variation')
% 
% %% 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');
% xlabel('time in s')
% ylabel('total volume in m^3')
% % error extra y-axis
% axisHandle = axes;
% plot(X, Y(2,:)./Y(1,:), 'k--', 'parent', axisHandle)
% set(axisHandle, 'color', 'none', 'yaxisLocation', 'right', ...
%     'xTick', [])
% ylabel('coefficient of variation')
% 
% %% plot surface over time
% figure,hold on
% [X, Y] = extractDataFromPbeResult('integralProperty', 'surface', result, pbeDefinition);
% plot(X, Y(1,:), 'k-'), plot(X, Y(1,:)-Y(2,:), 'r.'), plot(X, Y(1,:)+Y(2,:), 'r.')
% title('surface');
% xlabel('time in s')
% ylabel('total surface in m^2')
% % error extra y-axis
% axisHandle = axes;
% plot(X, Y(2,:)./Y(1,:), 'k--', 'parent', axisHandle)
% set(axisHandle, 'color', 'none', 'yaxisLocation', 'right', ...
%     'xTick', [])
% ylabel('coefficient of variation')
% 
% %% 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(X1, 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 = 0; binUpper = 2e-13; binSpacing = 3e-15;
% % binLower = 1e-5; binUpper = 9e-5; binSpacing = 3e-6;
% boundDefinitionMatrix = [binLower:binSpacing:(binUpper - binSpacing);...
%         (binLower + binSpacing):binSpacing:binUpper];
% [T, Y] = extractDataFromPbeResult('propertyDistribution', ...
%     {'volume', 'volume', boundDefinitionMatrix}, result, pbeDefinition);
% %%
% pivotPositionVector = mean(boundDefinitionMatrix, 1);
% 
% nPivot = length(pivotPositionVector);
% figure,hold on
% for iTime = 1:5:(length(result))
%     
% %     plot(result(iTime).hcMatrix(1,1), ...
% %          result(iTime).hcMatrix(2,1), ...
% %           'm.', 'MarkerSize', 20)
%     switch mod((iTime-1)/5, 4)
%         case 0, markerString = 'kx-';
%         case 1, markerString = 'bx-';
%         case 2, markerString = 'rx-';
%         case 3, markerString = 'mx-';
%     end
%     plot(pivotPositionVector, Y(1:nPivot, iTime), markerString, ...
%         'lineWidth', 1.5)
%     % 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
% title('particle volume distr. density over volume')
% ylabel('particle volume distr. density in m^3/m^3')
% xlabel('particle volume in m^3')
% % set(gca,'xscale', 'log')
