function pbeDefinition = seededGrowthUsingDensity

%% Name and information
pbeDefinition.info = ['seeded growth (one initial population), ' ...
    'using a growth field that accelerates, decellerates ' ...
    'and goes against a boundary'];
pbeDefinition.source = 'generic';

%% Crystal setup
% crystal
pbeDefinition.cdef = obtainCrystal('generic_potash_alum', 1);

% decomposition data
pbeDefinition.dec = obtainDec(pbeDefinition.cdef.A, 'mapping', ...
    pbeDefinition.cdef.groupMappingMatrix, 'verbose', 0);
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'volume');
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'surface');
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'meanWidth');



%% Initial distribution
rng(1);
pbeDefinition.seeding.time = 0; % seconds
std = 0.5*[3 3]'*1e-6;   % m
disD = obtainDistributionDiscrete('normal', ...
    'nSample', 100, ...
    'mean', [10 10]'*1e-6, ...  %m
    'cov', diag(std.^2));
pbeDefinition.seeding.positionMatrix = disD.positionMatrix;       % m
pbeDefinition.seeding.numberVector = disD.probabilityVector*0 + 1e10; % #
pbeDefinition.seeding.probabilityVector = disD.probabilityVector;
pbeDefinition.seeding.probabilityVector

%% growth and properties
% shortcut
dec = pbeDefinition.dec;

% growth rate
% crystals start at mean [1 1]*1e-5 and grow with [1 1]*5e-9.
% For hC1 > 2e-5 the rate in hC1 direction doubles until 3e-5.
% For hC2 >
% 2.5e-9, the rate halfs down until 4e-5.
factorOne = @(H)(1.0 + ...
    (H(1,:) > 2e-5).*(H(1,:) < 2.5e-5).*((H(1,:)-2e-5)/0.5e-5)*4.0 + ...
    (H(1,:) >= 2.5e-5)*4.0...
    );
% factorOne(pbeDefinition.seeding.positionMatrix)
factorTwo = @(H)(1.0 - ...
    (H(2,:) > 1.5e-5).*(H(2,:) < 1.8e-5).*((H(2,:)-1.5e-5)/0.3e-5)*0.75 - ...
    (H(2,:) >= 1.8e-5)*0.75...
    );
% factorOne(pbeDefinition.seeding.positionMatrix)

pbeDefinition.growthRate = @(H, P, y, t)(...
    [5; 10]*1e-9 * ones(1, size(H, 2)) .* ...
    [factorOne(H); ones(1, size(H, 2))] .* ...
    [ones(1, size(H, 2)); factorTwo(H)] );

% properties
pbeDefinition.pivotProperty = {'volume'};
pbeDefinition.bulkProperty = @(P, N, x, t)(P*N');

% time span for solution
pbeDefinition.timeSpan = 0:30:3.5e3;          %s

% options
pbeDefinition.optionStruct.distanceLimit = 3e-7; %m
pbeDefinition.optionStruct.density = 1;