function derivativeStruct = genericDerivative(solverState, time, ...
    pbeDefinition)
% compute derivative of states from structures
%
% THIS IS NO USER FUCNTION

% 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/>.

%% verbosity
% if verbose
%     time/pbeDefinition.timeSpan(end)*100
% end

%! DEV: dbugging pbe lactose case
% if time > 6000
%     warning('time is up')
% end
%/

%% initialize
time = time
nDecCoordinate = pbeDefinition.nDecCoordinate;
nExtraCoordinate = pbeDefinition.nExtraCoordinate;

%% ensure hcVectors in cone and proper upVector
% Note that pivotUnifiedPartitionVector can contain -1 when solver steps
%   have become too large. Fore measure calculation, this is corrected to
%   arbitrary partition=1, the growth rate will then be set to 0, a warning
%   displayed and the soler should correct the stepsize cause of the sudden
%   change in growth.
[pivotPositionMeasureMatrix pivotUnifiedPartitionVector fractionCorrected] = ...
    ensureProperDecData(solverState.hcMatrix, solverState.upVector, ...
      pbeDefinition);

%% pivot properties
pivotPropertyMatrix = pbeDefinition.pivotProperty(pivotPositionMeasureMatrix, ...
    max(pivotUnifiedPartitionVector, 1));

%% bulk properties
bulkPropertyVector = pbeDefinition.bulkProperty(pivotPropertyMatrix, ...
    solverState.numberVector, solverState.bulkStateVector, time);

%% physical rates
bulkRate = pbeDefinition.bulkRate(bulkPropertyVector, time);
nucleationRate = pbeDefinition.nucleationRate(bulkPropertyVector, time);
growthRateMatrix = pbeDefinition.growthRate(pivotPositionMeasureMatrix, ...
    pivotPropertyMatrix, bulkPropertyVector, time);

%% adjust rates to boundary
% The growthRateMatrix might be given in the proper embedding for improper
%   representations (users choice). This detail is handled by
%   limitGrowthRate to avoid uneccesary projections during this algorithm.
if pbeDefinition.optionStruct.limitGrowthRate
    positionDerivativeMatrix = [...
            limitGrowthRate(growthRateMatrix(1:(end-nExtraCoordinate), :), ...
            solverState.hcMatrix(1:nDecCoordinate, :), ...
            pivotUnifiedPartitionVector, ...
            pbeDefinition.dec, ...
            pbeDefinition.optionStruct.distanceLimit);...
        growthRateMatrix((end-nExtraCoordinate+1):end, :)];
else
    positionDerivativeMatrix = growthRateMatrix;
end

%% assign derivatives as solver struct
derivativeStruct.bulkStateVector = bulkRate;
derivativeStruct.hcMatrix = positionDerivativeMatrix;
derivativeStruct.flagDensity = solverState.flagDensity;
if ~isempty(positionDerivativeMatrix)
    derivativeStruct.numberVector = positionDerivativeMatrix(1,:)*0;
else
    derivativeStruct.numberVector = zeros(1, 0);
end
if solverState.flagNucleation
    derivativeStruct.numberVector(1) = nucleationRate;
    derivativeStruct.hcMatrix(:,1) = 0.5*derivativeStruct.hcMatrix(:,1);
end

%% rate of probability density
% to provide the rate of the probability density, we need to calculate the
%   divergence of the growth rate.
if solverState.flagDensity
%     disp('DENSITY')
    divergenceVector = zeros(1, length(solverState.probabilityVector));
    for iDim = 1:size(solverState.hcMatrix, 1)
        % get a state somewhat distant in h-space
        forwardSolverState = solverState;
        deltaVector = min(...
            pbeDefinition.optionStruct.gradientRelativeDelta * ...
              forwardSolverState.hcMatrix(iDim, :), ...
            pbeDefinition.optionStruct.gradientAbsoluteDelta + ...
              0*forwardSolverState.hcMatrix(iDim, :));
               
        forwardSolverState.hcMatrix(iDim, :) = ...
            forwardSolverState.hcMatrix(iDim, :) + deltaVector;
        % remove density flag for that we only need the growth rate
        forwardSolverState.flagDensity = 0;
        forwardDerivativeStruct = genericDerivative(forwardSolverState, time, ...
            pbeDefinition);
        
        % accumulate divergence
        divergenceVector = divergenceVector + ...
            (forwardDerivativeStruct.hcMatrix(iDim, :) - ...
             derivativeStruct.hcMatrix(iDim, :)...
             ) ./ deltaVector;
    end
    derivativeStruct.probabilityVector = -1*...
        solverState.probabilityVector .* divergenceVector;
end

%% asign info for unified partition updating
derivativeStruct.fractionCorrected = fractionCorrected;
derivativeStruct.partitionIndexVector = pivotUnifiedPartitionVector;

end