function pbeDefinition = snyderCuboidGrowthDissolution
clear all
%% Basic information
pbeDefinition.info = ['seeded growth AND dissolution of cuboids with ' ...
    'new facets evolving on edges (stepped faces)'];
tSwitch = 0.06;
fracSwitch = 0.425;
tEnd = 30*tSwitch;

% obtain decomposition
dec = obtainSynderCuboidDec2;
pbeDefinition.source = dec.cdef.source;
pbeDefinition.dec = dec;

%% initial condition
% % NORMAL distribution for cuboid facets. The other facets are exectly
% %   disappeared, giving some dirac functions that would be multiplied.
% pbeDefinition.seeding.mean = [1 2 3]*500*1e-6;
% pbeDefinition.seeding.std = [0.5 1 2]'*50*1e-6;
% % pbeDefinition.seeding.std = [0.5 1 2]'*1e-6;
% pbeDefinition.seeding.cov = diag(pbeDefinition.seeding.std.^2);
% pbeDefinition.seeding.totalValue = 1e10;
% pbeDefinition.seeding.nSample = 1;
% 
% % we sample the cuboid coordinates directly..
% rng(0)
% disTmp = obtainDistributionDiscrete('normal', ...
%     'nSample',pbeDefinition.seeding.nSample, ...
%     'mean', pbeDefinition.seeding.mean, ...
%     'cov', pbeDefinition.seeding.cov);
% disD.positionMatrix = disTmp.positionMatrix;
% % disD.positionMatrix = disTmp.positionMatrix(:, 20:end);
% 
% % check standard deviation
% % checkStd = std(disD.positionMatrix, [], 2)'
% % ..construct the mathing additional coordinates..
% additionalCoordinateMatrix = [zeros(4,3) eye(4)]*...
%     mapToValidityCone([disD.positionMatrix; ...
%                        100*ones(4,1)*max(disD.positionMatrix,[], 1)], ...
%                        dec, 1e-12);
% % ..add a small numerical correction (validity mapping is not perfect)
% additionalCoordinateMatrix = bsxfun(@minus, additionalCoordinateMatrix, ...
%     [1 1 1 2]'*1e-15);
% % ..and finally add the information for the initial distribution:
% pbeDefinition.seeding.positionMatrix = [disD.positionMatrix; ...
%     additionalCoordinateMatrix];
% pbeDefinition.seeding.numberVector = disD.positionMatrix(1,:)*0 + ....
%     pbeDefinition.seeding.totalValue/pbeDefinition.seeding.nSample;

% for initialization, we use a single particle
pbeDefinition.seeding.positionMatrix = dec.cdef.relativeGrowthOne;
% pbeDefinition.seeding.positionMatrix = [1 1 1 10 10 10];
pbeDefinition.seeding.numberVector = 1;
% ..correct validity NOW .. .. this would otherwise happen in the pbe
% solver
pbeDefinition.seeding.positionMatrix = mapToConfinementCone(...
    pbeDefinition.seeding.positionMatrix, pbeDefinition.dec, 0);

%% growth rates
% the relative growth rates are already specified in the crystal definition
%   where two versions are available.
% pbeDefinition.growthRate = @(H, P, y, t)(...
%     (1 - 2*(mod(t, tSwitch) < fracSwitch*tSwitch))*...
%     dec.cdef.relativeGrowthOne*...
%     ones(1, size(H,2))...
%     );
pbeDefinition.growthRate = @(H, P, y, t)((...
    (mod(t, tSwitch) >= fracSwitch*tSwitch)*dec.cdef.relativeGrowthTwo - ...
    (mod(t, tSwitch) < fracSwitch*tSwitch)*dec.cdef.relativeDissolutionTwo)*...
    ones(1, size(H,2)));

%% pivot properties
% we need the unified partitions in P for growth rate projections
pbeDefinition.pivotProperty = @(H, up)(up);

%% time span for solution
pbeDefinition.timeSpan = 0:0.05*tSwitch:tEnd;

% options
pbeDefinition.optionStruct.distanceLimit = 1e-4; %m
pbeDefinition.optionStruct.relTol = 1e-6;
pbeDefinition.optionStruct.absTolNumber = 1e-30;
pbeDefinition.optionStruct.absTolPosition = 1e-30;
% pbeDefinition.optionStruct.updateLimit = 2;
% pbeDefinition.optionStruct.updateTime = 1;

% pbeDefinition.optionStruct.odeSolver = @ode113; % 113, 45