
%% rounding fraction computation
qRounding = @(mu)(par.kAttritionMax.*(...
    exp(par.attritionSlope.*mu) - exp(par.attritionSlope)...
    ) ./ (1 - exp(par.attritionSlope)));
% qRounding = @(mu)(0.01)
                       
gAttRounding = @(mu, t)(gAttrition(t) .* qRounding(mu) ./ (1 + qRounding(mu)));

gAttConstShape = @(mu, t)(gAttrition(t) .* 1 ./ (1 + qRounding(mu)));

%% directional derivative
dirDiffMeanWidth = @(H, P)(...
    (P(4,:) - kernelMeanWidthWhileEnforceValid(...
       H - par.finiteDiff, P(1,:)...
       )...
    ) ./ par.finiteDiff);
% dirDiffMeanWidth = @(hB, up)(-1*mwKernel(ones(4, length(up)), up));

%% shape coordinate derivative components
dLamRounding = @(H, P, t, dirDiff)(bsxfun(@times, ...
    gAttRounding(P(10,:), t),  1 ./ (2 - dirDiff)));
dHcRounding = @(H, P, t, dirDiff)(-1*ones(nC, 1)*dLamRounding(H, P, t, dirDiff));

dLamConstShape = @(H, P, t)(min(0, bsxfun(@times, ...
    gAttConstShape(P(10,:), t), 1 ./ (P(4,:) ./ H(nC+1,:) + 2) ...
    )));
% we use min to eliminate nan values
dHcConstShape = @(H, P, t)(min(0, bsxfun(@times, ...
    H(1:3,:), dLamConstShape(H, P, t) ./ H(4,:)) ...
    ));

%% growth components
dHcGrowth = @(H, P, y, t)([...
    physicalRateGrowth(H, P, y, t); ...
    zeros(1, size(H, 2))...
    ]);
dLamGrowth = @(H, P, y, t, oneMw, growthMw)(...
    par.qGrowth ./ (oneMw - 2) .* growthMw ...
    );

rateGrowth = @(H, P, y, t, oneMw, growthMw)(...
    dHcGrowth(H, P, y, t) + ...
    [1; 1; 1; -1]*dLamGrowth(H, P, y, t, oneMw, growthMw) ...
    );

%% assembling the geometrical state derivative
% lam>0 is one constraint that must be considered. Note that attrition is
%   not problematic, here since it shifts to a shape preserving rate down
%   to zero, by definition. The only concerning thing is the derounding via
%   growth. Here, we have to "slow down" the negative lamda rate down to
%   zero BUT we also have to remove the same portions from the hC rate!
% intermediate total rate
rateTotal = @(H, P, y, t, dirDiff, oneMw, growthMw)(...
    [dHcRounding(H, P, t, dirDiff) + dHcConstShape(H, P, t); ...
     dLamRounding(H, P, t, dirDiff) + dLamConstShape(H, P, t)...
    ] + rateGrowth(H, P, y, t, oneMw, growthMw));
% the limiter function
limiter = @(rate, lam)(rate - [-1 -1 -1 1]'*(...
    ((rate(nC+1,:) >= 0)) .* rate(nC+1,:) * 0 + ...
    ((rate(nC+1,:) < 0) & (lam >= par.lamMin)) .* rate(nC+1,:) * 0 + ...
    ((rate(nC+1,:) < 0) & (lam < par.lamMin) & (lam >= 0)) .* ...
      (par.lamMin-lam)./par.lamMin .* rate(nC+1,:) + ...
    ((rate(nC+1,:) < 0) & (lam < 0)) .* rate(nC+1,:)...
    ))
% adding the growth constraint preventing lam<0
pbeDefinition.growthRate = @(H, P, y, t)(...
    limiter(rateTotal(H, P, y, t, ...
        dirDiffMeanWidth(H, P), ...
        oneMeanWidth(H, P(1,:)), ...
        computeMeasureDec(pbeDefinition.dec, 'hC', ...
          {P(1,:), physicalRateGrowth(H, P, y, t)}, 'meanWidth', [], 0) ...
        ), ...
      H(nC+1,:) ...
      ) ...
    );

