function [sigma3n_twin_variants, twin_variants_label, disorient_and_indcs_struct] = gtTwinCalcSigma3nVariantsQuat(N, space_group, symm, latticepar)
% GTTWINCALCSIGMA3NVARIANTSQUAT  Function to calculate the sigma 3n twin
%                                variants and their labels in TRD
%
%     [sigma3n_twin_variants, twin_variants_label] = gtTwinCalcSigma3nVariantsQuat(N)
%     ----------------------------------------------------------------------------
%     INPUT:
%       N                     = <double>   Order of sigma value.
%
%     OUTPUT:
%       sigma3n_twin_variants = <doubel>   The quaternions of the sigma 3n
%                                          twin variants. 4x2(3^N - 1)
%                                          array.
%       twin_variants_label   = <doubel>   The labels of twin variants.
%                                          Nx2(3^N - 1) array.
%       disorient_and_indcs_struct = <struct> Containing rotation angles
%                                             and orientation families of
%                                             rotation axes and the
%                                             corresponding column indices 
%                                             in "sigma3n_twin_variants" 
%                                             and "twin_twin_variants_label"
%
%     Verison 001 23-02-2021 by ZLiu, zheheng.liu@esrf.fr
%
%     Reference
%     [1] Bryan W Reed, Roger W Minich, Robert E Rudd, and Mukul Kumar. The structure of the cubic coincident site lattice rotation group. Acta Crystallographica Section A: Foundations of Crystallography,60(3):263–277, 2004.
%     [2] Heinz, A., & Neumann, P. (1991). Representation of orientation and disorientation data for cubic, hexagonal, tetragonal and orthorhombic crystals. Acta crystallographica section a: foundations of crystallography, 47(6), 780-789.
    if ~nargin || isempty(N) || N < 1
        N = 1;
    end
    if nargin < 2
        space_group = 225;
    end
    if nargin < 3
        symm = [];
    end
    fcc_space_group = [196 202 203 209 210 216 219 225 226 227 228];
    bct_space_group = 123;
    flag_isfcc = any(space_group == fcc_space_group);
    flag_isbct = any(space_group == bct_space_group);
    if flag_isfcc || flag_isbct || true
        nof_s3 = 4;
    end
    % total number of twin variants from sigma 3 to sigma 3^N equals to 4 + 4*3 + ... +4*3^(N-1) = 2(3^N - 1)
    nof_s3n = nof_s3 / (nof_s3 - 2) * ((nof_s3 - 1) ^ N - 1);
    sigma3n_twin_variants = zeros(4, nof_s3n);
    twin_variants_label = uint8(zeros(N, nof_s3n)); % label the path of sigma 3 twin variants
    % sigma 3 twin variants. Use rotating 180 degrees to substitute rotating
    % 60 degrees.
    if flag_isfcc || flag_isbct
        twin_variants_label(1, 1:4) = [1 2 3 4]; % labels of sigma 3 twin variants
        if flag_isfcc
            sigma3n_twin_variants(:, 1:4) = [[0, 1, 1, 1]', [0, 1, 1, -1]', [0, 1, -1, 1]', [0, -1, 1, 1]'] / sqrt(3); % Reed et al. (2004)
            indcs_twin_variants = [2, 3, 4; 1, 3, 4; 1, 2, 4; 1, 2, 3].'; % indices of next three possible sigma 3 twin variants for a given twin variant
        else
            misor_tri = [latticepar(1), latticepar(3)] / norm([latticepar(1), latticepar(3)]);
            sigma3n_twin_variants(:, 1:4) = [misor_tri(2) * ones(1, nof_s3); ...
                [1, -1, 0, 0; 0, 0, 1, -1] * misor_tri(1); ...
                zeros(1, nof_s3)];
            indcs_twin_variants = [1, 3, 4; 2, 3, 4; 1, 2, 3; 1, 2, 4].'; % indices of next three possible sigma 3 twin variants for a given twin variant
        end
    end
    % loop for n
    sigma3n_inds = 1:nof_s3; % sigma 3
    for n = 2:N
        % loop for the twin variants of the last sigma value to calculate the twin variants of current sigma value.
        % e.g. if the last sigma value is 3^(n-1), there are 4*3^(n-2) twin
        % variants of last sigma value. Since the number of twin variants of all the sigma values before last sigma value is
        % 2(3^(n-2) - 1), the indices of twin variants of last sigma value are from 2(3^(n-2) - 1) + 1 to
        % 2(3^(n-2) - 1) + 4*3^(n-2), i.e., from 2*3^(n-2) - 1 to 2(3^(n-1) - 1)
        % Start to calculate the twin variants of current sigma value based on
        % each twin variant of last sigma value.
        % The indices of next three possible sigma 3 twin variants of
        % this twin variant.
        sub_twin_variants_indcs = indcs_twin_variants(:, twin_variants_label(n - 1, sigma3n_inds)); % sigma 3^(n-1)
        sub_twin_variants_indcs = sub_twin_variants_indcs(:);
        sigma3n_ind_dup_3x = ones(nof_s3 - 1, 1) * sigma3n_inds;
        sigma3n_ind_dup_3x = sigma3n_ind_dup_3x(:);
        % The indices of twin variants of current sigma value are from
        % 2*3^(n-1) - 1 to 2(3^n - 1) and last_order_sigma3n_inds(2) is
        % 2*3^(n-2). 
        % Then, calculate the corresponding three twin variants and
        % store them at corresponding indices.
        nof_sigma3n_inds = length(sigma3n_inds);
        sigma3n_inds_1 = (sigma3n_inds(nof_sigma3n_inds) + 1):(sigma3n_inds(nof_sigma3n_inds) + nof_sigma3n_inds * (nof_s3 - 1)); % sigma 3^n
        sigma3n_twin_variants(:, sigma3n_inds_1) = gtMathsQuatProduct(... % q_(3^n) = q_(3^(n-1)) * q_3
            sigma3n_twin_variants(:, sigma3n_ind_dup_3x), ...
            sigma3n_twin_variants(:, sub_twin_variants_indcs));
        % Copy the old path and add the new path.
        twin_variants_label(1:n, sigma3n_inds_1) = ...
            [twin_variants_label(1 : n - 1, sigma3n_ind_dup_3x); ...
            twin_variants_label(1, sub_twin_variants_indcs)];
        sigma3n_inds = sigma3n_inds_1;
    end
    
    % Rotation angles and orientation families of rotation axes of twin
    % variants
    if nargout > 2
        disorient_and_indcs_struct = struct(...
            'rot_ang', [], ...
            'rot_axis', [], ...
            'indices', []);
        nof_TwinVar = size(sigma3n_twin_variants, 2);
        [dis_ang_TwinVar, dis_axis_TwinVar] = gtDisorientation(...
            [ones(1, nof_TwinVar); zeros(3, nof_TwinVar)], ...
            sigma3n_twin_variants, symm, 'input', 'quat', 'constrain_to_sst', 'true', 'sort', 'descend', 'mode', 'passive');
        disorient_quat = [cosd(dis_ang_TwinVar(:)'/2); bsxfun(@times, sind(dis_ang_TwinVar(:)'/2), dis_axis_TwinVar')];
        tol = 1e-2;
        indcs_TwinVar = 1:nof_TwinVar;
        nof_struct = 0;
        while ~isempty(indcs_TwinVar)
            nof_struct = nof_struct + 1;
            disorient_and_indcs_struct(nof_struct).rot_ang = dis_ang_TwinVar(indcs_TwinVar(1));
            disorient_and_indcs_struct(nof_struct).rot_axis = dis_axis_TwinVar(indcs_TwinVar(1), :);
            if numel(indcs_TwinVar) > 1
                dev_ori = gtDisorientation(...
                    disorient_quat(:, indcs_TwinVar(1)) * ones(1, numel(indcs_TwinVar) - 1), ...
                    disorient_quat(:, indcs_TwinVar(2:end)), symm, 'input', 'quat');
                tmp_indcs = [true; abs(dev_ori) < tol];
            else
                tmp_indcs = true;
            end
            disorient_and_indcs_struct(nof_struct).indices = indcs_TwinVar(tmp_indcs);
            indcs_TwinVar = indcs_TwinVar(~tmp_indcs);
        end
    end
end
