function pbeDefinition = growthCommon(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.kAttritionMax = 1;% high implies fast round-off
parStandard.kAttrition = 0e-12;
parStandard.attritionSlope = 1e-4;
parStandard.qGrowth = 1;   % high implies fast grow out of rounding
parStandard.finiteDiff = 2e-8;
parStandard.lamMin = 1e-7;
% time span
parStandard.timeSpan = [0:0.01:(1+10/60)]*60*60;     %s
if ~exist('par', 'var') || isempty(par)
    par = parStandard;
else
    par = mergeStruct(par, parStandard);
end
pbeDefinition.timeSpan = par.timeSpan      %s


%% vessel
temperatureCelciusStart = 50;                         % °C
temperatureCelciusEnd   = 30;                         % °C
temperatureSlopeFinished = 60*60;                     % s
temperatureKelvin = @(t)(max(temperatureCelciusEnd, ...
    temperatureCelciusStart + ...
    (temperatureCelciusEnd-temperatureCelciusStart)*...
    t/temperatureSlopeFinished) + 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 = @(t)(10^5*exp(12.1/temperatureKelvin(t) + 10.47*log(temperatureKelvin(t)) - 65.73)); % mol/m^3
initialSoluteConcentration = solubility(0)*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
nC = pbeDefinition.dec.nC;
% measuree data
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'volume');
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'surface');
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'mean width');

initialSoluteVolume = vesselVolume*initialSoluteConcentration * ...
    cdef.molarMass / cdef.density %m^3
finalSoluteVolume = vesselVolume*solubility(3600) * ...
    cdef.molarMass / cdef.density %m^3

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

% get a distribution
meanSize = [0.01 0.01 0.01 1]'*1e-6;   % m
pbeDefinition.nDecCoordinate = 3;
pbeDefinition.nExtraCoordinate = 1;
pbeDefinition.seeding.time = 0; % seconds
pbeDefinition.seeding.positionMatrix = meanSize * logspace(-1, 1, 100);
pbeDefinition.seeding.numberVector = 2e7 + ...
    pbeDefinition.seeding.positionMatrix(1,:)*0; % #

initialPositiontiMatrix = pbeDefinition.seeding.positionMatrix

%% pivot measures
run partPivotMeasures
% 1 - unified partition for kernel crystal
% 2 - kernel crystal volume
% 3 - kernel crystal surface
% 4 - kernel crystal mean width
% 5 - particle volume
% 6 - particle surface
% 7 - particle mean width
% 8 - muB1
% 9 - muB2
% 10 - muB3
pbeDefinition.pivotProperty = @(H, up)(...
    [up; ...
     pivotMeasures(H, up) ...
     ]);
% ..and their initial values
initialP = pbeDefinition.pivotProperty(...
    pbeDefinition.seeding.positionMatrix, ...
    computePartitionIndexDec(pbeDefinition.dec, ...
      pbeDefinition.seeding.positionMatrix(1:(end-1), :), elkZerotol, 1))
% the number above is assumed to be a DENSITY so that the number estimate
%   for the pivot is:
initialN = pbeDefinition.seeding.numberVector

initialCrystalVolume = initialP(5, :)*initialN'
initialCrystalMass = cdef.density * initialCrystalVolume  % kg
initialCrystalVolumeFraction = initialCrystalVolume / vesselVolume

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

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

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

%% rates: physical base
% crystal growth: k1*s^k2
givenRateGrowth = @(P, super)([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
physicalRateGrowth = @(H, P, y, t)(...
    givenRateGrowth(P, y(1)) ...
    )
% attrition rate (fixed to match experiments well)
gAttrition = @(t)(-1*par.stirrerSpeed * par.kAttrition) % m/s

run partGrowthDynamic

%% bulk properties
% 1 - supersaturation
% 2 - total volume
% 3 - total surface
% 4 - total mean width
% 5 - total number
% 6 - total muB1
% 7 - total muB2
% 8 - total muB3
% 9 - solute concentration
% 10 - "mean" unified partition
% 11 - temperature
pbeDefinition.bulkProperty = @(P, N, x, t)([...
    supersaturation(P, N, t); ...
    P(5, :) * N'; ...
    P(6, :) * N'; ...
    P(7, :) * N'; ...
    sum(N); ...
    P(8, :) * N'; ...
    P(9, :) * N'; ...
    P(10, :) * N'; ...
    soluteConcentration(P, N); ...
    (P(1,:)*N')/sum(N); ...
    temperatureKelvin(t) - 273.15 ...
    ]);   

% test fucntion
initialBulkProperty = pbeDefinition.bulkProperty(initialP, initialN, [], 0)

% some initial value printout
initialGrowthRate = pbeDefinition.growthRate(...
    initialP(1:(nC+1),:), initialP, initialBulkProperty, 0)

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

%% solver options
pbeDefinition.optionStruct.relTol = 1e-5;
pbeDefinition.optionStruct.odeSolver = @ode45;
pbeDefinition.optionStruct.distanceLimit = par.lamMin;


