function rayMatrix = computeExtremeRaysFromCone(A, zerotol)
% compute extreme rays of a cone
%
% THIS IS NO USER FUNCTION

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

%% input verification
if rank(A, zerotol) < size(A, 2)
    error('elk:decomposition:numericalError', ['Rank of input should be ' ...
        'full column rank. This should not happen.']);
end

%% closing ray
% according to Davis, 1954: "Theory of Positive Linear Dependence", a
%   closing ray can be constructed from a linear basis of the Rn. 
%   Therefore, we first need to reduce the Matrix A
tmpA = A;
iRow = 1;
while size(tmpA, 1) > size(tmpA, 2)
    thisFilter = (1:size(tmpA,1)) ~= iRow;
    if rank(tmpA(thisFilter, :), zerotol) == size(A, 2)
        tmpA = tmpA(thisFilter, :);
    else
        iRow = iRow + 1;
    end
    if iRow > size(tmpA, 1)
        error('elk:decomposition:numericalError', 'This should not happen.');
    end
end
% Now, Davis states that a positive basis (and, hence, the missing closing
%   Ray) can be constructed by taking any negative combination of the input
%   vectors:
closeRay = -1*sum(tmpA, 1);
closeRay = closeRay / norm(closeRay, 2);

%% polytope operations
closedConeHrep = struct('A', [A; closeRay], ...
                        'h', [zeros(size(A,1), 1); 1]);
closedConeVrep = convertHrepToVrep(closedConeHrep, zerotol);

%% ray extraction
rayMatrix = eliminateRows(closedConeVrep.V, ...
                          zeros(1, size(closedConeVrep.V, 2)), ...
                          zerotol);