function [eventValue, eventTerminate, eventDirection] = ...
    genericEventFunction(state, time, pbeDefinition, solverStateOriginal)
% evaluate if events have to be triggered
%
% Event list:
%  1) dissolution          = distance to inner empty polytope cone
%  2) nucleation: new bin  = binLoadMax - binLoad
%  3) nucleation: on/off   = nucleation rate
%  4) confined vector      = test for any(partition index = unconfined)
%
% 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/>.

%! Debug:
%% Debug: visualize solver time stepping
% % console output
% disp(['Event.. ..with time: ' num2str(time, 16)])
% 
% % graphical output
% persistent iCall lastTime
% if isempty(iCall)
%     iCall = 0;lastTime = 0;
% else
%     iCall = iCall + 1;
% end
% figure(1),hold on,plot(time, time-lastTime, 'b.', 'markerSize', 5)
% lastTime = time;
%/

%% Info
%! Dev:
% Unfortunatelly, it is required to know the derivative. And to identify 
%   disappearing facets, the unpacking is mandatory.
% However, this function is planned for improvement.
%/
nTotalCoordinates = pbeDefinition.nDecCoordinate + ...
                   pbeDefinition.nExtraCoordinate;
nDecCoordinates = pbeDefinition.nDecCoordinate;

%% event configuration
eventTerminate = [ 1; 1; 1; 1; solverStateOriginal.numberVector'*0 + 1];
eventDirection = [ 0;-1; 0; 0; solverStateOriginal.numberVector'*0 + 1];
eventValue = nan(size(eventTerminate));

%% generate derivative WITH nucleation activated (always)
% to track the nucleation on/off event the nucleation must be switched on
%   at all times
eventSolverState = unpackState(state, time, solverStateOriginal);
if ~solverStateOriginal.flagNucleation
    
    % in this case where nucleation is switched off, the pivot and boundary
    %   position do NOT matter. The pivots and bins do exist only so that 
    %   the nucleation rate is calculated. Hence, they can simply be filled
    %   by zeros.
    eventSolverState = addNucleation(eventSolverState, pbeDefinition);
    
    % get derivative
    derivativeState = genericDerivative(eventSolverState, time, pbeDefinition);
else
    % get derivative from original state/stuff
    derivativeState = genericDerivative(eventSolverState, time, pbeDefinition);
end

%% start / stop of nucleation
eventValue(1) = derivativeState.numberVector(1);

%% new nucleation bin
if solverStateOriginal.flagNucleation
    eventValue(2) = pbeDefinition.optionStruct.binLoad ...
                          - state(nTotalCoordinates + 1);
else
    eventValue(2) = -1;
end

%% update unified partition information
eventValue(3) = derivativeState.fractionCorrected - ...
    pbeDefinition.optionStruct.updateUnifiedLimit;

%% confinement tracking
if pbeDefinition.optionStruct.limitGrowthRate
    if any(pbeDefinition.dec.unifiedIsConfinedVector(...
            derivativeState.partitionIndexVector)) || ...
            any(derivativeState.partitionIndexVector == 0)
        eventValue(4) = 1;
    else
        eventValue(4) = -1;
    end
else
    eventValue(4) = 1;
end


%% dissolution
% If the vector grows out, the boundary cannot be mapped properly so that
%   the correction value from mapToValidCone is necessarily +Inf or -Inf.
% The event simply jumps from -0.5 to 0.5 if the corresponding facet grows
%   out.
if pbeDefinition.dec.isComplete
    emptyCone = obtainFacetValidityCone(pbeDefinition.dec, pbeDefinition.dec.nH+1);
    distanceFromEmptyConeVector = max(emptyCone.A * ...
        eventSolverState.hcMatrix(1:nDecCoordinates, :) + ...
        10*pbeDefinition.optionStruct.distanceLimit, [], 1);
    if ~isempty(pbeDefinition.extraZeroBound)
        % the distance from the extra condition must also be considered..
        %   this is an <the above> OR <the below> must be fulfilled for a
        %   particle to remain existent
        distanceFromEmptyConeVector = min(distanceFromEmptyConeVector, ...
            max(pbeDefinition.extraZeroBound * ...
            eventSolverState.hcMatrix((nDecCoordinates+1):end), [], 1));
    end
else
    % there is no dissolution for incomplete decomposition data
    distanceFromEmptyConeVector = 1+0*eventSolverState.hcMatrix(1, :);
end

if solverStateOriginal.flagNucleation
    eventValue(5:end) = distanceFromEmptyConeVector;
elseif length(eventTerminate) > 3
    eventValue(5:end) = distanceFromEmptyConeVector(3:end);
elseif length(eventTerminate) == 3
    % in this case, nucleation was switched off and no pivots exist so far
else
    error('elk:populationBalance:internalError', 'This should not happen.');
end
%! DEV: improve speed for dissolution events
% eventValue(4:end) = sign(eventValue(4:end)).*abs(eventValue(4:end)).^0.1;