function pbeDefinition = constGrowth_common

%% Name and information
pbeDefinition.info = ['seeded growth 7-dimensional, ' ...
    'with constant growth rate for which an analytical solution' ...
    'is provided'];
pbeDefinition.source = 'Kumar and Ramkrischna 1997 WITH generalization in technical report';

% We choose a rather high-dimensional case, basically to show off and
%   demonstrate that the dimension is NOT an issue using the characteristic
%   and MC-integration approach
% The decomposition is only computed once to save computation time. The
%   measure generation takes quite a while.
if exist('constGrowthDec.matx', 'file')
    load('constGrowthDec.matx', '-mat')
else
    pbeDefinition.dec = obtainCrystalDec('generic_paracetamol', 15, 0);
%     pbeDefinition.dec = obtainCrystalDec('generic_potash_alum', 3, 1);
    pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'volume');
    pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'surface');
    pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'meanWidth');
    pbeDefinition.dec
    save('constGrowthDec.matx', '-mat')
end
nDim = pbeDefinition.dec.nC;

%% initial condition
% We use a normal distribution that is truncated by the validity. This
%   introduces a small error when the distribution is sampled by the
%   corresponding normal distribution. This error is mainly controlled by
%   the standard deviation.
% distribution parameters
pbeDefinition.seeding.mean = ones(nDim,1)*1e-6;
pbeDefinition.seeding.std = ones(nDim,1)*1.5e-7;
pbeDefinition.seeding.cov = diag(pbeDefinition.seeding.std.^2);
pbeDefinition.seeding.margin = 2*pbeDefinition.seeding.std;
pbeDefinition.seeding.totalValue = 1e10;
pbeDefinition.seeding.nSample = 100;
disC = obtainDistributionContinuous('normal', ...
    'mean', pbeDefinition.seeding.mean, ...
    'cov', pbeDefinition.seeding.cov, ...
    'totalValue', pbeDefinition.seeding.totalValue);

% helpers to apply truncation by margin and validity cone
pbeDefinition.seeding.marginFilter = @(h)(double(all(...
    abs(bsxfun(@minus, h, pbeDefinition.seeding.mean)) <= ...
        pbeDefinition.seeding.margin(1))));

pbeDefinition.seeding.validityFilter = @(h)(...
    isValidDec(h, pbeDefinition.dec, 1e-16));

% handles for function evaluations (note that for measures, the
%   distribution density is automatically truncated by validity through the
%   measures)
pbeDefinition.seeding.normalHandle = disC.fhandle;
pbeDefinition.seeding.marginHandle = @(h)(disC.fhandle(h) .* ...
    pbeDefinition.seeding.marginFilter(h));
pbeDefinition.seeding.validityHandle = @(h)(disC.fhandle(h) .* ...
    pbeDefinition.seeding.validityFilter(h));

pbeDefinition.seeding.volumeHandle = @(h)(disC.fhandle(h) .* ...
    pbeDefinition.seeding.marginFilter(h) .* ...
    computeMeasureDec(pbeDefinition.dec, 'hC', h, 'volume', 1e-16, [], -1));

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

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

pbeDefinition.seeding.moment3Handle = @(h)(disC.fhandle(h) .* ...
    pbeDefinition.seeding.marginFilter(h) .* ...
    pbeDefinition.seeding.validityFilter(h) .* ...
    sum(h.^3, 1));
pbeDefinition.seeding.moment2Handle = @(h)(disC.fhandle(h) .* ...
    pbeDefinition.seeding.marginFilter(h) .* ...
    pbeDefinition.seeding.validityFilter(h) .* ...
    sum(h.^2, 1)); 
pbeDefinition.seeding.moment1Handle = @(h)(disC.fhandle(h) .* ...
    pbeDefinition.seeding.marginFilter(h) .* ...
    pbeDefinition.seeding.validityFilter(h) .* ...
    sum(h, 1)); 

% final continuous particle distribution density (margin AND validity)
pbeDefinition.seeding.fhandle = @(h)(disC.fhandle(h) .* ...
    pbeDefinition.seeding.marginFilter(h) .* ...
    isValidDec(h, pbeDefinition.dec, 1e-16));

% integration bounds for the disC
intMin = max(pbeDefinition.seeding.mean - pbeDefinition.seeding.margin, 0);
intMax = pbeDefinition.seeding.mean + pbeDefinition.seeding.margin;
intMinExtra = max(pbeDefinition.seeding.mean - 2.5*pbeDefinition.seeding.std, 0);
intMaxExtra = pbeDefinition.seeding.mean + 2.5*pbeDefinition.seeding.std;
  
%% evaluate integrals for the IC
% commonDivonneOptions = {'algorithm', 'Divonne', ...
%     'epsrel', 1e-2, 'maxeval', 2e5, ...
%     'key1', 7, 'key2', -1, 'key3', 0, ...
%     'flags', 16+3, ...
%     'xgiven', pbeDefinition.seeding.mean};
thisVegasOptions = {'algorithm', 'vegas', ...
    'epsrel', 1e-4, 'maxeval', 5e8, ...
    'nBatch', 5e4, 'flags', 16+1};
