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', 'active');
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', 'mode', 'active');
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