function [g_twins] = gtTwinOrientations(r_vectors, twin_angle, twin_axis, phaseid, plotOn)
%     [g_twins] = gtTwinOrientations(r_vectors, twin_angle, twin_axis, phaseid, plotOn)
%     -------------------------------------------------------------------------------------------
%     use r_vectors.mat output from gtReadGrains.m
% 
%     pass in an R vector, return the g (crystal to sample transformation) plus
%     all the twin variants


grainID = r_vectors(1);
Rvec    = r_vectors(2:4);
load('parameters.mat');
spacegroup = parameters.cryst(phaseid).spacegroup;
lp = parameters.cryst(phaseid).latticepar;

% come back to Miller notation for direction
if size(twin_axis,2) == 4
    u = twin_axis(1)-twin_axis(3);
    v = twin_axis(2)-twin_axis(3);
    w = twin_axis(4);
    twin_axis = [u v w];
end

% the basic g matrix
cry2sam = Rod2g(Rvec);
g_twin = [];

% Just look for sigma 3 case (cubic)
%twin_angle=60
%twin_axis=[1 1 1]
% Just the tensile twin case (hexagonal)
%twin_angle=85
%twin_axis=[1 1 -2 0]

% Normalize twin_axis
twin_axis = twin_axis/norm(twin_axis);
% Symmetry equivilents
twin_axis = gtGetReflections(twin_axis, spacegroup);

if spacegroup==225 % fcc case

    % in crystal system, parent to twin transformation
    for ii=1:size(twin_axis, 1)
        
        rotcomp = gtMathsRotationMatrixComp(twin_axis(ii, :)','col');
        
        cry2twin = gtMathsRotationTensor(twin_angle, rotcomp);
        
        g_twins(ii).g = cry2sam*cry2twin;
        
        %    i.e.  normally,  pl = g*[h k l]'    where g comes from Rod2g
        %    here, pl = g*(twin_operation*[h k l]')
        %    so our "new" g is g*twin_operation
        
    end
    
elseif spacegroup==194 % hexagonal case
    
    % convert from hex to cart - A matrix for a UVW DIRECTION (A matrix for a
    % UVW direction)
    
    % Miller notation for directions <UVW>
    [~,Amat] = gtCrystHKL2CartesianMatrix(lp);
    
    for ii=1:size(twin_axis, 1)
        twin_axis_cart(ii, :) = Amat*twin_axis(ii, :)';
        twin_axis_cart(ii, :) = twin_axis_cart(ii, :)/norm(twin_axis_cart(ii, :));
        
        rotcomp = gtMathsRotationMatrixComp(twin_axis_cart(ii, :)','col');
        
        cry2twin = gtMathsRotationTensor(twin_angle, rotcomp);
        
        g_twins(ii).g = cry2sam*cry2twin;
        g_twins(ii).R_vector = gtMathsRotMat2Rod( g_twins(ii).g );
        
    end
    
end


if plotOn
    
    figure();
    
    for ii=1:size(g_twins,2)
        subplot(3,4,ii)
        hold on
        
        x=cry2sam*[1 0 0]';
        y=cry2sam*[0 1 0]';
        z=cry2sam*[0 0 1]';
        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:')
        
        
        x=g_twins(ii).g*[1 0 0]';
        y=g_twins(ii).g*[0 1 0]';
        z=g_twins(ii).g*[0 0 1]';
        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(gca, 'xlim', [-1 1])
        set(gca, 'ylim', [-1 1])
        set(gca, 'zlim', [-1 1])
        
    end
end

end % end of function