function [info, sigma] = gtCrystTwinTest(R1, R2, phaseID, varargin)
% GTCRYSTTWINTEST  Tests if the two orientations are from two twins
%     [info, sigma] = gtCrystTwinTest(R1, R2, [phaseID], varargin)
%     ------------------------------------------------------------
%     INPUT:
%       R1 and R2   = Rodriguez vectors of the two different orientations
%       phaseID     = current phase
%
%     OPTIONAL INPUT (parse by pairs):
%       symm        = symmetry operators from gtCrystGetSymmetryOperators {[]}
%       sigmas      = sigma values from gtCrystGetSigmas {[]}
%       merge_angle = a threshold to decide if the angle is small {5}
%       convention  = hcp crystal axes convention {'X'}
%
%     Version 001 13-11-2013 by LNervo

if (~exist('phaseID','var') || isempty(phaseID))
    phaseID = 1;
end

conf.symm   = [];
conf.sigmas = [];
conf.merge_angle = 5;
conf.convention = 'X';
conf.parameters = [];
conf = parse_pv_pairs(conf, varargin);

if (isempty(conf.parameters))
    conf.parameters = gtLoadParameters();
end

crystal_system = conf.parameters.cryst(phaseID).crystal_system;
latticepar = conf.parameters.cryst(phaseID).latticepar;
spacegroup = conf.parameters.cryst(phaseID).spacegroup;

if (isempty(conf.symm))
    conf.symm = gtCrystGetSymmetryOperators(crystal_system, spacegroup);
end
if (isempty(conf.sigmas))
    conf.sigmas = gtCrystGetSigmas(crystal_system, latticepar, conf.convention);
end

sigmas = conf.sigmas;

sigma = []; % default sigma
switch(crystal_system)
    case 'cubic'
        [mis_angle, mis_axis, mis_info] = gtDisorientation(R1.', R2.', conf.symm, ...
            'latticepar', latticepar, 'sort', 'descend', varargin{:});
        if (mis_angle > 90)
            mis_angle = 180 - mis_angle;
        end

        sigmaOr = [];
        if (~isempty(mis_axis))
            % determine sigma type
            test_angle = abs(mis_angle - sigmas(:, 2)) < sigmas(:, 6);
            test_axis  = abs(acosd(dot(repmat(mis_axis,size(sigmas,1),1), sigmas(:, 3:5), 2))) < sigmas(:, 6);
            sigmaOr = sigmas(test_angle | test_axis, :);
            sigma   = sigmas(test_angle & test_axis, 1);
        end
    case 'hexagonal'
        % In the case of hexagonal crystals
        [mis_angle, mis_axis, mis_info] = gtDisorientation(R1.', R2.', conf.symm, ...
            'latticepar', latticepar, 'convention', conf.convention, varargin{:});
        if (mis_angle > 90)
            mis_angle = 180 - mis_angle;
        end

        sigmaOr = [];
        if (~isempty(mis_axis))
            % determine sigma type
            test_angle = abs(mis_angle - sigmas(:, 2)) < sigmas(:, 7);
            test_axis  = abs(acosd(dot(repmat(mis_axis,size(sigmas,1),1), sigmas(:, 8:10), 2))) < sigmas(:, 7);
            sigmaOr = sigmas(test_angle | test_axis, :);
            sigma   = sigmas(test_angle & test_axis, :);
        end
    otherwise
        disp('Crystal system is not currently supported... Quitting')
        return
end

% add low angle (merge grains) criteria
if (mis_angle < conf.merge_angle)
    sigma = -1;
end

info.R1           = R1;
info.R2           = R2;
info.mis_angle    = mis_angle;
info.mis_axis     = mis_axis;
info.sigmaOr      = sigmaOr;
info.sigmaAnd     = sigma;
info.sigmas       = conf.sigmas;
info.merge_angle  = conf.merge_angle;
info.convention   = conf.convention;
info = gtAddMatFile(info, mis_info, true, true);

end