pbeDefinition.vegasOptions =  {'algorithm', 'vegas', ...
    'epsrel', 1e-3, 'maxeval', 5e8, ...
    'nBatch', 5e4, 'flags', 16+1};
pbeDefinition.integrateMin = intMin;
pbeDefinition.integrateMax = intMax;

thisFileName = 'constGrowthIntegrals_0.matx';
if ~exist(thisFileName, 'file')
    % template of integral information
    integralTemplate = struct('normal', [], 'normalMarginApplied', [], ...
        'normalValidityApplied', [], 'number', [], ...
        'volume', [], 'surface', [], 'meanWidth', [], ...
        'moment1', [], 'moment2', [], 'moment3', []);
    integral = integralTemplate;
    
    % number integrals
    integral.normal = integrateCuba(pbeDefinition.seeding.normalHandle, ...
        intMinExtra, intMaxExtra, thisVegasOptions{:});
    intNormal = integral.normal
    
    integral.normalMarginApplied = integrateCuba(...
        pbeDefinition.seeding.marginHandle, ...
        intMin, intMax, thisVegasOptions{:})
    intMargin = integral.normalMarginApplied
    
    integral.normalValidityApplied = integrateCuba(...
        pbeDefinition.seeding.validityHandle, ...
        intMin, intMax, thisVegasOptions{:});
    intValid = integral.normalValidityApplied
    
    % important integrals
    integral.number = integrateCuba(pbeDefinition.seeding.fhandle, ...
        intMin, intMax, thisVegasOptions{:});
    intNumber = integral.number
    
    integral.volume = integrateCuba(pbeDefinition.seeding.volumeHandle, ...
        intMin, intMax, thisVegasOptions{:});
    intVolume = integral.volume
    
    integral.surface = integrateCuba(pbeDefinition.seeding.surfaceHandle, ...
        intMin, intMax, thisVegasOptions{:});
    intSurface = integral.surface
    
    integral.meanWidth = integrateCuba(pbeDefinition.seeding.meanWidthHandle, ...
        intMin, intMax, thisVegasOptions{:});
    intMeanWidth = integral.meanWidth
    
    % one moment integral
    integral.moment3 = integrateCuba(pbeDefinition.seeding.moment3Handle, ...
        intMin, intMax, thisVegasOptions{:});
    intMoment3 = integral.moment3
    
    integral.moment2 = integrateCuba(pbeDefinition.seeding.moment2Handle, ...
        intMin, intMax, thisVegasOptions{:});
    intMoment2 = integral.moment2
    
    integral.moment1 = integrateCuba(pbeDefinition.seeding.moment1Handle, ...
        intMin, intMax, thisVegasOptions{:});
    intMoment1 = integral.moment1
    
    marginFraction = intMargin.value / intNormal.value
    validFraction = intValid.value / intNormal.value
    numberFraction = intNumber.value / intNormal.value
              
    save(thisFileName, '-mat', 'integral', 'integralTemplate')
else
    load(thisFileName, '-mat');
    if false
        intNormal = integral.normal
        intMargin = integral.normalMarginApplied
        intValid = integral.normalValidityApplied
        intNumber = integral.number
        intVolume = integral.volume
        intSurface = integral.surface
        intMeanWidth = integral.meanWidth
        
        intMoment3 = integral.moment3
        intMoment2 = integral.moment2
        intMoment1 = integral.moment1
        
        marginFraction = intMargin.value / intNormal.value
        validFraction = intValid.value / intNormal.value
        numberFraction = intNumber.value / intNormal.value
    end
end

%% temp-test
% tic
% % commonCuhreOptions = {'algorithm', 'cuhre', ...
% %     'epsrel', 1e-3, 'maxeval', 1e5, 'key', 9, ...
% %     'flags', 16+3};
% commonSuaveOptions = {'algorithm', 'suave', ...
%     'epsrel', 1e-2, 'maxeval', 1e5, ...
%     'flatness', 2, 'nnew', 1000, ...
%     'flags', 16+3};
% commonVegasOptions = {'algorithm', 'vegas', ...
%     'epsrel', 1e-3, 'maxeval', 1e6, ...
%     'flags', 16+1};
% % commonVegasOptions{8} = 2e5;
% intVolume2 = integrateCuba(pbeDefinition.seeding.fhandle, ...
%     intMin, intMax, commonVegasOptions{:})
% toc
% intVolume

%% growth and analytical solution
% we choose the growth constant to be 0.01 micrometers per second in all
%   directions
pbeDefinition.growthConstant = 1e-8; 
pbeDefinition.growthRate = @(H, P, y, t)(H*0 + pbeDefinition.growthConstant);
pbeDefinition.analyticSolution = @(h, t)(...
    pbeDefinition.seeding.fhandle(h - t*pbeDefinition.growthConstant));

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

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