function mis = gtMathsMisorientation(r_vectors1, r_vectors2, symm, biunivocal, cellarray, varargin)
% GTMATHSMISORIENTATION  Computes misorientation angle and axis between all the
%                        grains from the two datasets, given the corresponding r_vectors.
%
%     mis = gtMathsMisorientation(r_vectors1, r_vectors2, [symm], [biunivocal], [cellarray], varargin)
%     ------------------------------------------------------------------------------------------------
%     INPUT:
%       r_vectors1  = <double Nx3>      R_vectors for dataset 1
%       r_vectors2  = <double Mx3>      R_vectors for dataset 2
%       symm        = <struct array>    symmetry operators list from gtCrystGetSymmetryOperators
%       biunivocal  = <logical>         if true only biunivocal matches are
%                                       performed, by means of same index into
%                                       the two datasets {false}
%       cellarray   = <logical>         output as cell array if true,
%                                       otherwise it is one structure with field arrays
%       varargin    =                   Options for gtDisorientation
%
%     OUTPUT:
%       mis = <cell>       Length of all the combinations between grains of dataset 1
%                          and 2
%          .dis_angle = <double>     misorientation angle
%          .dis_axis  = <double>     misorientation axis
%          .id1       = <double>     grain index from list 1
%          .id2       = <double>     grain index from list 2
%
%
%     Version 002 04-06-2013 by LNervo
%       Updated version


if (~exist('symm','var') || isempty(symm))
    parameters = gtLoadParameters();
    symm = gtCrystGetSymmetryOperators(parameters.cryst.crystal_system);
end
if (~exist('biunivocal','var') || isempty(biunivocal))
    biunivocal = false;
end
if (~exist('cellarray','var') || isempty(cellarray))
    cellarray = false;
end
    
if (biunivocal && (length(r_vectors1) ~= length(r_vectors2)))
    gtError('gtMathsMisorientation:wrongArgument','Different length of r_vector lists... Wrong value for the flag ''biunivocal''.')
end

if (biunivocal)
    [mis.dis_angle, mis.dis_axis, ~] = arrayfun(...
        @(num)gtDisorientation(r_vectors1(num, :).', r_vectors2(num, :).', symm, varargin{:}), ...
        1:length(r_vectors1), 'UniformOutput', false);
    mis.dis_angle = [mis.dis_angle{:}];
    mis.dis_axis  = reshape([mis.dis_axis{:}], 3, []);

    mis.id1 = 1:length(mis.dis_angle);
    mis.id2 = 1:length(mis.dis_angle);
    mis.biunivocal = biunivocal;
    mis.n_grain1   = length(r_vectors1);
    mis.n_grain2   = length(r_vectors2);
else
    if (cellarray)
        mis = cell(0,1);
        for ii = 1:size(r_vectors1, 1)
            % calculate the misorientation angle and axis between each grain from 1 with
            % all the grains from 2
            for jj = 1:size(r_vectors2, 1)
                [dis_angle, dis_axis, ~] = gtDisorientation(r_vectors1(ii,:).', r_vectors2(jj,:).', symm, varargin{:});

                mis{end+1}.dis_angle = dis_angle;
                mis{end}.dis_axis    = dis_axis;
                mis{end}.id1         = ii;
                mis{end}.id2         = jj;
            end
        end
    else
        mis.dis_angle = [];
        mis.dis_axis  = [];
        mis.id1       = [];
        mis.id2       = [];

        for ii = 1:length(r_vectors1)
            % calculate the misorientation angle and axis between each grain from 1 with
            % all the grains from 2
            [dis_angles, dis_axes, ~] = arrayfun(@(num) ...
             gtDisorientation(r_vectors1(ii,:).', r_vectors2(num,:).', symm, varargin{:}),...
                              1:size(r_vectors2,1), 'UniformOutput', false);
            dis_angles = [dis_angles{:}];
            dis_axes   = reshape([dis_axes{:}],3,[]);

            mis.dis_angle(1,:,ii)  = dis_angles;
            mis.dis_axis(3,:,ii)   = dis_axes;
            mis.id1(1,:,ii)        = repmat(ii,1,length(r_vectors2));
            mis.id2(1,:,ii)        = 1:length(r_vectors2);
        end
        
        mis.biunivocal  = biunivocal;
        mis.n_grain1    = length(r_vectors1);
        mis.n_grain2    = length(r_vectors2);
    end
end

end