function pbeDefinition = caseParacetamolSingleCharacteristic

%% description
% This case intends to reproduce the case study by Ma and Wang in
%   "Morphological Population Balance Model in Principal Component Space".
%   Explicit references are given for almost al parameters that are used,
%   here.
pbeDefinition.name = 'paracetamolSingleCharacteristic';
pbeDefinition.source = ['Borchert and Sundmacher 2012: "Efficient ' ...
    'formulation of crystal shape evolution equations"'];

%% crystal definition
% the crystal definition already provides molar mass and density as a
%   property of the crystal compound
pbeDefinition.cdef = obtainCrystal('paracetamol');
cdef = pbeDefinition.cdef;

%% growth rate parameters
% The following coefficients are required for the growth rate model.
%   Overall, it is a piecewise polynomial model of the supersaturation so 
%   that each vector represents the coefficients: first entry (power 1),
%   second entry (power 2) and so on. Power 0 does not exist.
% The vectors are defined for regions of supersaturation, each so that the
%   naming is g_xyz_i for {xyz} representing the growing crystal form and i
%   the region to use. The upper bound for each region is ub_xyz_i.

% form 1 (Borchert) and form 1 in library
g_110_1 = 1.041e-8;
ub_110_1 = 0.1076;
g_110_2 = [-4.336e-10 -2.142e-7 2.125e-6];

% form 2 (Borchert) and form 4 in library
g_201_1 = [-1.131e-9 1.847e-6 -5.883e-6];

% form 3 (Borchert) and form 2 in library
g_011_1 = [-5.940 -7.461 -2.238];

% form 4 (Borchert) and form 3 in library
g_001_1 =[2.431e-8 -1.272e-8];
ub_001__1 = 0.0533;
g_001_2 = [-5.975e-8 1.867e-6 -5.699e-6];


%% growth rate
growthRate = @(sigma, T)(...
    bsxfun(@times, growthConstant*0+1, ...
      bsxfun(@power, ... 
        [sigma + (P(2, :)*0); ...
         sigmaEff(sigma, ...
                  min([P(2,:); P(3,:)], [], 1)); ...
         sigmaEff(sigma, ...
                  min([P(3,:); P(4,:); P(5,:)], [], 1))], ...
        growthExponent*0+1 ...   
      ) ...
    ));
% verify implementation:
% growthRate(0.0659, [nan nan nan; ...
%                  1   1  4; ...
%                  1   3  3; ...
%                  1   3  4; ...
%                  1   3  4]*1e-6)

% pbe growth rate
pbeDefinition.growthRate = @(H, P, x, y, t)(...
    growthRate(y(1), P));   % m/s


%% decomposition data and polytope properties
% decomposition data
pbeDefinition.dec = obtainDec(pbeDefinition.cdef.A, 'mapping', ...
    pbeDefinition.cdef.groupMappingMatrix, 'verbose', 0);
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'volume');
% edge 1 - 100/111
% edge 2 - 100/110 = 110/100
% edge 3 - 110/110
% edge 4 - 110/111
% edge 100/100 does not exist as such (equals 100/110)
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'edgeLength', ...
    [3 1 1; 3 1 0; 2 3 1; 3 3 1]);

% pivot properties
pbeDefinition.pivotProperty = @(H, up)(...
    computeMeasureDec(pbeDefinition.dec, 'hC', {up, H}, ...
      {'volume', 'edgeLength', [3 1 1; 3 1 0; 2 3 1; 3 3 1]}, ...
      [], 0));

%% bulk states and bulk properties
% bulk state
pbeDefinition.bulkRate = @(x, y, t)([]);
pbeDefinition.bulkInitial = [];

% bulk properties
% 1) supersaturation
pbeDefinition.bulkProperty = @(P, N, x)([...
    P(1,:)*N' ...
    ]);

%% initial condition
% initial distribution: defined as in Ma (2008, appendix) and verified
%   against figure 5 in Wang (2009)
% Note that the standard deviation is such that 2% of the population is
%   non-existent if only one coordinate is observed and the geometrical
%   constraints are ignored. A much larger part of this population will
%   reside at the boundary of the valid cone.
par.mean = [9 8 8.5]'*1e-6;  % m
par.std = [4 4 4]'*1e-6;     % m
par.totalValue = 2200;       % #
disC = obtainDistributionContinuous('normal', par);

%! ToDo: to be done; best is to enforce a regular grid, here.
% pbeDefinition.seeding{1} = struct(...
%     'time', 0, ...      % s
%     'disC', disC);
% pbeDefinition.seeding{1}.optionStructDiscretization = struct(...
%     'nPivot', 1e3, ...
%     'method', 'box uniform', ...
%     'regionMin', disC.par.mean - 0.1*disC.par.std, ...
%     'regionMax', disC.par.mean + 0.1*disC.par.std);
% pbeDefinition.seeding{1}.optionStructReduction = struct(...
%     'nPivot', 1e2, ...
%     'nRerun', 10);


%% time span for solution
% Wang 2009, state 57 time steps which gives an odd spacing of 58.9286
%   seconds, while a reasonable spacing of 60 seconds gives 56 time steps.
%   Hence, a 60 seconds spacing with 56 time steps is chosen.
pbeDefinition.timeSpan = 0:60:3300;       % s