%% Snyder Cuboid Dissolution
% This one runs all simulations and stores the results (and shows some
% figure).
clc
dynamicCaseName = 'ureaGenericGrowth';

%% Run original simulation(s)
% run simulation
executeLocalPbeCase(dynamicCaseName);
% run simulation without limits
executeLocalPbeCase([dynamicCaseName '_noGrowthLimit']);

% load data (no limits)
[result pbeDef] = loadLocalPbeResult([dynamicCaseName '_noGrowthLimit']);

%% Generate (handmade) approximation
spanMatrix = [pbeDef.growthRate(ones(7,1), [], [], []) ...
              pbeDef.seeding.positionMatrix];
% compute a basis for the range of spanMatrix
spanMatrix = orth(spanMatrix);
% generate data file
thisFile = ['./database/' dynamicCaseName '_4_handmade.matx'];
if ~exist(thisFile, 'file')
    dec = projectDec(pbeDef.dec, spanMatrix, 'zerotol', 1e-8);
    dataApprox.dec = dec;
    dataApprox.cdef = pbeDef.dec.cdef;
    save(thisFile, '-mat', 'dataApprox');
% else
%     data = load(thisFile, '-mat');
end

%% Run approximation simulation
thisCase = [dynamicCaseName '_4_handmade'];
executeLocalPbeCase(thisCase)
[resultTest pbeDefTest] = loadLocalPbeResult(thisCase);

%% Print some figures

% number over time
[X, Y] = extractDataFromPbeResult('integralProperty', 'number', ...
    result, pbeDef);
[X, Yapprox] = extractDataFromPbeResult('integralProperty', 'number', ...
    resultTest, pbeDefTest);
figure,hold on
plot(X, Y(1,:), 'b-x');
plot(X, Yapprox(1,:), 'r.-');
title('number')

% meanWidth over time
[X, Y] = extractDataFromPbeResult('integralProperty', 'meanWidth', ...
    result, pbeDef);
[X, Yapprox] = extractDataFromPbeResult('integralProperty', 'meanWidth', ...
    resultTest, pbeDefTest);
figure,hold on
plot(X, Y(1,:), 'b-x');
plot(X, Yapprox(1,:), 'r.-');
% plot(X, Y(1,:) + Y(2,:), 'r.');
ylabel('meanWidth')
figure,plot(X, Y(1,:)-Yapprox(1,:))