clear all;
caseName = 'growthDensitySampleUniform'; % 0.16

%% load result
[result pbeDefinition] = loadCase(caseName);
close all;

%% parameters
m = 1.125e-6;
off = 0.7*m;
d = m/100;

%% initial distribution
figure(1),hold on
set(1, 'renderer', 'zbuffer')
gridSpacing = (m-off):d:(m+off);
[X Y] = meshgrid(gridSpacing);
Z = pbeDefinition.seeding.fhandle([X(:)'; Y(:)']);
Z = reshape(Z, size(X));
hSurf = surfc(X, Y, Z, 'edgeColor', 'none')

% add validity cone
hrep.A = [pbeDefinition.dec.validityCone.A; ...
          1*eye(2)];
hrep.h = [pbeDefinition.dec.validityCone.h;...
          m+off; m+off];
hrep.A = [hrep.A zeros(length(hrep.h), 1); ...
          0 0 1; 0 0 -1];
hrep.h = [hrep.h;0; 0];
pHandle = viewPolytope(hrep, 'color', [0.5 1 0.5], 'shade', 0.2);
for iPatch = 1:length(pHandle)
    set(pHandle(iPatch), 'zData', ...
        get(pHandle(iPatch), 'zData') + 1e21) 
end

set(gca, 'PlotBoxAspectRatioMode', 'auto', ...
         'DataAspectRatioMode', 'auto')

xlim([m-off, m+off]);ylim([m-off, m+off])
zlim([0 2e23])
view([10 25])
box off
axis off

% print
print -dpng -r300 figure/uniform/initialDistribution.png

% replace distribution with sampled distribution
delete(hSurf)
[X, Y, Z] = computePatchDataFromPointCloud(...
        [result(1).hcMatrix; ...
         result(1).numberVector .* result(1).probabilityVector]);
patch(X, Y, Z, Z, 'lineWidth', 0.5)

% add sampling area (uniform)
vrep = pbeDefinition.seeding.vrep;
vrep.V = [vrep.V zeros(size(vrep.V, 1), 1); ...
          vrep.V 1e-6*ones(size(vrep.V, 1), 1)];
pHandle = viewPolytope(vrep, 'shade', 0, 'edgeColor', [1 0 0], ...
    'lineWidth', 3);
for iPatch = 1:length(pHandle)
    set(pHandle(iPatch), 'zData', ...
        get(pHandle(iPatch), 'zData') + 2e21) 
end

% add sampling points
plot3(result(1).hcMatrix(1,:), ...
      result(1).hcMatrix(2,:), ...
      result(1).numberVector .* result(1).probabilityVector, ...
      'rd', 'markerSize', 5, 'markerEdgeColor', [1 0.9 0.9], ...
      'markerFaceColor', [1 0 0], 'lineWidth', 1);
  
set(gca, 'PlotBoxAspectRatioMode', 'auto', ...
      'DataAspectRatioMode', 'auto')
view([10 25])

% print
print -dpng -r300 figure/uniform/initialDistributionSampled.png
close all
  
%% View 2D trajectories
figure, hold on
% 
hc1max = 0;hc2max = 0;
for iTime = 1:7:(length(result))
    
%     plot(result(iTime).hcMatrix(1,1), ...
%          result(iTime).hcMatrix(2,1), ...
%           'm.', 'MarkerSize', 20)
    switch mod((iTime-1)/7, 3)
        case 0, markerString = 'k.';
        case 1, markerString = 'b.';
        case 2, markerString = 'r.';
        case 3, markerString = 'm.';
    end

    plot(result(iTime).hcMatrix(1,:), ...
         result(iTime).hcMatrix(2,:), ...
         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];
viewPolytope(hrep, 'color', [0.5 1 0.5], 'shade', 0.2, 'edgeColor', 'none');

% axis scaling
set(gca,'FontSize',20)
xlim([0 hc1max]*1.05)
ylim([0 hc2max]*1.05)
xlabel('h_{C,1} in m');ylabel('h_{C,2} in m')

% print
% print -dpng -r300 figure/uniform/growthOverview.png
% close all

