function dec = addValidityCone(dec, validityConeNormalMatrix, zerotol)
% data to identify valid hC-vectors
%
% THIS FUNCTION 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/>.

%% generate validityCone (irredundant H-rep)
% we only need unique rows
validityConeNormalMatrix = unique(validityConeNormalMatrix, 'rows');
% hrep data
%! DEV: reduceHrep did not work out for one case, it reduced 4D A matrix to
%    only 3 rows even though when computing the extreme rays from the full
%    matrix, it results in 4 properly spanning rays.
%  Hence, the reduce step is ignored which should not do much harm in
%  applications.
dec.validityCone =  struct('A', validityConeNormalMatrix, ...
           'h', zeros(size(validityConeNormalMatrix, 1), 1)...
           );
% dec.validityCone = reduceHrep(dec.validityCone);
dec.validityCone = normalizeHrep(dec.validityCone);

% vrep data (rays)
rayMatrix = computeExtremeRaysFromCone(dec.validityCone.A, zerotol);

if isempty(rayMatrix)
    %% empty validity cone
    % this happens only in the context of projectDec, where it is expected
    %   to initialize the following variables with proper sizes.
    dec.se = zeros(0, dec.nC);
    dec.nS = 0;
    dec.validityCone.seIndexVector = [];
    dec.centerRay = zeros(1, dec.nC);
    dec.seScaledToClose = dec.se;
    dec.validityCone.volume = 0;
    
else
    % enter more SE data
    dec.se = scaleToReal(rayMatrix, dec, zerotol);
    dec.nS = size(rayMatrix, 1);
    dec.validityCone.seIndexVector = 1:dec.nS;
    
    % center ray
    dec.centerRay = mean(dec.se, 1)';
    dec.centerRay = dec.centerRay / norm(dec.centerRay, 2);
    % calculate SEs as they actually close the validity cone
    dec.seScaledToClose = scaleToClose(rayMatrix, dec);
    
    % volume computation is based on V-representation so we use the structuring
    %   elements and the origin as generating points.
    dec.validityCone.volume = computeVolumeHrep(closeCone(dec.validityCone, ...
        dec.centerRay));
end
