function pbeDefinition = commonPbeDefinitionOct(par)

%% Name and information
pbeDefinition.name = 'commonPbeDefinition';
pbeDefinition.info = ['potash alum grows with simplified kinetics and '...
    'constant attrition is applied'];
pbeDefinition.source = ['Ma et al. 2008: "Morphological Population ' ...
    'Balance for Modeling Crystal Growth in Face Directions"'];

%% From the source (growth) 
% concentrations: mol/m^3
% supersaturation: c/c* - 1
% we adopt here to use only the octahedron kernel, because growth already
%   gives the right behavior for de-rounding since the facets would not be
%   present, otherwise. The attrition however does not consider small
%   facets as almost-round parts of the crystal.

%% sensible parameters

% the stirrer speed is set here and some (arbirary) kinetic is scaled to that,
%   afterwards
parStandard.stirrerSpeed = 750 / 60; % turns per second (scaled from rpm)
% linking factors for the rounding effect of attrition
parStandard.kAttrition = 6.5;% high implies fast round-off
parStandard.qGrowth = 1.4;   % high implies fasst grow out of rounding
% time span
parStandard.timeSpan = [0:20 40:20:800]*60       %s
if ~exist('par', 'var') || isempty(par)
    par = parStandard;
else
    par = mergeStruct(par, parStandard);
end
stirrerSpeed = par.stirrerSpeed;
kAttrition = par.kAttrition;
qGrowth = par.qGrowth;
pbeDefinition.timeSpan = par.timeSpan      %s
clear par parStandard

% attrition rate (fixed to match experiments well)
physicalRateAttrition = -1*stirrerSpeed / (500/60) * 7*10^-9 % m/s


%% vessel
temperatureCelcius = 30;                         % °C
temperatureKelvin = temperatureCelcius + 273.15; % K
% if the solubility from the main source is compared to data by Nollet
%   (2006), then it becomes evident that the equation is for mol per 100g
%   H2O water. Here, we use mol/m^3.
solubility = 10^5*exp(12.1/temperatureKelvin + 10.47*log(temperatureKelvin) - 65.73); % mol/m^3
initialSoluteConcentration = solubility*1.00 % mol/m^3

% the vessel volume is assumed to be constant, even under crystallization
%   and temperature effects
vesselVolume = 1 / 10^3; % m^3 (scaled from liters)

%% crystal
pbeDefinition.dec = obtainCrystalDec('potash_alum_Ma');
cdef = pbeDefinition.dec.cdef;
pbeDefinition.dec = obtainCrystalDec('generic_potash_alum', 0);

% decomposition data
% pbeDefinition.dec = obtainDec(pbeDefinition.cdef.A, 'mapping', ...
%     pbeDefinition.cdef.groupMappingMatrix, 'verbose', 0);
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'volume');
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'surface');
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'mean width');

%% initial distribution
% coordinates 1-3: kernel crystal
% coordinate 4: lambda

% get a distribution
rng(1);
meanSize = 36*[1 1.45]'*1e-6;   % m
stdSize = 0.0005*meanSize; % m
% disD = obtainDistributionDiscrete('normal', ...
%     'nSample', 10, ...
%     'mean', meanSize, ...  %m
%     'cov', diag(stdSize.^2));
% assign to pbeDefinition
pbeDefinition.nDecCoordinate = 1;
pbeDefinition.nExtraCoordinate = 1;
pbeDefinition.seeding.time = 0; % seconds
% pbeDefinition.seeding.positionMatrix = disD.positionMatrix;       % m
pbeDefinition.seeding.positionMatrix = meanSize;
% pbeDefinition.seeding.probabilityVector = disD.probabilityVector;
% pbeDefinition.seeding.numberVector = disD.probabilityVector*0 + ...
%     1e6 / numel(disD.probabilityVector); % #
pbeDefinition.seeding.numberVector = 1e6;

initialPositionMatrix = pbeDefinition.seeding.positionMatrix

%% some intermediate functions
% .. mixed with initial value computation / verification
% we already define the pivot properties:
% 1 - hC1
% 2 - lambda
% 3 - kernel crystal volume
% 4 - kernel crystal surface
% 5 - kernel crystal mean width
% 6 - mean width of one based on kernel partition
% 7 - mean width
% 8 - relative additive roundness
% 9 - unified partition for kernel crystal


% base distribution functional
pbeDefinition.pivotProperty = @(H, up)(...
    [H(1:2, :); ...
     computeMeasureDec(pbeDefinition.dec, 'hC', {up, H(1, :)}, {'volume', 'surface', 'meanWidth'}, [], 0); ...
     computeMeasureDec(pbeDefinition.dec, 'hC', {up, H(1, :)*0+1}, {'meanWidth'}, [], 0); ...
     2*H(2,:) + computeMeasureDec(pbeDefinition.dec, 'hC', {up, H(1, :)}, {'meanWidth'}, [], 0); ...
     2*H(2,:) ./ (...
     2*H(2,:) + computeMeasureDec(pbeDefinition.dec, 'hC', {up, H(1, :)}, {'meanWidth'}, [], 0)); ...
     up ]);
 
