function meanWidth = computeMeanWidthCore(hrep, vrep, edgeStruct)
% compute the mean width of a polytope given H-rep and V-rep
%
% Syntax: meanWidth = computeMeanWidthCore(hrep, vrep, edgeStruct)
%
% This function is no user function! It is called by computeMeanWidthVrep
%   and computeMeanWidthHrep

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

%% Algorithm
% cycle through edges and sum up the mean width. The basic equation to do
%   this can be found in [Schneider, 2008: "Convex Bodies: the
%   Brunn-Minkowski Theory", equation4.2.30]. It calculates the intrinsic
%   volume. For the three dimensional case it can be interpreted as
%   calculating 3xV(K,K,B) (a mixed volume) so that calculating the mean 
%   width multiplied by 2pi can be reformulated to:
%     Sum up over all edges: [length of the edge] x [angle between the
%       adjacent facets] / [full angle] x [area of the unit disc]
%   which then is well interpretable with the behavior of mixed volumes or
%   the Steiner formula, respectively. For the general n-dimensional case
%   the formula reads:
%     Sum over all edges: [length of edge] x [external angle of edge]
%   where the external angle is defined in a way that does not easily
%   extend in higher dimensions as calculations on ball-sections are
%   required (Schneider, 2008 - page 100)
meanWidth = 0;
for iEdge = 1:size(edgeStruct.vertexIndexMatrix, 1)
    thisEdgeLength = norm(vrep.V(edgeStruct.vertexIndexMatrix(iEdge, 1), :) - ...
                          vrep.V(edgeStruct.vertexIndexMatrix(iEdge, 2), :) ...
                          );
    % scalar product
    scalarProduct = hrep.A(edgeStruct.facetIndexMatrix(iEdge, 1), :) * ...
                    hrep.A(edgeStruct.facetIndexMatrix(iEdge, 2), :)';
    % prevent imaginary numbers from acos
    if scalarProduct > 1
        scalarProduct = 1;
    elseif scalarProduct < -1
        scalarProduct = -1;
    end
    % opening anlge
    thisAngle = acos(scalarProduct);
    % [length] * [angle] / 2 / pi * pi:
	meanWidth = meanWidth + thisEdgeLength * thisAngle / 2;
end
meanWidth = meanWidth / 2 / pi;