clear all
%% generate or load results
tmp = executeCase('init');
caseList = {'linearGrowth_metroNumber', 'linearGrowth_metroVolume', ...
    'linearGrowth_metroMix', 'linearGrowth_highRes'};
nameList = {'MetroNum', 'MetroVol', 'MetroMix', 'HighRes'};
for iCase = 1:length(caseList)
    thisCaseName = [caseList{iCase}]
    try
        [result pbeDef] = loadCase(thisCaseName);
    catch
        pbeDef = obtainData(thisCaseName);
        [result pbeDef] = executeCase(thisCaseName);
    end
    eval(['result' nameList{iCase} '=result;'])
    eval(['pbeDef' nameList{iCase} '=pbeDef;'])
end

%% analytic results (must be generated beforehand)
% anaStruct = prepareConstGrowth(1, 'volume');

%% selection
% result = resultMetroNum;pbeDef = pbeDefMetroNum;
% result = resultMetroVol;pbeDef = pbeDefMetroVol;
result = resultMetroMix;pbeDef = pbeDefMetroMix;
result = resultHighRes;pbeDef = pbeDefHighRes;

%% prepare data
% numerical data
[time, Y] = extractDataFromPbeResult('integralProperty', 'number', ...
    result, pbeDef);
numberNum = Y(1,:);numberStd = Y(2,:);
[time, Y] = extractDataFromPbeResult('integralProperty', 'volume', ...
    result, pbeDef);
volumeNum = Y(1,:);volumeStd = Y(2,:);
[time, Y] = extractDataFromPbeResult('integralProperty', 'surface', ...
    result, pbeDef);
surfaceNum = Y(1,:);surfaceStd = Y(2,:);
[time, Y] = extractDataFromPbeResult('integralProperty', 'meanWidth', ...
    result, pbeDef);
meanWidthNum = Y(1,:);meanWidthStd = Y(2,:);

numberAna = pbeDef.analyticNumber(time);
volumeAna = pbeDef.analyticVolume(time);
surfaceAna = pbeDef.analyticSurface(time);
meanWidthAna = pbeDef.analyticMeanWidth(time);

%% plot values
figure,hold on
plot(time, numberAna, 'k-.')
plot(time, volumeAna, 'k-o')
plot(time, surfaceAna, 'k-^')
plot(time, meanWidthAna, 'k-+')
% plot(time, log(numberAna)./pbeDef.growthConstant, 'k-.')
% plot(time, log(volumeAna)./pbeDef.growthConstant, 'k-o')
% plot(time, log(surfaceAna)./pbeDef.growthConstant, 'k-^')
% plot(time, log(meanWidthAna)./pbeDef.growthConstant, 'k-+')

ylim([1e-5 2e10])
set(gca, 'yscale', 'log')
% set(gca, 'xscale', 'log')
xlabel('time in s')
ylabel('total volume in m^3')

%% error plot
figure,hold on
plot(time, numberStd./numberNum, 'k--.')
plot(time, volumeStd./volumeNum, 'k--o')
plot(time, surfaceStd./surfaceNum, 'k--^')
plot(time, meanWidthStd./meanWidthNum, 'k--+')

plot(time, abs(numberNum./numberAna - 1), 'k-.')
plot(time, abs(volumeNum./volumeAna - 1), 'k-o')
plot(time, abs(surfaceNum./surfaceAna - 1), 'k-^')
plot(time, abs(meanWidthNum./meanWidthAna - 1), 'k-+')

numberCVar = max(numberStd./numberNum)
volumeCVar = max(volumeStd./volumeNum)
surfaceCVar = max(surfaceStd./surfaceNum)
meanWidthCVar = max(meanWidthStd./meanWidthNum)

numberRelErr = max(abs(numberNum./numberAna - 1))
volumeRelErr = max(abs(volumeNum./volumeAna - 1))
surfaceRelErr = max(abs(surfaceNum./surfaceAna - 1))
meanWidthRelErr = max(abs(meanWidthNum./meanWidthAna - 1))

xlabel('time in s')
ylabel('error measures')

%% density error
errorVector = time*nan;
for iTime = 1:length(time)
    anaDensity = pbeDef.analyticSolution(result(iTime).hcMatrix, time(iTime));
    numDensity = length(result(iTime).numberVector)*...
        result(iTime).numberVector .*...
        result(iTime).probabilityVector;
    errorVector(iTime) = (anaDensity-numDensity)/anaDensity;
end
figure
plot(time, errorVector)



% axisHandle = axes;hold on
% plot(time, abs(volumeNum./volumeAna - 1), 'r--d', 'parent', axisHandle)
% % plot([0 X(end)], [0 0], 'm-', 'parent', axisHandle, 'lineWidth', 3)
% set(axisHandle, 'color', 'none', 'yaxisLocation', 'right', ...
%     'xTick', [])
% % legend('coef. of variation (right axis)', 'Location', 'NorthEast')
% ylabel('coefficient of variation')


