function pbeDefinition = paracetamolGrowth
% clear all
%% Basic information
pbeDefinition.info = 'seeded growth of paracetamol based on Borchert 2009';
% The Miller indices are defined in the following order:
%   [1  1  0;  >> four facets orthogongal to xy-plane
%    0  1  1;  >> four facets diagonal
%    0  0  1;  >> two facets between {0 1 1}
%    2  0 -1]  >> the other two facets
% In contrast, Borchert uses the following order:
%   {110} (1), {20-1} (4), {011} (2), {001} (3) so that we apply a mapping
%   matrix from the Borchert order to the order applied here:
mapGrowthRate = [1 0 0 0;...
                 0 0 1 0;...
                 0 0 0 1;...
                 0 1 0 0];
% mapGrowthRate = [1/sqrt(2) 0 0 0;...
%                  0 0 1/sqrt(3) 0;...
%                  0 0 0 1/sqrt(2);...
%                  0 1 0 0];

%% decomposition
fileNameDecomposition = 'paracetamolGrowth_dec.matx';
if exist(fileNameDecomposition, 'file')
    dec = load(fileNameDecomposition, '-mat');
    dec = dec.dec;
else
    % Load case (contains source information)
    dec = obtainCrystalDec('generic_paracetamol',1);
    % 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 mole or mass fraction:
%   [kg paracetamol] / [kg paracetamol + kg solvent]
tStart = 40;
paracetamolCrystalDensity = 1.263*1e3/1e6; % kg/m^3
molarMassParacetamol = 151.17;             % kg/kmol [Borchert, 2012]
molarMassWater = 18.02;                    % kg/kmol [Borchert, 2012]

solubilityParameters = [-298.59288 10495.9 45.11344];
solubility = @(tKelvin)(exp(...
    solubilityParameters * [1 1/tKelvin log(tKelvin)]'));
saturation = @(xPara, tKelvin)(xPara / solubility(tKelvin) - 1);

% total/initial masses
totalMassWater = 1;         % kg (1 liter)
initialCrystalMass = 1e-3;  % kg
totalMassParacetamol = initialCrystalMass + ...
    solubility(273.15+tStart)/(1-solubility(273.15+tStart))*...
    molarMassParacetamol/molarMassWater*totalMassWater; % kg

% mass calculations
paracetamolConcentration = @(V)(...
    (totalMassParacetamol - V*paracetamolCrystalDensity) / ...
    (totalMassParacetamol - V*paracetamolCrystalDensity + ...
     molarMassParacetamol/molarMassWater*totalMassWater));
 
% s = x/xSat - 1
% x = mPar / (mPar + molarPar/molarWat*mWat)
%
% s = 0 >> x = xSat
% mPar = xSat * (mPar + molarPar/molarWat*mWat)
% mPar = xSat*molarPar/molarWat*mWat/(1-xSat)

%% initial condition
pbeDefinition.seeding.mean = [1 1 1 1]'*5e-5;
pbeDefinition.seeding.std = (pbeDefinition.seeding.mean*5e-2);
pbeDefinition.seeding.cov = diag(pbeDefinition.seeding.std.^2);
pbeDefinition.seeding.nSample = 1000;
rng(1)
disD = obtainDistributionDiscrete('normal', ...
    'nSample',pbeDefinition.seeding.nSample, ...
    'mean', pbeDefinition.seeding.mean, ...
    'cov', pbeDefinition.seeding.cov);
% isValidNumber = sum(isValidDec(disD.positionMatrix, pbeDefinition.dec))

% 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 / paracetamolCrystalDensity);
% pbeDefinition.seeding.numberVector
% sum(pbeDefinition.seeding.numberVector .* ...
%     computeMeasureDec(pbeDefinition.dec, 'hc', pbeDefinition.seeding.positionMatrix, 'volume'))

%% growth and analytical solution
g1_bnd = 0.1076;
g1_lt = [1.041e-8 0];
g1_ge = [2.125e-6 -2.142e-7 -4.366e-10];
g2 = [-5.883e-6 1.847e-6 -1.131e-9 0];
g3 = [-2.238e-6 7.461e-7, 5.940e-8 0];
g4_bnd = 0.0533;
g4_lt = [-1.272e-8 2.431e-8 0];
g4_ge = [-5.699e-6 1.867e-6 -5.975e-8 0];

growthRate = @(s)(sign(s)*[...
    (s <  g1_bnd)*(abs(s).^((length(g1_lt)-1):-1:0)*g1_lt') + ...
    (s >= g1_bnd)*(abs(s).^((length(g1_ge)-1):-1:0)*g1_ge'); ...
                  (abs(s).^((length(g2)-1):-1:0)*g2'); ...
                  (abs(s).^((length(g3)-1):-1:0)*g3'); ...
    (s <  g4_bnd)*(abs(s).^((length(g4_lt)-1):-1:0)*g4_lt') + ...
    (s >= g4_bnd)*(abs(s).^((length(g4_ge)-1):-1:0)*g4_ge')]);


% THIS IS DONE TO PROVOKE GROWING OUT FACETS:
% coolingRateOne = 0.5/60;  % K/s
% coolingRateTwo = 0.3*0.13/60; % K/s
% temperature = @(t)(273.15+tStart - ...
%     min(t,3/coolingRateOne)*(coolingRateOne-coolingRateTwo) - ...
%     min(t,17/coolingRateTwo)*coolingRateTwo);
%/
coolingRate = 0.15/60; % K/s
temperature = @(t)(273.15+tStart - min(t,10/coolingRate)*coolingRate);

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

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

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

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