% ..and their initial values
initialP = pbeDefinition.pivotProperty(...
    pbeDefinition.seeding.positionMatrix, ...
    computePartitionIndexDec(pbeDefinition.dec, ...
    pbeDefinition.seeding.positionMatrix(1:(end-1), :), elkZerotol, 1));
initialN = pbeDefinition.seeding.numberVector

% crystal volume functionals
pivotCrystalVolume = @(P)(P(3,:) + ...
    P(4, :) .* P(2, :) + P(5, :).*P(2, :).^2 + ...
    4/3*pi*P(2, :).^3);
pivotCrystalVolume(initialP)
totalCrystalVolume = @(P, N)(pivotCrystalVolume(P)*N');
% ..and their initial values
initialCrystalVolume = totalCrystalVolume(initialP, initialN)
initialCrystalVolumeFraction = initialCrystalVolume / vesselVolume

% suspension functionals
soluteConcentration = @(P, N)(...
    (initialSoluteConcentration * vesselVolume + ...
     cdef.density/cdef.molarMass*...
       (initialCrystalVolume - totalCrystalVolume(P, N))...
    ) / vesselVolume);
supersaturation = @(P, N)(soluteConcentration(P, N)/solubility - 1);
% ..and their initial values
% initialSoluteConcentration = soluteConcentration(initialP, initialN)
initialSupersaturation = supersaturation(initialP, initialN)

% ..and their initial values
initialMeanWidth = initialP(7, :) * initialN' / sum(initialN)
initialMeanRelativeRoundness = initialP(8, :) * initialN' / sum(initialN)

% initial bulk states
bulkInitial = []; % concentration in mol/m^3

%% rates: physical base
% crystal growth: k1*s^k2
% ratio to fastest growing: 2.2495
% ration to second fastest: 1.4498
physicalRateGrowth = @(P, super, t)(7.753 * 10^-7 * super^1.5 * ...
                       ones(1, size(P, 2)) * (super>0)); % m/s

%% rates: attrition
qRounding = @(muB)(kAttrition*(1-muB));

rateRoundingConstH = @(P)(...
    [-1; 1] * (...
    qRounding(P(8, :)) ./ (1 + qRounding(P(8, :))) .* ...
    physicalRateAttrition ./ (2 - P(6, :)))...
    );

rateRoundingConstShape = @(P)(...
    (ones(2, 1) * ...
     (1 ./ (1 + qRounding(P(8, :))) .* ...
      physicalRateAttrition ./ (2 + P(5, :)./P(2,:)) ...
    )) .* ...
    [P(1,:)./(ones(1,1)*P(2,:)); ones(1, size(P, 2))]);

%% rates: growth
rateGrowth = @(P, s, t)(...
    [physicalRateGrowth(P, s, t); zeros(1, size(P, 2))] + ...
    [1; -1]*...
    (qGrowth ./ (P(6,:) - 2) .* ...
     computeMeasureDec(pbeDefinition.dec, 'hC', ...
      {P(9,:), physicalRateGrowth(P, s, t)}, {'meanWidth'}, [], 0) ...
    ));

%% totalRate

rateTotal = @(P, p, t)(...
    rateRoundingConstH(P) + rateRoundingConstShape(P) + ...
    rateGrowth(P, p(1), t));

pbeDefinition.growthRate = @(H, P, y, t)(...
    rateTotal(P, y, t) .* ...
    [ones(1, size(P, 2)); ...
     ([0 1] * rateTotal(P, y, t) > 0) | ...
     (P(2,:) > 1e-14)] ...
    );

%% bulk properties
% 1 - supersaturation
% 2 - total volume
% 3 - solute concentration
% 4 - mean width (mean)
% 5 - relative additive roundness (mean)
pbeDefinition.bulkProperty = @(P, N, x, t)([...
    supersaturation(P, N); ...       %1
    totalCrystalVolume(P, N); ...    %2
    soluteConcentration(P, N); ...   %3
    P(7, :) * N' / sum(N); ... 
    P(8, :) * N'/ sum(N) ...
    ]);   

% test fucntion
initialBulkProperty = pbeDefinition.bulkProperty(initialP, initialN, [])
% lambdaRateFactor(initialP(5, :))
% lambdaRate(initialSupersaturation, initialP(5, :), 0)
% crystalRate(initialSupersaturation, initialP(5, :), 0)
initialGrowthRate = pbeDefinition.growthRate([], initialP, initialBulkProperty, 0)


pbeDefinition.bulkRate = @(y, t)([]);
pbeDefinition.bulkInitial = bulkInitial;