function pbeDefinition = roundingRatioCommon(par)

%% Name and information
pbeDefinition.name = 'commonPbeDefinition';
pbeDefinition.info = ['potash alum crystal with normalized mean width is ' ...
    'abraded while different parameters for the rounding ratio apply'];
pbeDefinition.source = ['Ma et al. 2008: "Morphological Population ' ...
    'Balance for Modeling Crystal Growth in Face Directions"'];

%% sensible parameters
if ~exist('par', 'var'),par = struct();end
if ~isfield(par, 'kMax')
    par.kMax = 10;
end
if ~isfield(par, 'aSlope')
    par.aSlope = -10;
end

%% initial condition
pbeDefinition.dec = obtainCrystalDec('generic_potash_alum', 0);
cdef = pbeDefinition.dec.cdef;
% par.hC = [1 1.2 1.1]';
par.hC = [1]';
hrep.A = cdef.A;
hrep.h = cdef.groupMappingMatrix * par.hC;
% viewPolytope(hrep)
par.hC = par.hC / computeMeanWidthHrep(hrep);

% decomposition data
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'volume');
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'surface');
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'mean width');

% assign to pbeDefinition
pbeDefinition.nDecCoordinate = 1;
pbeDefinition.nExtraCoordinate = 1;
pbeDefinition.seeding.time = 0;
pbeDefinition.seeding.positionMatrix = [par.hC; 0];
pbeDefinition.seeding.numberVector = 1;


%% governing rate and timing
par.gAttrition = -1;
pbeDefinition.timeSpan = [0:0.01:0.99 0.991:0.001:0.999];

%% some geometric functions based on hB
% 1 - hC1
% 2 - hC2
% 3 - hC3
% 4 - lam

% kernel measures
volKernel = @(hB, up)(computeMeasureDec(pbeDefinition.dec, ...
    'hC', {up, hB(1, :)}, 'volume', [], 0, 1));
surfKernel = @(hB, up)(computeMeasureDec(pbeDefinition.dec, ...
    'hC', {up, hB(1, :)}, 'surface', [], 0, 1));
mwKernel = @(hB, up)(computeMeasureDec(pbeDefinition.dec, ...
    'hC', {up, hB(1, :)}, 'meanWidth', [], 0, 1));
mwKernel2 = @(hB, up)(computeMeasureDec(pbeDefinition.dec, ...
    'hC', hB(1, :), 'meanWidth', [], 0, 1));

% total measures
volTotal = @(hB, up)(volKernel(hB, up) + ...
                     surfKernel(hB, up).*hB(2, :) + ...
                     2.*pi.*mwKernel(hB,up).*hB(2, :).^2 + ...
                     4.*pi/3.*hB(2,:).^3);
surfTotal = @(hB, up)(surfKernel(hB, up) + ...
                      4.*pi.*mwKernel(hB,up).*hB(2,:) + ...
                      4.*pi.*hB(2,:).^2);
                  
mwTotal = @(hB, up)(mwKernel(hB,up) + 2.*hB(2,:));

% roundness values
% .. the max eliminates any nan that 0/0 produces
muB1 = @(hB, up)(min(1, 1 - mwKernel(hB, up) ./ mwTotal(hB, up)));
muB2 = @(hB, up)(min(1, 1 - surfKernel(hB, up) ./ surfTotal(hB, up)));
muB3 = @(hB, up)(min(1, 1 - volKernel(hB, up) ./ volTotal(hB, up)));

