%% Development example
% Run solver and display results for an initially selected example
clear all
close all

data = struct([]);
% data(end+1).caseName = 'growthCommon';
data(end+1).caseName = 'growthNoDerounding';
data(end+1).caseName = 'growthSlowDerounding';
data(end+1).caseName = 'growthFastDerounding';
obtainData('init')

%% execute
for iCase = 1:length(data)
    pbeDefinition = obtainData(data(iCase).caseName);
    result = solvePbe(pbeDefinition);
    save(['data' filesep data(iCase).caseName '.mat'], 'result');
end

%% load
for iCase = 1:length(data)
    iCase
    data(iCase).pbeDefinition = obtainData(data(iCase).caseName);
    result = load(['data' filesep data(iCase).caseName '.mat']);
    data(iCase).result = result.result;
end

%% summarize/evaluate
for iCase = 1:length(data)
    % size measures
    [data(iCase).time, data(iCase).mw] = extractDataFromPbeResult(...
        'bulkProperty', 4, ...
        data(iCase).result, data(iCase).pbeDefinition);
    [data(iCase).time, data(iCase).surf] = extractDataFromPbeResult(...
        'bulkProperty', 3, ...
        data(iCase).result, data(iCase).pbeDefinition);
    [data(iCase).time, data(iCase).vol] = extractDataFromPbeResult(...
        'bulkProperty', 2, ...
        data(iCase).result, data(iCase).pbeDefinition);
    % roundness measures
    [data(iCase).time, data(iCase).muB1] = extractDataFromPbeResult(...
        'bulkProperty', 6, ...
        data(iCase).result, data(iCase).pbeDefinition);
    [data(iCase).time, data(iCase).muB2] = extractDataFromPbeResult(...
        'bulkProperty', 7, ...
        data(iCase).result, data(iCase).pbeDefinition);
    [data(iCase).time, data(iCase).muB3] = extractDataFromPbeResult(...
        'bulkProperty', 8, ...
        data(iCase).result, data(iCase).pbeDefinition);
    % number
    [data(iCase).time, data(iCase).number] = extractDataFromPbeResult(...
        'bulkProperty', 5, ...
        data(iCase).result, data(iCase).pbeDefinition);
    % supersaturation
    [data(iCase).time, data(iCase).supersaturation] = extractDataFromPbeResult(...
        'bulkProperty', 1, ...
        data(iCase).result, data(iCase).pbeDefinition);
    % temperature
    [data(iCase).time, data(iCase).temperature] = extractDataFromPbeResult(...
        'bulkProperty', 11, ...
        data(iCase).result, data(iCase).pbeDefinition);
    data(iCase).timeMin = data(iCase).time/60;
    % large particle samples
    for iTime = 1:length(data(iCase).result)
        data(iCase).largeMuB1(iTime) = data(iCase).result(iTime).pivotPropertyMatrix(8, end);
        data(iCase).largeMuB3(iTime) = data(iCase).result(iTime).pivotPropertyMatrix(10, end);
        data(iCase).largeWidth(iTime) = data(iCase).result(iTime).pivotPropertyMatrix(7, end);
        data(iCase).largeVolume(iTime) = data(iCase).result(iTime).pivotPropertyMatrix(5, end);
        data(iCase).mediumMuB1(iTime) = data(iCase).result(iTime).pivotPropertyMatrix(8, round(0.9*end));
        data(iCase).mediumMuB3(iTime) = data(iCase).result(iTime).pivotPropertyMatrix(10, round(0.9*end));
        data(iCase).mediumWidth(iTime) = data(iCase).result(iTime).pivotPropertyMatrix(7, round(0.9*end));
        data(iCase).mediumVolume(iTime) = data(iCase).result(iTime).pivotPropertyMatrix(5, round(0.9*end));
        data(iCase).smallMuB1(iTime) = data(iCase).result(iTime).pivotPropertyMatrix(8, 1);
        data(iCase).smallMuB3(iTime) = data(iCase).result(iTime).pivotPropertyMatrix(10, 1);
        data(iCase).smallWidth(iTime) = data(iCase).result(iTime).pivotPropertyMatrix(7, 1);
        data(iCase).smallVolume(iTime) = data(iCase).result(iTime).pivotPropertyMatrix(5, 1);
    end
    % final size distribution
    delta = 1e-5;
    sizeRange = 0:delta:8e-5;
    sizeMatrix = [sizeRange; sizeRange + delta; sizeRange + 0.5*delta];
    [X Y] = extractDataFromPbeResult(...
        'propertyDistribution', {'volume', 'meanWidth', sizeMatrix}, ...
        data(iCase).result, data(iCase).pbeDefinition);
    data(iCase).finalDistX = sizeMatrix(3,:)
    data(iCase).finalDistY = Y(1:(end/2), end);
    
    % plot distribution
