function [seeding seedingTimeVector] = prepareSeeding(seeding, pbeDef, ...
    optionStruct)
% ensure proper seeding structure
%
% THIS IS NO USER FUNCTION

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU General Public License as published by the 
%   Free Software Foundation, either version 3 of the License, or (at your 
%   option) any later version.
%
% This program is distributed in the hope that it will be useful, but 
%   WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
%   General Public License for more details.
%
% You should have received a copy of the GNU General Public License along 
%   with this program.  If not, see <http://www.gnu.org/licenses/>.

seedingTimeVector = nan(1, length(seeding));
flagDensity = pbeDef.optionStruct.density;
dec = pbeDef.dec;

%! Dev: Multiple seed populations cannot be supported since the
%    distribution density might be wrong after a second seeding event. This
%    does also affect the error estimates for integral properties, though,
%    the integrals itself should be properly calculated.
% It should be straight forward to allow multiple seeding in this
%   implementation, but it is not required in this work, so far.
% If you attempt to add seeding, please consider that multiple seeding 
%   might imply having multiple crystal types in the system that might use
%   different decomposition data. This change is straight forward, but much
%   more complicated since it affects most of the functions in this folder.
%/
if length(seeding) > 1
    error('elk:pbeSolver:wrongInput', ['Handling of multiple seed ' ...
        'populations is currently not properly supported. View notes in ' ...
        'this function for details']);
end
%/

for iSeed = 1:length(seeding)
    %% seeding time
    if isfield(seeding, 'time') && ~isempty(seeding(iSeed).time)
        seedingTimeVector(iSeed) = seeding(iSeed).time;
    elseif iSeed == 1
        seedingTimeVector(iSeed) = 0;
    else
        error('elk:pbeSolver:wrongInput', ['You need to provide the seeding ' ...
            'time with pbeDefinition.seeding(iSeed).time.']);
    end
    
    % check case when user provides disC or fhandle and disD
    if isfield(seeding, 'disD') && ...
            (isfield(seeding, 'disC') || isfield(seeding, 'fhandle'))
        seeding(iSeed).positionMatrix = seeding(iSeed).disD.positionMatrix;
        seeding(iSeed).probabilityVector = seeding(iSeed).disD.probabilityVector;
        if isfield(seeding, 'fhandle')
            thisHandle = seeding(iSeed).fhandle;
        else
            thisHandle = seeding(iSeed).disC.fhandle;
        end
        seeding(iSeed).numberVector = ...
            1/length(seeding(iSeed).probabilityVector) * ...
            thisHandle(seeding(iSeed).positionMatrix) ./ ...
            seeding(iSeed).probabilityVector;
    end
    
    % check probability vector
    if ~isfield(seeding, 'probabilityVector') || ...
            isempty(seeding(iSeed).probabilityVector')       
        % throw error if density should have been provided
        if flagDensity
            error('elk:pbeSolver:wrongInput', ['To provide ' ...
             'distribution density information, you need to provide the ' ...
             'probability density of selecting the pivots in seeding ' ...
             'population ' num2str(iSeed)]);
        end;
        % Dev: we change the probability vector to nan to ensure it cannot 
        %   be used
        seeding(iSeed).probabilityVector = ...
            seeding(iSeed).numberVector*0 + nan;
    end
    
    %% validate positionMatrix, numberVector and probabilityVector
    % enforce proper sizing
    seeding(iSeed).numberVector = seeding(iSeed).numberVector(:)';
    seeding(iSeed).probabilityVector = seeding(iSeed).probabilityVector(:)';
    % check dimension
    if size(seeding(iSeed).positionMatrix, 1) - pbeDef.nExtraCoordinate ...
            ~= pbeDef.dec.nC
        error('elk:populationBalance:wrongInput', ['Number of rows in ' ...
            ' positionMatrix of seeding(' num2str(iSeed) ') does not ' ...
            'match the dimension of the problem.']);
    end

    nSample = size(seeding(iSeed).positionMatrix, 2);
    if length(seeding(iSeed).numberVector) ~= nSample || ...
            length(seeding(iSeed).probabilityVector) ~= nSample
        error('elk:populationBalance:wrongInput', ['Elements in in ' ...
            ' numberVector or probabilityVector of seeding(' ...
             num2str(iSeed) ') does not match the number of sample points']);
    end
           
    %% ensure positions are valid (ignored for incompletes)
    if dec.isComplete
    toleranceValue = optionStruct.distanceLimit;
    if (dec.isOrdinary && ...
        any(any(dec.confinementCone.A * seeding(iSeed).positionMatrix(1:dec.nC, :) > toleranceValue)) ...
       ) || (~dec.isOrdinary && ...
        any(any(dec.confinementConeHull.A * seeding(iSeed).positionMatrix(1:dec.nC, :) > toleranceValue)) ...
       )
        error('elk:populationBalance:wrongInput', ...
            ['Not all pivots of seeding population ' num2str(iSeed) ...
             ' reside in the confinement cone']);
    end
    end
    
end