function pbeDefinition = lactoseGrowth
clear all
%% Basic information
pbeDefinition.info = ['seeded growth of lactose based on growth rate '...
    'data by Dincer, 2009'];
% The Miller indices are defined in the following order:
%   [0  1  0;  >> base              >> measured
%    0 -1  0;  >> tip               >> ASSUME NO GROWTH >> FIX
%    1  0  0;  >> side (middle)     >> measured
%    0 -1 -1;  >> large flat facet  >> measured
%    1  1  0;  >> side (at base)    >> measured
%   -1 -1  0]  >> side (at tip)     >> NOT MEASURED >> set to {???}

%% decomposition
fileNameDecomposition = 'lactoseGrowth_dec.matx';
if exist(fileNameDecomposition, 'file')
    dec = load(fileNameDecomposition, '-mat');
    dec = dec.dec;
else
    % Load case (contains source information)
    cdef = obtainCrystal('lactose');
    dec = obtainDec(cdef.A, 'mapping', cdef.groupMappingMatrix(:, [1 3:end]), ...
                    'allowFixPosition', 1);
    dec.cdef = cdef;
    % generate measure information
    dec = addMeasureDataDec(dec, 'volume');
    dec = addMeasureDataDec(dec, 'surface');
    dec = addMeasureDataDec(dec, 'meanWidth');
    dec = addMeasureDataDec(dec, 'support', dec.A([1 2 3 5 7 9], :));
    % save file
    save(fileNameDecomposition, '-mat', 'dec');
end
% copy to pbeDef
pbeDefinition.dec = dec;

%% bulk properties
% CONCENTRATION is defined by:
%    [mass of water-free lactose] / [100g water]

% solubility (factors by C.Choscz based on Dincer[2000])
factorF = @(tCelcius)(0.0016*tCelcius + 0.0048);
factorKm = @(tCelcius)(-0.002284*tCelcius + 1.627);
solubility = @(tCelcius)(sum(...
    [11.23, 0.3375, 1.236e-3, 172.57e-6, 0.5188e-6] .* ...
    tCelcius.^([0,1,2,3,4])));
saturation = @(C, tCelcius)(C ./ (solubility(tCelcius) - ...
    factorF(tCelcius)*factorKm(tCelcius)*(C - solubility(tCelcius)) ));

% total/initial masses
totalMassWater = 1;         % kg (1 liter)
totalMassLactose = 40*1e-2; % kg
initialCrystalMass = 1e-3;  % kg

% mass calculations
lactoseCrystalDensity = 1525; % kg/m^3 [Dincer, 2000]
molarMassWater = 18.01532;    % kg/mol
molarMassLactose = 342.2979;  % kg/mol
lactoseWaterMassFraction = molarMassWater/(molarMassLactose+molarMassWater);
lactoseMassFraction = 1-lactoseWaterMassFraction;
lactoseConcentration = @(V)(100*...
    (totalMassLactose - V*lactoseCrystalDensity*lactoseMassFraction) / ...
    (totalMassWater - V*lactoseCrystalDensity*lactoseWaterMassFraction));

%% initial condition
% NORMAL distribution for cuboid facets. The other facets are exectly
%   disappeared, giving some dirac functions that would be multiplied.
pbeDefinition.seeding.mean = [2 0.3 0.05 0.8 0.1]'*1e-5;
pbeDefinition.seeding.std = (pbeDefinition.seeding.mean*4e-2);
pbeDefinition.seeding.cov = diag(pbeDefinition.seeding.std.^2);
% pbeDefinition.seeding.totalValue = 1e10;
pbeDefinition.seeding.nSample = 1;

% we sample the cuboid coordinates directly..
rng(1)
disD = obtainDistributionDiscrete('normal', ...
    'nSample',pbeDefinition.seeding.nSample, ...
    'mean', pbeDefinition.seeding.mean, ...
    'cov', pbeDefinition.seeding.cov);
% compute total volume
% max(max(dec.validityCone.A * disD.positionMatrix))
sampleVolume = sum(computeMeasureDec(dec, 'hc', disD.positionMatrix, 'volume'));
pbeDefinition.seeding.positionMatrix = disD.positionMatrix;
pbeDefinition.seeding.numberVector = ones(1, pbeDefinition.seeding.nSample) * ...
    (initialCrystalMass / sampleVolume / lactoseCrystalDensity);
% sum(pbeDefinition.seeding.numberVector .* ...
%     computeMeasureDec(dec, 'hc', pbeDefinition.seeding.positionMatrix, 'volume'))

%% growth and analytical solution
% taken from Dincer, 2009 (without dead zone) in the order:
%  {0 1 0} (1), {1 1 0} (4), {1 0 0} (2) , {0 -1 -1} (3)
% missing is then {-1 -1 0} (5) set equal to {
n = [2.6 2.8 3.08 2.5]';                             % exponent
k = [4.2*10^5 4.2*10^5 2.6*10^8 1.0*10^5]'*1e-6/60;  % factor m/s
E = [9.5 9.7 13.7 10.2]'*4186.8;                     % energy J/mol
R = 8.314472;                                        % J/mol/K
growthRate = @(s, tKelvin)(k.*exp(-1*E./R./tKelvin).*abs(s-1).^n*sign(s-1));
mapGrowthRate = [1 0 0 0;...
                 0 0 1 0;...
                 0 0 0 1;...
                 0 1 0 0;...
                 0 0 1 0];
coolingRateOne = 1/60;  % K/s
coolingRateTwo = 0.005/60; % K/s
temperature = @(t)(273.15+40 - ...
    min(t,10/coolingRateOne)*(coolingRateOne-coolingRateTwo) - ...
    min(t,10/coolingRateTwo)*coolingRateTwo);
% temperature = @(t)(273.15+40 - ...
%     min(t,10/coolingRateOne)*coolingRateOne);

%% function handles
% we need the unified partitions in P for growth rate projections
pbeDefinition.pivotProperty = @(H, up)([...
    up;
    computeMeasureDec(dec, 'hc', H, 'volume')]);
% bulk property computes:
%  1 - total crystal volume (m^3)
%  2 - temperature (Kelvin)
%  3 - concentration
%  4 - saturation
pbeDefinition.bulkProperty = @(P, N, x, t)(...
    [N*P(2,:)'; ...
     temperature(t); ...
     lactoseConcentration(N*P(2,:)'); ...
     saturation(lactoseConcentration(N*P(2,:)'),...
                temperature(t)-273.15);...
     ]);

pbeDefinition.growthRate = @(H, P, y, t)(...
    (mapGrowthRate * growthRate(y(4), y(2))) * ones(1, size(H,2)));

%% time span for solution
% pbeDefinition.timeSpan = [0:60:(20*1e3) (20*1e3+300):300:1e5] ; %s
pbeDefinition.timeSpan = 0:120:(2*24*60*60);
% pbeDefinition.timeSpan = 0:20:100;

% options
pbeDefinition.optionStruct.distanceLimit = 1e-10; %m
pbeDefinition.optionStruct.relTol = 1e-5;
pbeDefinition.optionStruct.absTolNumber = 1e-30;
pbeDefinition.optionStruct.absTolPosition = 1e-30;
% pbeDefinition.optionStruct.odeSolver = @ode113; % 113, 45