%     figure,plot(data(iCase).finalDistX,  data(iCase).finalDistY)
    
    % adjust muBx to mean (by number)
    data(iCase).muB1 = data(iCase).muB1 ./ data(iCase).number;
    data(iCase).muB2 = data(iCase).muB2 ./ data(iCase).number;
    data(iCase).muB3 = data(iCase).muB3 ./ data(iCase).number;
    data(iCase).filter = data(iCase).mw > 0;

    % plot overview
    figure,hold on,set(gcf, 'name', data(iCase).caseName)
    subplot(4, 1, 1)
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).temperature(data(iCase).filter))
    subplot(4, 1, 2)
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).supersaturation(data(iCase).filter))
    subplot(4, 1, 3)
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).mw(data(iCase).filter), 'k-')
    subplot(4, 1, 4)
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).vol(data(iCase).filter), 'k-')
     
    % plot particle roundness
    figure, hold on,title(data(iCase).caseName)
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).largeMuB1(data(iCase).filter), 'k--', 'lineWidth', 3)
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).largeMuB3(data(iCase).filter), 'k-', 'lineWidth', 3)
     
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).mediumMuB1(data(iCase).filter), 'k--', 'lineWidth', 2)
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).mediumMuB3(data(iCase).filter), 'k-', 'lineWidth', 2)
     
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).smallMuB1(data(iCase).filter), 'k--', 'lineWidth', 1)
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).smallMuB3(data(iCase).filter), 'k-', 'lineWidth', 1)

    % plot particle mean width
    figure, hold on,title(data(iCase).caseName)
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).largeWidth(data(iCase).filter), 'k-', 'lineWidth', 3)
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).mediumWidth(data(iCase).filter), 'k-', 'lineWidth', 2)
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).smallWidth(data(iCase).filter), 'k-', 'lineWidth', 1)
     
    % plot particle volume
%     figure, hold on,title(data(iCase).caseName)
%     plot(data(iCase).time(data(iCase).filter), ...
%          data(iCase).largeVolume(data(iCase).filter), 'k-', 'lineWidth', 3)
%     plot(data(iCase).time(data(iCase).filter), ...
%          data(iCase).mediumVolume(data(iCase).filter), 'k-', 'lineWidth', 2)
%     plot(data(iCase).time(data(iCase).filter), ...
%          data(iCase).smallVolume(data(iCase).filter), 'k-', 'lineWidth', 1)
    
    % plot some phase
    figure,hold on, title(data(iCase).caseName)
    color = [1 0 0];
    for iTime = [1:1:length(data(iCase).result) length(data(iCase).result)]
        if iTime == 1
            thisColor = [0 0 0];
        else
            thisColor = color*(1-iTime/length(data(iCase).result)) + ...
                       [1 1 1]*0.5*iTime/length(data(iCase).result);
        end
        plot(data(iCase).result(iTime).hcMatrix(1 ,:), ...
             data(iCase).result(iTime).hcMatrix(4 ,:), ...
             'color', thisColor, 'lineStyle', 'none', ...
             'marker', 'x')
    end
    % plot some phase
%     figure,hold on, title(data(iCase).caseName)
%     color = [1 0 0];
%     for iTime = [1:1:length(data(iCase).result) length(data(iCase).result)]
%         if iTime == 1
%             thisColor = [0 0 0];
%         else
%             thisColor = color*(1-iTime/length(data(iCase).result)) + ...
%                        [1 1 1]*0.5*iTime/length(data(iCase).result);
%         end
%         plot(data(iCase).result(iTime).pivotPropertyMatrix(4 ,:), ...
%              data(iCase).result(iTime).pivotPropertyMatrix(8 ,:), ...
%              'color', thisColor, 'lineStyle', 'none', ...
%              'marker', 'x')
%     end
end

save('data/growthSummary.mat', 'data')

