function pbeDef = linearGrowth_common

%% Name and information
pbeDef.info = ['linear growth for a cuboid (nC=3), ' ...
    'including evaluation of distribution densities and analytical ' ...
    'comparison.'];
pbeDef.source = 'generic';

% In this example, we choose a simplified 2D case for potash alum. In doing
%   so, we should be able to print some density distribuitions, which would
%   be impussible, otherwise.
pbeDef.dec = obtainCrystalDec('cuboid');
pbeDef.dec = addMeasureDataDec(pbeDef.dec, 'volume');
pbeDef.dec = addMeasureDataDec(pbeDef.dec, 'surface');
pbeDef.dec = addMeasureDataDec(pbeDef.dec, 'meanWidth');

%% initial condition
% We use an exponential distribution that is truncated by the validity
%   cone.
pbeDef.seeding.nSample = 100;
pbeDef.seeding.meanSize = 1e-5;
pbeDef.seeding.totalNumber = 1e10;
% the resulting number density is a multivariate exponential distribution;
% note that it must be truncated for hC < 0
pbeDef.seeding.fhandle = @(hC)(...
    pbeDef.seeding.totalNumber * ...
    (1/pbeDef.seeding.meanSize)^3 * ...
    prod(exp(-1/pbeDef.seeding.meanSize * hC), 1) .* ...
    all(hC >= 0, 1));

% moments
pbeDef.seeding.firstMoment = @(h)(...
    bsxfun(@times, pbeDef.seeding.fhandle(h), h));

pbeDef.seeding.thirdMoment = @(h)(...
    pbeDef.seeding.fhandle(h) .* prod(h, 1));

% measures
pbeDef.seeding.volumeHandle = @(h)(...
    pbeDef.seeding.fhandle(h) .* ...
    computeMeasureDec(pbeDef.dec, 'hC', h, 'volume', 1e-16, [], -1));

pbeDef.seeding.surfaceHandle = @(h)(...
    pbeDef.seeding.fhandle(h) .* ...
    computeMeasureDec(pbeDef.dec, 'hC', h, 'surface', 1e-16, [], -1)); 

pbeDef.seeding.meanWidthHandle = @(h)(...
    pbeDef.seeding.fhandle(h) .* ...
    computeMeasureDec(pbeDef.dec, 'hC', h, 'meanWidth', 1e-16, [], -1));

%% maxima
% hInit = 1./pbeDefinition.seeding.spreadConstant'
% numberCenter = fminsearch(...
%     @(x)(-1*pbeDefinition.seeding.fhandle(x)), hInit)
% volumeCenter = fminsearch(...
%     @(x)(-1*pbeDefinition.seeding.volumeHandle(x)), hInit)
% surfaceCenter = fminsearch(...
%     @(x)(-1*pbeDefinition.seeding.surfaceHandle(x)), hInit)
% meanWidthCenter = fminsearch(...
%     @(x)(-1*pbeDefinition.seeding.meanWidthHandle(x)), hInit)

%% reference integrals
% integration bounds for the disC
intMin = 0*[1 1 1]';
intMax = 5*pbeDef.seeding.meanSize*[1 1 1]';

