function g_twins = gtTwinOrientations(R_vector, twin_angle, twin_axis, spacegroup, lp, axis_in_cart, plot_variants)
% GTTWINORIENTATIONS
%     g_twins = gtTwinOrientations(R_vector, twin_angle, twin_axis, spacegroup, lp, [plot_variants])
%     ------------------------------------------------------------------------------------------------
%     By passing in the parent R_vector, return the g (sample to crystal transformation) plus
%     all the twin variants
%
%     INPUT:
%       R_vector      = <double>   Rodrigues vector (1x3)
%       twin_angle    = <double>   twin rotation angle
%       twin_axis     = <double>   hkl of twin direction (Miller/Miller
%                                  Bravais indexes)
%       spacegroup    = <double>   current crystal spacegroup
%       lp            = <double>   current crystal lattice parameters
%       plot_variants = <logical>  plot figures with parent/twin crystal axes
%                                  for each variant {false}
%
%     OUTPUT:
%       g_twins       = <struct>   twin variants info
%
%     Version 002 10-01-2013 by LNervo

    if (~exist('axis_in_cart', 'var') || isempty(axis_in_cart))
        axis_in_cart = true;
    end
    if (~exist('plot_variants', 'var') || isempty(plot_variants))
        plot_variants = false;
    end

    % come back to Miller notation for direction
    [~, crystal_system] = gtReadSpaceGroup(spacegroup);

    % Compute the orientation matrix g = sam2cry
    sam2cry = gtMathsRod2OriMat(R_vector');

    % Symmetry equivalents
    symm = gtCrystGetSymmetryOperators(crystal_system, spacegroup);

    if (strcmpi(crystal_system, 'hexagonal') && axis_in_cart)
        twin_axis = gtCart2Hex(twin_axis, lp);
    end
    [twin_axes, ~, num_axes] = gtCrystSignedHKLs(twin_axis, symm);
    num_axes_comps = size(twin_axes, 2);

    if (strcmpi(crystal_system, 'hexagonal'))
        twin_axes_hex = twin_axes;
        % convert from hex to cart
        if (num_axes_comps == 4)
            twin_axes = gtHex2Cart(gtCrystFourIndexes2Miller(twin_axes_hex, 'direction'), lp);
        else
            twin_axes = gtHex2Cart(twin_axes_hex, lp);
        end
        twin_axes_hex = gtMathsNormalizeVectorsList(twin_axes_hex);
    else
        twin_axes_hex = zeros(num_axes, 0);
    end
    % Normalize twin_axis
    twin_axes = gtMathsNormalizeVectorsList(twin_axes);

    xc = gtVectorCryst2Lab([1 0 0], sam2cry); % sam2cry = g
    yc = gtVectorCryst2Lab([0 1 0], sam2cry); % sam2cry = g
    zc = gtVectorCryst2Lab([0 0 1], sam2cry); % sam2cry = g

    for ii = num_axes:-1:1
        rotcomp = gtMathsRotationMatrixComp(twin_axes(ii, :)', 'col');
        cry2twins(:, :, ii) = gtMathsRotationTensor(twin_angle, rotcomp);
        gs(:, :, ii) = cry2twins(:, :, ii) * sam2cry;
    end

    orimats = gtMathsAngleAxis2OriMat(twin_angle(1, ones(num_axes, 1)), twin_axes');
    r_vecs = gtMathsOriMat2Rod(gs)';
    euler_vecs = gtMathsOriMat2Euler(gs)';

    xt = gtVectorCryst2Lab([1 0 0], gs);
    yt = gtVectorCryst2Lab([0 1 0], gs);
    zt = gtVectorCryst2Lab([0 0 1], gs);
    xyzt = cat(1, ...
        reshape(xt', 1, 3, []), ...
        reshape(yt', 1, 3, []), ...
        reshape(zt', 1, 3, []) );

    for ii = num_axes:-1:1
        % i.e.  normally,  [h k l] = g * pl', where g comes from
        % gtMathsRod2OriMat here, [h k l] = g * (twin_operation * pl'), so
        % our "new" g is g*twin_operation
        g_twins(ii) = struct( ...
            'axis', twin_axes(ii, :), ...
            'axis_hex', twin_axes_hex(ii, :), ...
            'angle', twin_angle, ...
            'g', gs(:, :, ii), ...
            'orimat', orimats(:, :, ii), ...
            'R_vector', r_vecs(ii, :), ...
            'eulers', euler_vecs(ii, :), ...
            'cry2twin', cry2twins(:, :, ii), ...
            'sam2cry', sam2cry, ...
            'xyz', xyzt(:, :, ii), ...
            'abc', [xc; yc; zc] );
    end

    if (plot_variants)
        hf = figure();
        ha = [];
        for ii = num_axes:-1:1
            ha(ii) = subplot(3, 4, ii);
            hold(ha(ii), 'on');

            plot3([0 xc(1)], [0 xc(2)], [0 xc(3)], 'r:');
            plot3([0 yc(1)], [0 yc(2)], [0 yc(3)], 'g:');
            plot3([0 zc(1)], [0 zc(2)], [0 zc(3)], 'b:');

            x = xt(ii, :);
            y = yt(ii, :);
            z = zt(ii, :);
            plot3([0 x(1)], [0 x(2)], [0 x(3)], 'r');
            plot3([0 y(1)], [0 y(2)], [0 y(3)], 'g');
            plot3([0 z(1)], [0 z(2)], [0 z(3)], 'b');

            set(ha(ii), 'XLim', [-1 1]);
            set(ha(ii), 'YLim', [-1 1]);
            set(ha(ii), 'ZLim', [-1 1]);
            set(get(ha(ii), 'Title'), 'String', num2str(ii), 'FontWeight', 'bold');
            set(ha(ii), 'UserData', g_twins(ii));
        end

        h = rotate3d(hf);
        set(h, 'Enable', 'on', 'RotateStyle', 'Box');

        h_link = linkprop(ha, 'View');
        setappdata(hf, 'linkprop', h_link);
        setappdata(hf, 'g_twins', g_twins);
    end
end