function pbeDefinition = casePotashAlumGrd

%% description
% This case intends to reproduce the case study by Ma and Wang in
%   "Morphological Population Balance Model in Principal Component Space".
%   Explicit references are given for almost al parameters that are used,
%   here.
pbeDefinition.name = 'potashAlumGrd';

%% crystal definition
% the crystal definition already provides molar mass and density as a
%   property of the crystal compound
pbeDefinition.cdef = obtainCrystal('potash_alum_Ma');
cdef = pbeDefinition.cdef;

%% growth rate dispersion kernel
% This equation equals equation 4 of "Crystal Growth Rate Dispersion
%   Modeling Using Morphological Population Balance" by C. Y. Ma and X. Z.
%   Wang, 2008. Figure 3 in that paper is not exactly reproduced with the
%   following parameters, but no better set of parameters was found.
% The parameters are mainly taken from Lacmann (1999) who took them from
%   the thesis of U. Zacher ("Die Kristallwachstumsdispersion in einem
%   kontinuierlichen Suspensionskristallisator")

% fixed temperature
temperature = 30 + 273.15;          % K

% constants
gasConstant = 8314*10^-3;           % J/mol/K
shearModule = 6.7*10^9;             % N/m^2
poissonRatio = 0.33;                % -
mosaicLength = 0.2*10^-6;           % m
burgerLength = 1.22*10^-9;          % m
radiusDislocation = 3*burgerLength; % m
% molar volume equals the mass of the crystal per mol of KAl(SO4)2 that is
%   bound in the crystal divided by the mass of the crystal per unit
%   volume. Note that per molecule KAl(SO4)2, there are 12 Water molecules
%   bounded in the crystal.
molarVolume = (cdef.molarMass+cdef.molarExtraMass) / cdef.density; % m^3/mol
% currently unused
boltzmannConstant = 1.3806488*10^-23 % J/K
molecularVolume = 10^-28 % m^3
molecularVolume = 0.5*10^-27 % m^3

meRatio = molarVolume / gasConstant
meRatioFixed = 1.2*meRatio
testRatio = molecularVolume / boltzmannConstant

% parameters474.39/(1.725*10^6)  % m^3/mol
mosaicAngle = 1.0*pi/180;           % radian
coeff1 = 3*shearModule*burgerLength*molarVolume*mosaicAngle / ...
    (4*pi*(1-poissonRatio)*gasConstant*temperature)/sqrt(mosaicLength);
coeff2 = 3/4/(1-poissonRatio) - ...
         log(4*pi*radiusDislocation/burgerLength);
coeff3 = 1/mosaicAngle/sqrt(mosaicLength);
sigmaEff = @(sigma, L)(sigma - coeff1./sqrt(L).*(coeff2 + log(coeff3*sqrt(L))));

% verification by graph:
% figure,hold on
% L = (1:50)*10^-6;
% plot(L*10^6, sigmaEff(0.0659, L));
% ylim([0 0.0659])

%% growth rate
growthConstant = [7.753*10^-7; ...
                  1.744*10^-6; ...
                  1.124*10^-6];
growthExponent = [1.5; 1.5; 1.5];
% the characteristic lengths are directly calculated as edge lengths, as
%   written by Ma (2008, p. 2324 bottom left). The provided equations 10
%   and 11 do not fit because:
%     * they evaluate to negative edge lengths
%     * for 100, there are three different equations for edge lengths,
%       dependent on the partition.
%     * for 110 there are even 5 equations to consider
growthRate = @(sigma, P)(...
    bsxfun(@times, growthConstant*0+1, ...
      bsxfun(@power, ... 
        [sigma + (P(2, :)*0); ...
         sigmaEff(sigma, ...
                  min([P(2,:); P(3,:)], [], 1)); ...
         sigmaEff(sigma, ...
                  min([P(3,:); P(4,:); P(5,:)], [], 1))], ...
        growthExponent*0+1 ...   
      ) ...
    ));
% verify implementation:
% growthRate(0.0659, [nan nan nan; ...
%                  1   1  4; ...
%                  1   3  3; ...
%                  1   3  4; ...
%                  1   3  4]*1e-6)

% pbe growth rate
pbeDefinition.growthRate = @(H, P, x, y, t)(...
    growthRate(y(1), P));   % m/s


%% decomposition data and polytope properties
% decomposition data
pbeDefinition.dec = obtainDec(pbeDefinition.cdef.A, 'mapping', ...
    pbeDefinition.cdef.groupMappingMatrix, 'verbose', 0);
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'volume');
% edge 1 - 100/111
% edge 2 - 100/110 = 110/100
% edge 3 - 110/110
% edge 4 - 110/111
% edge 100/100 does not exist as such (equals 100/110)
pbeDefinition.dec = addMeasureDataDec(pbeDefinition.dec, 'edgeLength', ...
    [3 1 1; 3 1 0; 2 3 1; 3 3 1]);

% pivot properties
pbeDefinition.pivotProperty = @(H, up)(...
    computeMeasureDec(pbeDefinition.dec, 'hC', {up, H}, ...
      {'volume', 'edgeLength', [3 1 1; 3 1 0; 2 3 1; 3 3 1]}, ...
      [], 0));

%% bulk states and bulk properties
% bulk state
pbeDefinition.bulkRate = @(x, y, t)([]);
pbeDefinition.bulkInitial = [];

% bulk properties
% 1) supersaturation
pbeDefinition.bulkProperty = @(P, N, x)([...
    P(1,:)*N' ...
    ]);

%% initial condition
% initial distribution: defined as in Ma (2008, appendix) and verified
%   against figure 5 in Wang (2009)
% Note that the standard deviation is such that 2% of the population is
%   non-existent if only one coordinate is observed and the geometrical
%   constraints are ignored. A much larger part of this population will
%   reside at the boundary of the valid cone.
par.mean = [9 8 8.5]'*1e-6;  % m
par.std = [4 4 4]'*1e-6;     % m
par.totalValue = 2200;       % #
disC = obtainDistributionContinuous('normal', par);

%! ToDo: to be done; best is to enforce a regular grid, here.
% pbeDefinition.seeding{1} = struct(...
%     'time', 0, ...      % s
%     'disC', disC);
% pbeDefinition.seeding{1}.optionStructDiscretization = struct(...
%     'nPivot', 1e3, ...
%     'method', 'box uniform', ...
%     'regionMin', disC.par.mean - 0.1*disC.par.std, ...
%     'regionMax', disC.par.mean + 0.1*disC.par.std);
% pbeDefinition.seeding{1}.optionStructReduction = struct(...
%     'nPivot', 1e2, ...
%     'nRerun', 10);


%% time span for solution
% Wang 2009, state 57 time steps which gives an odd spacing of 58.9286
%   seconds, while a reasonable spacing of 60 seconds gives 56 time steps.
%   Hence, a 60 seconds spacing with 56 time steps is chosen.
pbeDefinition.timeSpan = 0:60:3300;       % s