%% sampled distributions
timeVector = [1 150 250 396];
disp('starting')
off = 2e-6;
for iTime = timeVector
    thisTime = result(iTime).time
    thisSimFraction = iTime / length(result)
    
    [X, Y, Z] = computePatchDataFromPointCloud(...
        [result(iTime).hcMatrix; ...
        result(iTime).numberVector .* result(iTime).probabilityVector]);
    mHere = mean(result(iTime).hcMatrix, 2);
    
    hFigure = figure, set(hFigure, 'renderer', 'zbuffer');
    hold on
    patch(X, Y, Z, Z, 'lineWidth', 0.5)
    plot3(result(iTime).hcMatrix(1,:), result(iTime).hcMatrix(2,:), ...
          result(iTime).numberVector .* result(iTime).probabilityVector, ...
          'rd')
    xlim([mHere(1)-off, mHere(1)+off]);
    ylim([mHere(2)-off, mHere(2)+off]);
    zlim([0 3e23]);set(gca, 'clim', [0 3e23]);
    view([0 90])
%     axis off
    box on
    
%     print('-dpng', '-r300', ['figure/uniform/density_' num2str(iTime) '.png']);
%     close all
end

%% 3D distribution
figure, hold on
hc1max = 0;hc2max = 0;

for iTime = 1:15:(length(result))
    result(iTime).time
    switch mod((iTime-1)/15, 4)
        case 0, markerString = 'k.';
        case 1, markerString = 'b.';
        case 2, markerString = 'r.';
        case 3, markerString = 'm.';
    end

    plot3(result(iTime).hcMatrix(1,:), ...
          result(iTime).hcMatrix(2,:), ...
          result(iTime).numberVector .* result(iTime).probabilityVector)
%     result(iTime).probabilityVector
    hc1max = max(hc1max, result(iTime).hcMatrix(1,:));
    hc1max = max(hc1max);
    hc2max = max(hc2max, result(iTime).hcMatrix(2,:));
    hc2max = max(hc2max);
end

% axis scaling
set(gca,'FontSize',20)
xlim([0 hc1max]*1.05)
ylim([0 hc2max]*1.05)
zlim([0 3e23])
% set(gca, 'zScale', 'log')
xlabel('h_{C,1} in m');ylabel('h_{C,2} in m')
view([20 35])

% % print
print -dpng -r300 figure/uniform/growthDensityOverview.png
close all

%% Detail of final distribution
plot3(result(end).hcMatrix(1,:), ...
      result(end).hcMatrix(2,:), ...
      result(end).numberVector .* result(end).probabilityVector, 'b.', ...
      'markerSize', 15)
% axis scaling
set(gca,'FontSize',20)
xlabel('h_{C,1} in m');ylabel('h_{C,2} in m');
zlabel({'particle number'; 'density in #/m^2'})
view([45 25])

% print
print -dpng -r300 figure/uniform/growthDensityFinal.png
close all



%% plot volume over time
figure,hold on
set(1, 'paperPosition', [0 0 30 15])

[X, Y] = extractDataFromPbeResult('integralProperty', 'volume', result, pbeDefinition);

plot(X, Y(1,:), 'k-', 'lineWidth', 3)
plot(X, Y(1,:)-Y(2,:), 'r.'), plot(X, Y(1,:)+Y(2,:), 'r.')
set(gca,'FontSize',20)
box off
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, 'lineWidth', 3)
ylim([0.15 0.18])
set(gca,'FontSize',20)
box off
set(axisHandle, 'color', 'none', 'yaxisLocation', 'right', ...
    'xTick', [])
legend('coef. of variation (right axis)', 'Location', 'NorthEast')
ylabel('coefficient of variation')

% print
print -dpng -r300 figure/uniform/volumeOverTime.png
close all

%% volume distribution
figure(1),hold on
set(1, 'paperPosition', [0 0 30 15])

% volume distribution density, represented by volume bins
binLower = 0; binUpper = 2e-13; binSpacing = 3e-15;
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);
for iTime = 1:10:(length(result))
    
    markerString = 'k-';
    plot(pivotPositionVector, Y(1:nPivot, iTime), markerString, ...
        'lineWidth', 2)
    % error
    plot(pivotPositionVector, ...
         Y(1:nPivot, iTime) + Y(nPivot+(1:nPivot), iTime), ...
         'r.', 'markerSize', 12)
    plot(pivotPositionVector, ...
         Y(1:nPivot, iTime) - Y(nPivot+(1:nPivot), iTime), ...
         'r.', 'markerSize', 12)
end
set(gca,'FontSize',20)
ylabel('particle volume distr. density in m^3/m^3')
xlabel('particle volume in m^3')

% % print
print -dpng -r300 figure/uniform/volumeDensityDistr.png
close all
