function [outX outY] = obtainTimeEvolution(resultStruct, functionHandle, ...
    varargin)
% generate time evolution from arbitrary function handle
%
% The function handle is supposed to generate a scalar from a structure 
%   that contains the variables of the solver in a time-dependent state. 
%   These variables are: H, N, up, P, x, y, time

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2013 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/>.

% we can assume that result struct is NOT empty so that we can evaluate the
%   first entry and determine the dimension.
outX = resultStruct.time;
outY = evaluateHandle(functionHandle, resultStruct, 1, varargin{:});

% stretch output to foll size
nTime = length(resultStruct);
outX(nTime) = 0;
outY(1, nTime) = 0;

for iTime = 2:nTime  
    outX(iTime) = resultStruct(iTime).time;
    outY(:, iTime) = evaluateHandle(functionHandle, resultStruct, iTime, ...
        varargin{:});
end

end

function outY = evaluateHandle(functionHandle, resultStruct, iTime, varargin)
    info.i = iTime;
    info.H =  resultStruct(iTime).hcMatrix;
    info.N =  resultStruct(iTime).numberVector;
    info.up = resultStruct(iTime).upVector;
    info.P =  resultStruct(iTime).pivotPropertyMatrix;
    
    info.x =  resultStruct(iTime).bulkStateVector;
    info.y =  resultStruct(iTime).bulkPropertyVector;
    info.time = resultStruct(iTime).time;
    
    info.nPivot = size(info.H, 2);
    info.nDim = size(info.H, 1);
    if length(resultStruct(iTime).dissolvedPivotVector) >= 1
        info.nZero = resultStruct(iTime).dissolvedPivotVector(1);
    else
        info.nZero = 0;
    end
    info.nSample = info.nPivot + info.nZero;
    
    % some extra information that can be precomputed to save time in
    %   functionHandle(info) if these values are required multiple times.
    info.extra = struct([]);
    for iArgin = 1:length(varargin)
        info.extra(iArgin).data = feval(varargin{iArgin}, info);
    end
    outY = functionHandle(info);
end