% muB1([0 0 0 0]', 1)
% muB2([0 0 0 0]', 1)
% muB3([0 0 0 0]', 1)
% 
% volTotal([0 0 0 0]', 1)
% surfTotal([0 0 0 0]', 1)
% mwTotal([0 0 0 0]', 1)
% 
% volKernel([0 0 0 0]', 1)
% surfKernel([0 0 0 0]', 1)
% mwKernel([0 0 0 0]', 1)


%% dynamic functions

% rounding fraction
% qRounding = @(mu)(par.kMax.*(exp(par.aSlope.*mu) - exp(par.aSlope)) ./ ...
%                            (1 - exp(par.aSlope)));
qRounding = @(mu)((mu < 0.9)*par.kMax + 0.001)
                       
gRounding = @(mu)(par.gAttrition .* qRounding(mu) ./ (1 + qRounding(mu)));
gConstShape = @(mu)(par.gAttrition .* 1 ./ (1 + qRounding(mu)));

par.finiteDiff = 2e-8;
dirDiffMeanWidth = @(hB, up)(...
    (mwKernel2(hB - par.finiteDiff, up) - mwKernel(hB, up)) ./ ...
     par.finiteDiff...
     );
% dirDiffMeanWidth = @(hB, up)(-1*mwKernel(ones(2, length(up)), up));

dLamRounding = @(hB, up)(bsxfun(@times, ...
    gRounding(muB3(hB, up)), ...
    1 ./ (2 + dirDiffMeanWidth(hB, up))));
dHcRounding = @(hB, up)(-1*dLamRounding(hB, up));

dLamConstShape = @(hB, up)(min(0, bsxfun(@times, ...
    gConstShape(muB3(hB, up)), ...
    1 ./ (mwKernel(max(0, bsxfun(@times, hB(1,:), 1./hB(2,:))), up) + 2))));
% again, we use min to eliminate nan values
dHcConstShape = @(hB, up)(min(0, bsxfun(@times, ...
    bsxfun(@times, hB(1,:), 1 ./ hB(2,:)), ...
    dLamConstShape(hB, up))));

pbeDefinition.growthRate = @(H, P, s, t)([...
    dHcRounding(P(1:2,:), P(3,:)) + dHcConstShape(P(1:2,:), P(3,:)); ...
    dLamRounding(P(1:2,:), P(3,:)) + dLamConstShape(P(1:2,:), P(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 - up
%  6 - total mean widht
%  7 - total surface area
%  8 - total volume
%  9 - muB1
% 10 - muB2
% 11 - muB3

% base distribution functional
pbeDefinition.pivotProperty = @(H, up)(...
    [H(1:2, :); ...
     up; ...
     mwTotal(H, up); ...
     surfTotal(H, up); ...
     volTotal(H, up); ...
     muB1(H, up); ...
     muB2(H, up); ...
     muB3(H, up)]);
 
%% testing
% hB0 = [0 0]';
% hBI = pbeDefinition.seeding.positionMatrix;
% hB1 = [1 1]';
% P0 = pbeDefinition.pivotProperty(hB0, 1)
% PI = pbeDefinition.pivotProperty(hBI, 1)
% P1 = pbeDefinition.pivotProperty(hB1, 1)
% 
% g0 = pbeDefinition.growthRate([], P0, [], [])
% gI = pbeDefinition.growthRate([], PI, [], [])
% g1 = pbeDefinition.growthRate([], P1, [], [])
% 
% dLamRounding(hB0, 1)
% dHcRounding(hB0, 1)
% dLamConstShape(hB0, 1)
% dHcConstShape(hB0, 1)
% [volKernel(hB0, 1) volTotal(hB0, 1) muB3(hB0, 1) ...
%     gRounding(muB3(hB0, 1)) gConstShape(muB3(hB0, 1))]
% 
% 
% dLamRounding(hBI, 1)
% dHcRounding(hBI, 1)
% dLamConstShape(hBI, 1)
% dHcConstShape(hBI, 1)
% [volKernel(hBI, 1) volTotal(hBI, 1) muB3(hBI, 1) ...
%     gRounding(muB3(hBI, 1)) gConstShape(muB3(hBI, 1))]
% 
% dLamRounding(hB1, 1)
% dHcRounding(hB1, 1)
% dLamConstShape(hB1, 1)
% dHcConstShape(hB1, 1)
% [volKernel(hB1, 1) volTotal(hB1, 1) muB3(hB1, 1) ...
%     gRounding(muB3(hB1, 1)) gConstShape(muB3(hB1, 1))]
% 
% hB = [1 1]';
% d = 1e-2;
% [mwKernel(hB, 1) mwKernel(hB-1e-2, 1) ...
%  (mwKernel(hB-d, 1) - mwKernel(hB, 1)) / d ...
%  dirDiffMeanWidth(hB, 1)]
% hB = [1 9   9   1]';
% [mwKernel(hB, 1) mwKernel(hB-1e-2, 1) ...
%  (mwKernel(hB-d, 1) - mwKernel(hB, 1)) / d ...
%  dirDiffMeanWidth(hB, 1)]
% hB = [9 1.2 9   1]';
% [mwKernel(hB, 1) mwKernel(hB-1e-2, 1) ...
%  (mwKernel(hB-d, 1) - mwKernel(hB, 1)) / d ...
%  dirDiffMeanWidth(hB, 1)]
% hB = [9 9   1.1 1]';
% [mwKernel(hB, 1) mwKernel(hB-1e-2, 1) ...
%  (mwKernel(hB-d, 1) - mwKernel(hB, 1)) / d ...
%  dirDiffMeanWidth(hB, 1)]


%% bulk properties
% 1 - mw total
% 2 - surf total
% 3 - vol total
% 4 - muB1
% 5 - muB2
% 6 - muB3
pbeDefinition.bulkProperty = @(P, N, x, t)([...
    sum(P(4:9,:), 2) ...
    ]);   

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

pbeDefinition.optionStruct.relTol = 1e-5;
pbeDefinition.optionStruct.odeSolver = @ode113;
pbeDefinition.optionStruct.distanceLimit = 1e-14;