function momentData = simulateParacetamolMoments(pbeDef, dataApprox)

dec = dataApprox.dataApprox.dec;
dec = addMonomialDataDec(dec, 3);
decExt = dataApprox.dataApprox.dec.properData;

%% convert initial condition
hcMatrix = pbeDef.seeding.positionMatrix;
numberVector = pbeDef.seeding.numberVector;
heMatrix = dataApprox.meanMappingMatrix * hcMatrix;
haMatrix = dataApprox.dataApprox.dec.properData.mappingProperToNew * heMatrix;
monoVector = convertHcMatrixToMonomials(haMatrix, numberVector, dec);

%% measure handles
% volume
totalVolume = @(monoVec)(dec.measureData.volume.tensorHc{1}(...
    dec.monomialTensorIndexMatrix(...
      dec.monomialTensorIndexMatrix(:, 3) ~= 0, 3))' * ...
    monoVec(dec.monomialTensorIndexMatrix(:, 3) ~= 0, :));
refVolume = numberVector * computeMeasureDec(decExt.originalData, 'hc', ...
    hcMatrix, 'volume')';
% [totalVolume(monoVector) refVolume]

% surface
totalSurface = @(monoVec)(dec.measureData.surface.tensorHc{1}(...
    dec.monomialTensorIndexMatrix(...
      dec.monomialTensorIndexMatrix(:, 2) ~= 0, 2))' * ...
    monoVec(dec.monomialTensorIndexMatrix(:, 2) ~= 0, :));
refSurface = numberVector * computeMeasureDec(...
    decExt.originalData, 'hc', pbeDef.seeding.positionMatrix, 'surface')';
% [totalSurface(monoVector) refSurface]

% mean width
totalWidth = @(monoVec)(dec.measureData.meanWidth.tensorHc{1}(...
    dec.monomialTensorIndexMatrix(...
      dec.monomialTensorIndexMatrix(:, 1) ~= 0, 1)) * ...
    monoVec(dec.monomialTensorIndexMatrix(:, 1) ~= 0, :));
refSurface = numberVector * computeMeasureDec(...
    decExt.originalData, 'hc', pbeDef.seeding.positionMatrix, 'meanWidth')';
% [totalWidth(monoVector) refSurface]

% support 40
supportMatrix = nan(decExt.nH, dec.nC);
for iFacet = 1:decExt.nH
    supportMatrix(iFacet, :) = dec.measureData.support.tensorHc{iFacet}{1};
end
support = @(monoVec)(supportMatrix(:, ...
    dec.monomialTensorIndexMatrix(...
      dec.monomialTensorIndexMatrix(:, 1) ~= 0, 1)) * ...
    monoVec(dec.monomialTensorIndexMatrix(:, 1) ~= 0, :) / ...
    monoVec(end));
refSupport = mean(decExt.mappingReducedToFull*...
    dataApprox.meanMappingMatrix*hcMatrix, 2);
% [support(monoVector) refSupport]'

% mean shape
meanShape = @(monoVec)(...
    dec.properData.mappingNewToProper * ...
    monoVec((end-1):-1:(end-dec.nC), :) / monoVec(end));
meanShapeRef = mean(heMatrix, 2);
% [meanShape(monoVector) meanShapeRef]

%% extract growth rate
% this part is muc more special to the paracetamol case. First, we note
%   from its source in database/paracetamolGrowth that the growth rate is 
%   computed by @(H, P, y, t) where y(4) is required while H is only a
%   dummy and P or t are irrelevant.
% The parameter y represents the bulk properties which are computed by some
%   @(P, N, x, t). Here, t is used for the temperature profile, x is unused
%   and N*P(:,2) is used for the total volume. Hence, this function can
%   totally be tricked:
bulkProperties = @(t, monoVec)(pbeDef.bulkProperty([0; totalVolume(monoVec)], 1, ...
    [], t));
% and then the growth rate is given by
growth = @(t, monoVec)([dec.properData.mappingProperToNew*...
    dataApprox.meanMappingMatrix*...
    (pbeDef.growthRate(1, [], bulkProperties(t, monoVec), []) + 0e-8); ...
    0] );
monomialDerivativeIndexMatrix = dec.monomialDerivativeIndexMatrix;
monomialDerivativeIndexMatrix(monomialDerivativeIndexMatrix==0) = ...
    dec.nC+1;
monoDeriv = @(growth, monoVec)((dec.monomialDerivativeFactorMatrix.*...
    growth(monomialDerivativeIndexMatrix)) * monoVec);
rhs = @(t, monoVec)(monoDeriv(growth(t, monoVec), monoVec));

[momentData.timeVector, momentData.monoMatrix] = ...
    ode45(rhs, pbeDef.timeSpan, monoVector, odeset(...
    'abstol', 1e-20, 'reltol', 1e-6));

momentData.timeVector = momentData.timeVector';
momentData.monoMatrix = momentData.monoMatrix';
momentData.volumeVector = totalVolume(momentData.monoMatrix);
momentData.surfaceVector = totalSurface(momentData.monoMatrix);
momentData.widthVector = totalWidth(momentData.monoMatrix);
momentData.heMeanMatrix = meanShape(momentData.monoMatrix);
momentData.supportMatrix = support(momentData.monoMatrix);
momentData.bulkMatrix = nan(4, length(momentData.timeVector));
for iTime = 1:length(momentData.timeVector)
    momentData.bulkMatrix(:, iTime) = bulkProperties(...
        momentData.timeVector(iTime), ...
        momentData.monoMatrix(:, iTime));
end

save('data/paracetamolGrowth_4u.mat', 'momentData');