function cdef = extendCdef(cdef, varargin)
% extend crystal definition according to symmetry
%
% Syntax: cdef = extendCdef(cdef)
%
% This function creates additional symmetry information (Herman-Mauguin 
%   symbol, Schonflies symbol, Laue class, crystal system), the complete 
%   list of Miller indices according to symmetry, the group mapping matrix 
%   and the matrix of facet normals in Cartesian space. The output 
%   structure has the fields:
%     cdef.sym.Schonflies - Schonflies symbol (symmetry)
%     cdef.sym.HermanMaugin - Herman-Mauguin symbol (symmetry)
%     cdef.sym.LaueClass - Laue class (symmetry, not unique)
%     cdef.sym.system - crystal system (symmetry, not unique)
%     cdef.millerComplete - the complete list of Miller indices
%     cdef.groupMappingMatrix - the group mapping matrix (hC -> h)
%     cdef.A - matrix of facet normals
%     cdef.mapMillerToReal - mapping matrix from Miller coordinates (row
%       vector) to real coordinates
%     cdef.mapRealToMiller - mapping matrix from real coordinates (row
%       vector) to Miller coordinates
%
% The input crystal definition must contain the fields:
%   cdef.millerDefined - the face groups of the crystal
%   cdef.sym.defined - a name of the crystal definition (Schoenflies, Laue,
%     International)
%   cdef.lat.a, cdef.lat.b, cdef.lat.c - the length of the crystallographic
%     axes
%   cdef.lat.alpha, cdef.lat.beta, cdef.lat.gamma - the angles of the
%     lattice coordinate system
%
% This function is an interface to functions provided by Christian Borchert
%   under the GNU GPL. See subfolder borchert-05-2010 for details.
%
% See also: obtainCrystal, viewCrystal, editCrystal

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

%% Parameters
borchertDir = 'borchert-05-2010';

%% Handle options
[normalize, checkDoubleFacet] = mapOption(varargin, ...
    'normalize', 1, ...
    'check double', 1);

%% Input handling
%! ToDo: Standard should be degree, no crystallographer reads radian
if strcmp(cdef.lat.unit_angle, 'degree')
    cdef.lat.unit_angle = 'radian';
    cdef.lat.alpha = cdef.lat.alpha / 180 * pi;
    cdef.lat.beta  = cdef.lat.beta  / 180 * pi;
    cdef.lat.gamma = cdef.lat.gamma / 180 * pi;
end
if ~strcmp(cdef.lat.unit_angle, 'radian')
    error(['angles in cdef.lat must be given in radian or degree (cdef.lat.unit_angle = ' cdef.lat.unit_angle]);
end

%% Add external path
% where am I
myPath = which('extendCdef');
tmp = find(myPath == filesep, 1, 'last');
myPath = myPath(1:tmp);
% where is geometry
borchertPath = [myPath borchertDir];
eval(['addpath ' borchertPath]);

%% Symmetry information
pg = pointgroup;
cdef.sym.Schoenflies = pg(cdef.sym.defined).SF;
cdef.sym.HermanMauguin = pg(cdef.sym.defined).HM;
cdef.sym.Laue = pg(cdef.sym.defined).LC;
cdef.sym.system = abbreviateCrystalSystem(pg(cdef.sym.defined).CF);
%cdef.sym.name = pg(cdef.sym.defined).TN; % implemented in german

%% Obtain complete set of Miller indices
[cdef.millerComplete, ~, facetTypeVector] = ...
    form2face(cdef.millerDefined, pg(cdef.sym.defined).HM);

%% Obtain matrix of facet normals
[transformationMatrix, ~] = transmatrix(cdef.lat.a,cdef.lat.b,cdef.lat.c,...
                            cdef.lat.alpha,cdef.lat.beta,cdef.lat.gamma);
transformationMatrixInverted = inv(transformationMatrix);
for iFacet = 1:size(cdef.millerComplete, 1)
    cdef.A(iFacet,:) = cdef.millerComplete(iFacet, :) * transformationMatrixInverted;
end
% safe transformation matrix
cdef.mapMillerToReal = transformationMatrixInverted';
cdef.mapRealToMiller = transformationMatrix';

%% Generate grouping matrix
cdef.groupMappingMatrix = zeros([size(cdef.A, 1) max(facetTypeVector)]);
for iFacet = 1:length(facetTypeVector)
    cdef.groupMappingMatrix(iFacet, facetTypeVector(iFacet)) = 1;
end

%% Remove required paths
eval(['rmpath ' borchertPath]);

%% Normalize facet normals
if normalize
    for i = 1:size(cdef.A, 1)
        cdef.A(i,:) = cdef.A(i,:) / norm(cdef.A(i,:));
    end
end

%% Check four double facets
if checkDoubleFacet
    uniqueMillerMatrix = unique(cdef.millerComplete, 'rows');
    if size(uniqueMillerMatrix, 1) ~= size(cdef.millerComplete, 1)
        error('Created Miller indices are not unique')
    end
end
