function pbeDefinition = growthHighDimCoupled

%% Name and information
pbeDefinition.info = ['seeded growth (one initial population), ' ...
    'using an arbitrary high faceted crystel ang growth rates, ' ...
    'coupled with mass balance (supersaturation)'];
pbeDefinition.source = 'generic';

%% Crystal setup

% decomposition data
% pbeDefinition.dec = obtainCrystalDec('generic_paracetamol', 15);
% pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'volume');
% pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'surface');
% pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'meanWidth');
% save tmpDec
load tmpDec

%% Initial distribution
% the particle distribution
distributionMean = [1 1 1 1 1 1 1]' * 1e-6;
distributionStd = [0.1 1 1.3 2.3 3.14 4.2 7]' * 1e-8;
particleDisC = obtainDistributionContinuous('normal', ...
    'mean', distributionMean, ...  %m
    'cov', diag(distributionStd.^2), ...
    'totalValue', 1e9);

% discrete sampling
rng(1)
disD = obtainDistributionDiscrete('normal', ...
    'nSample', 1e2, ...
    'mean', distributionMean, ...
    'cov', diag(distributionStd.^2));
% debugging - fraction outside domain
% sum(isValidDec(disD.positionMatrix, pbeDefinition.dec, 0))

% application to pbeDefinition
pbeDefinition.seeding.time = 0;
pbeDefinition.seeding.positionMatrix = disD.positionMatrix;       % m
pbeDefinition.seeding.numberVector = ...
    particleDisC.fhandle(disD.positionMatrix) ./ ...
    disD.probabilityVector; % #
pbeDefinition.seeding.probabilityVector = disD.probabilityVector;


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

% some extra variables
initialVolume = 1/length(pbeDefinition.seeding.numberVector)*...
    computeMeasureDec(dec, 'hc', ...
      pbeDefinition.seeding.positionMatrix, 'volume') * ...
    pbeDefinition.seeding.numberVector'
% inititalCubicMillimeters = inititalVolume * 1e6

% temperature profile (1h cooling)
timeEndSlope = 1*(60*60);
temperature = @(t)(50 - 30*t/timeEndSlope*(t<timeEndSlope) ...
    - 30*(t >= timeEndSlope));

%%
% some arbitrary solubility - directly given in volume and designed to
%   increase the total volume to about 25 cm^3
solubility = @(T)(4e-5 + 2e-9*T + 5*3.14159e-9*T.^2);
tempVector = 20:0.1:50;
% volumeInSupersaturationInCm3 = (solubility(50) - solubility(30)) * (10^2)^3
% figure
% plot(tempVector, solubility(tempVector));

%%

% some arbitrary supersaturation function
initialSoluteVolume = solubility(50)
supersaturation = @(V, t)(...
    (-1*solubility(temperature(t)) + ...
     (initialSoluteVolume - V + initialVolume)...
    ) ./ solubility(temperature(t)));

%% growth rate
% y(1) will be supersaturation
kFactor = [7 6 5 4 3 2 1]'*1e-8;
aExponent = ([1 2 3 4 5 6 7]'+7)/4;
pbeDefinition.growthRate = @(H, P, y, t)(...
    kFactor .* y(1).^aExponent * ones(1, size(H,2)));
% make a plot for reasonable supersaturation
sVector = 0:0.01:2;
for i = 1:length(sVector)
    gVector(:, i) = pbeDefinition.growthRate(1, [], sVector(i), []);
end
figure
plot(sVector, gVector*1e6)
legend({'1', '2', '3', '4', '5', '6', '7'})

%% properties
pbeDefinition.pivotProperty = {'volume'};
pbeDefinition.bulkProperty = @(P, N, x, t)([...
    supersaturation(1/length(N)*P*N', t); ...
    temperature(t); ...
    solubility(temperature(t)) / 1e6]);

% time span for solution
% pbeDefinition.timeSpan = 0:30:3.5e3;          %s
% pbeDefinition.timeSpan = 0:60:2e4;          %s
pbeDefinition.timeSpan = 0:180:3e4;          %s

% options
pbeDefinition.optionStruct.distanceLimit = 3e-7; %m
pbeDefinition.optionStruct.flagDensity = 0;
pbeDefinition.optionStruct.relTol = 1e-5;
pbeDefinition.optionStruct.absTolPosition = 1e-10;
pbeDefinition.optionStruct.absTolDensity = 1e8;
pbeDefinition.optionStruct.odeSolver = @ode45;