function pbeDefinition = common(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

%% 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

% this value pretty much seems to increase mean width growth with small
% values
parStandard.kqAttrition = 140;% high implies fast round-off
parStandard.kAttrition = 1*10^-10;
parStandard.qGrowth = 30;   % high implies fast grow out of rounding
% time span
parStandard.timeSpan = [0:20 40:5:1600]*60       %s
if ~exist('par', 'var') || isempty(par)
    par = parStandard;
else
    par = mergeStruct(par, parStandard);
end
stirrerSpeed = par.stirrerSpeed;
kqAttrition = par.kqAttrition;
kAttrition = par.kAttrition;
qGrowth = par.qGrowth;
pbeDefinition.timeSpan = par.timeSpan      %s
clear par parStandard


%% 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

% 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 = 100*[1 1.45 1.1 0.2]'*1e-6;   % m
stdSize = 0.005*meanSize; % m
disD = obtainDistributionDiscrete('normal', ...
    'nSample', 10, ...
    'mean', meanSize, ...  %m
    'cov', diag(stdSize.^2));
% assign to pbeDefinition
pbeDefinition.nDecCoordinate = 3;
pbeDefinition.nExtraCoordinate = 1;
pbeDefinition.seeding.time = 0; % seconds
% pbeDefinition.seeding.positionMatrix = disD.positionMatrix;       % m
% pbeDefinition.seeding.positionMatrix = [0.5*meanSize meanSize];
pbeDefinition.seeding.positionMatrix = meanSize * [0.10:0.05:1];
% pbeDefinition.seeding.probabilityVector = disD.probabilityVector;
% pbeDefinition.seeding.numberVector = disD.probabilityVector*0 + ...
%     3e6 / numel(disD.probabilityVector); % #
pbeDefinition.seeding.numberVector = 5e2*...
    pbeDefinition.seeding.positionMatrix(1,:).^-1/...
    size(pbeDefinition.seeding.positionMatrix, 2); % #
% pbeDefinition.seeding.numberVector = [1e7 4e6];

initialPositionMatrix = pbeDefinition.seeding.positionMatrix

%% pivot measures
kernelVolume = @(H, up)(computeMeasureDec(pbeDefinition.dec, 'hC', ...
    {up, H(1:3, :)}, 'volume', [], 0));
kernelSurface = @(H, up)(computeMeasureDec(pbeDefinition.dec, 'hC', ...
    {up, H(1:3, :)}, 'surface', [], 0));
kernelMeanWidth = @(H, up)(computeMeasureDec(pbeDefinition.dec, 'hC', ...
    {up, H(1:3, :)}, 'meanWidth', [], 0));
particleVolume = @(H, up)(...
    kernelVolume(H, up) + ...
    kernelSurface(H, up) .* H(4, :) + ...
    2*pi*kernelMeanWidth(H, up).*H(4, :).^2 + ...
    4/3*pi*H(4, :).^3);


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

% base distribution functional
pbeDefinition.pivotProperty = @(H, up)(...
    [H(1:4, :); ...
     kernelVolume(H, up)
     kernelSurface(H, up)
     kernelMeanWidth(H, up)
     computeMeasureDec(pbeDefinition.dec, 'hC', {up, H(1:3, :)*0+1}, {'meanWidth'}, [], 0); ...
     2*H(4,:) + kernelMeanWidth(H, up); ...
     2*H(4,:) ./ (2*H(4,:) + kernelMeanWidth(H, up)); ...
     up;...
     particleVolume(H, 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(5,:) + ...
    P(6, :) .* P(4, :) + ...
    2*pi*P(7, :).*P(4, :).^2 + ...
    4/3*pi*P(4, :).^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(9, :) * initialN' / sum(initialN)
initialMeanRelativeRoundness = initialP(10, :) * initialN' / sum(initialN)

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

%% rates: physical base
% crystal growth: k1*s^k2
physicalRateGrowth = @(P, super, t)([7.753 * 10^-7 * super^1.5; ...
                       1.744 * 10^-6 * super^1.5; ...
                       1.124 * 10^-6 * super^1.5] * ...
                       ones(1, size(P, 2)) * (super>0)); % m/s
% attrition rate (fixed to match experiments well)
physicalRateAttrition = -1*stirrerSpeed *kAttrition % m/s

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

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

rateRoundingConstShape = @(P)(...
    (ones(4, 1) * ...
     (1 ./ (1 + qRounding(P(10, :))) .* ...
      physicalRateAttrition ./ (2 + P(7, :)./P(4,:)) ...
    )) .* ...
    [P(1:3,:)./(ones(3,1)*P(4,:)); ones(1, size(P, 2))]);

%% rates: growth
rateGrowth = @(P, s, t)(...
    [physicalRateGrowth(P, s, t); zeros(1, size(P, 2))] + ...
    [1; 1; 1; -1]*...
    (qGrowth ./ (P(8,:) - 2) .* ...
     computeMeasureDec(pbeDefinition.dec, 'hC', ...
      {P(11,:), 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(3, size(P, 2)); ...
     ([0 0 0 1] * rateTotal(P, y, t) > 0) | ...
     (P(4,:) > 1e-14)] ...
    );

%% bulk properties
% 1 - supersaturation
% 2 - total volume
% 3 - solute concentration
% 4 - mean width
% 5 - relative additive roundness
% 6 - "mean" unified partition
% 7 - min distance to validity cone boundary
% 8 - total number
pbeDefinition.bulkProperty = @(P, N, x, t)([...
    supersaturation(P, N); ...       %1
    P(12, :) * N'; ...    %2
    soluteConcentration(P, N); ...   %3
    P(9, :) * N' / sum(N); ... 
    P(10, :) * N'/ sum(N); ...
    P(11,:) * N'/ sum(N); ...
    max(max(pbeDefinition.dec.validityCone.A * P(1:3,:))); ...
    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;


pbeDefinition.optionStruct.relTol = 1e-5;
pbeDefinition.optionStruct.odeSolver = @ode15s;