function pbeDefinition = commonExtraMassStorage(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
parStandard.kqAttrition = 500;% high implies fast round-off
parStandard.kAttrition = 10^-10;
parStandard.qGrowth = 1.2;   % high implies fasst grow out of rounding
% rate for the "small particle" storage
parStandard.kSmallParticle = 2e-4; % 1/s
% 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;
kqAttrition = par.kqAttrition;
kAttrition = par.kAttrition;
qGrowth = par.qGrowth;
kSmallParticle = par.kSmallParticle;
pbeDefinition.timeSpan = par.timeSpan      %s
clear par parStandard

%% 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');

%% vessel
% temperature
temperatureCelcius = 30;                         % °C
temperatureKelvin = temperatureCelcius + 273.15; % K

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

% 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
initialSoluteMass = initialSoluteConcentration * cdef.molarMass * vesselVolume
% initialExtraMolConcentration = solubility*0.05 % mol/m^3
% initialExtraMassConcentration = cdef.molarMass * ...
%     initialExtraMolConcentration % kg / m^3
initialSmallParticleMass = solubility*0.01*cdef.molarMass*vesselVolume



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

% get a distribution
rng(1);
meanSize = 41*[1 1.45 1.1 1.3]'*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 = meanSize;
% pbeDefinition.seeding.probabilityVector = disD.probabilityVector;
% pbeDefinition.seeding.numberVector = disD.probabilityVector*0 + ...
%     3e6 / numel(disD.probabilityVector); % #
pbeDefinition.seeding.numberVector = 3e6;

initialPositionMatrix = pbeDefinition.seeding.positionMatrix

pbeDefinition.bulkInitial = initialSoluteMass; % kg
initialState = pbeDefinition.bulkInitial

%% 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 - relative additive roundness
% 11 - unified partition for kernel crystal

% base distribution functional
pbeDefinition.pivotProperty = @(H, up)(...
    [H(1:4, :); ...
     computeMeasureDec(pbeDefinition.dec, 'hC', {up, H(1:3, :)}, {'volume', 'surface', 'meanWidth'}, [], 0); ...
     computeMeasureDec(pbeDefinition.dec, 'hC', {up, H(1:3, :)*0+1}, {'meanWidth'}, [], 0); ...
     2*H(4,:) + computeMeasureDec(pbeDefinition.dec, 'hC', {up, H(1:3, :)}, {'meanWidth'}, [], 0); ...
     2*H(4,:) ./ (...
       2*H(4,:) + ...
       computeMeasureDec(pbeDefinition.dec, 'hC', {up, H(1:3, :)}, {'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(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 = @(x)(x(1)/cdef.molarMass/vesselVolume);

smallParticleMass = @(vCry, x)(initialSmallParticleMass + ...
    (initialSoluteMass - x) + ...
     cdef.density * (initialCrystalVolume - vCry)...
    );

supersaturation = @(x)(soluteConcentration(x)/solubility - 1);
% ..and their initial values
% initialSoluteConcentration = soluteConcentration(initialP, initialN)
initialSupersaturation = supersaturation(initialState)

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

%% 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 state: total solute mass
pbeDefinition.bulkRate = @(p, t)(...
    kSmallParticle * max(smallParticleMass(p(2), p(7)), 0) - ...
    cdef.density*p(6));

%% bulk properties
% 1 - supersaturation
% 2 - total volume
% 3 - NAN -- (was: solute concentration)
% 4 - mean width
% 5 - relative additive roundness
% 6 - volume accumulation by growth (finite difference)
% 7 - bulk state = solute particle mass
% 8 - small pasrticle mass
deltaT = 0.1; %s
pbeDefinition.bulkProperty = @(P, N, x, t)([...
    supersaturation(x); ...       %1
    totalCrystalVolume(P, N); ...    %2
    nan; ...   %3
    P(9, :) * N' / sum(N); ... 
    P(10, :) * N'/ sum(N); ...
    (totalCrystalVolume(pbeDefinition.pivotProperty(P(1:4,:) + ...
          deltaT*rateGrowth(P, supersaturation(x), t), ...
          P(11,:)), N) - ...
     totalCrystalVolume(P, N))/deltaT; ...
    x(1);
    smallParticleMass(totalCrystalVolume(P, N), x(1)) ...
    ]);   

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