% commonDivonneOptions = {'algorithm', 'Divonne', ...
%     'epsrel', 1e-2, 'maxeval', 2e5, ...
%     'key1', 7, 'key2', -1, 'key3', 0, ...
%     'flags', 16+3, ...
%     'xgiven', 0.5*pbeDefinition.seeding.spreadConstant'};
commonVegasOptions = {'algorithm', 'vegas', ...
    'epsrel', 1e-3, 'maxeval', 1e5, ...
    'flags', 16+1};
% pbeDefinition.vegasOptions = commonVegasOptions;
% pbeDefinition.integrateMin = intMin;
% pbeDefinition.integrateMax = intMax;

thisFileName = 'linearGrowthIntegrals_0.matx';
if ~exist(thisFileName, 'file')
    % template of integral information
    integralTemplate = struct('number', [], ...
        'volume', [], 'surface', [], 'meanWidth', [], ...
        'firstMoment', [], 'secondMoment', [], 'thirdMoment', []);
    integral = integralTemplate;
        
    % number
    integral.number = integrateCuba(pbeDef.seeding.fhandle, ...
        intMin, intMax, commonVegasOptions{:});
    intNumber = integral.number
    
    % moments
    integral.firstMoment = integrateCuba(pbeDef.seeding.firstMoment, ...
        intMin, intMax, commonVegasOptions{:});
    intFirstMoment = integral.firstMoment
    
    integral.thirdMoment = integrateCuba(pbeDef.seeding.thirdMoment, ...
        intMin, intMax, commonVegasOptions{:});
    intThirdMoment = integral.thirdMoment
    
    % measures
    integral.volume = integrateCuba(pbeDef.seeding.volumeHandle, ...
        intMin, intMax, commonVegasOptions{:});
    intVolume = integral.volume
    
    integral.surface = integrateCuba(pbeDef.seeding.surfaceHandle, ...
        intMin, intMax, commonVegasOptions{:});
    intSurface = integral.surface
    
    integral.meanWidth = integrateCuba(pbeDef.seeding.meanWidthHandle, ...
        intMin, intMax, commonVegasOptions{:});
    intMeanWidth = integral.meanWidth
                  
    save(thisFileName, '-mat', 'integral', 'integralTemplate')
else
    load(thisFileName, '-mat');
    if false
        intNumber = integral.number
        intFirstMoment = integral.firstMoment
        intFirstMomentValue = integral.firstMoment.value'
        intFirstMomentError = integral.firstMoment.error'
        intFirstMomentProb = integral.firstMoment.prob'
        intThirdMoment = integral.thirdMoment
        intVolume = integral.volume
        intSurface = integral.surface
        intMeanWidth = integral.meanWidth
    end
end

%% growth and analytical solution
pbeDef.growthConstant = 1;
pbeDef.growthRate = @(H, P, y, t)(pbeDef.growthConstant*H);
pbeDef.analyticSolution = @(h, t)(...
    pbeDef.seeding.fhandle(h * exp(-t*pbeDef.growthConstant)) * ...
    exp(-pbeDef.dec.nC*t));

pbeDef.analyticNumber = @(t)(pbeDef.seeding.totalNumber * (t*0+1));
pbeDef.analyticVolume = @(t)(pbeDef.seeding.totalNumber * ...
    (2*pbeDef.seeding.meanSize * ...
     exp(pbeDef.growthConstant*t)).^3);

pbeDef.analyticSurface = @(t)(pbeDef.seeding.totalNumber * 6 * ...
     (2*pbeDef.seeding.meanSize * ...
     exp(pbeDef.growthConstant*t)).^2);

 pbeDef.analyticMeanWidth = @(t)(pbeDef.seeding.totalNumber * 3 * ...
     (pbeDef.seeding.meanSize * ...
     exp(pbeDef.growthConstant*t)));

%% time span for solution
% We simulate growth for half an hour, giving 17 micrometers.
% pbeDefinition.timeSpan = 0:60:(0.5*60*60); %s
pbeDef.timeSpan = 0:0.2:3; %s

% options
% pbeDefinition.optionStruct.distanceLimit = 3e-7; %m
pbeDef.optionStruct.density = 1;
% pbeDefinition.optionStruct.gradientAbsoluteDelta = 1e-10;
pbeDef.optionStruct.gradientRelativeDelta = 1e-3;
pbeDef.optionStruct.relTol = 1e-7;
pbeDef.optionStruct.absTolDensity = 1e-30;
pbeDef.optionStruct.absTolNumber = 1e-30;
% pbeDefinition.optionStruct.odeSolver = @ode113; % 113, 45