gtCrystVector2SST.m

function [rgb, Vsst, Vsst_p, Vc_symm] = gtCrystVector2SST(Vc, crystal_system, symm, saturate)
% GTCRYSTVECTOR2SST  Moves crystal vectors to SST zone.
%
%    [rgb, Vsst, Vsst_p, Vc_symm] = gtCrystVector2SST(Vc, crystal_system, symm[, saturate])
%    --------------------------------------------------------------------------------------
%
%    INPUT:
%      Vc             = <double>     Row vectors expressed in crystal frame
%      crystal_system = <string>     Crystal system type ('cubic'|'hexagonal')
%
%    OPTIONAL INPUT:
%      symm           = <cell>       Symmetry operators (recomputed if empty)
%      saturate       = <logical>    Saturate or not the RGB output {true}
%
%    OUTPUT:
%      rgb            = <double>     RGB colors in the SST triangle
%      Vsst           = <double>     Vc vectors moved to the SST triangle
%      Vsst_p         = <double>     Vc vectors moved to the SST triangle,
%                                    projected on the equatorial plane
%      Vc_symm        = <double>     Vc vectors in the crystal frame using
%                                    symmetry
%
%    Version 002 18-01-2013 by LNervo
%
%    Version 001 15-10-2012 by YGuilhem

    if (~exist('saturate', 'var') || isempty(saturate))
        saturate = true;
    end

    % Get symmetry operators and store g3 operators into an array
    if (~exist('symm', 'var') || isempty(symm))
        symm = gtCrystGetSymmetryOperators(crystal_system);
    end
    symm = reshape(cell2mat({symm.g3}), 3, 3, []);

    % All equivalents by symm
    symm = cat(3, symm, -symm);

    % Here we want to do a vectorized call to get Vc_symm = symm * Vc
    % It is simportant to get symmetry operator index as first index
    Vc_symm = gtVectorLab2Cryst(Vc, symm, [2 3 1]);

    x = Vc_symm(:, :, 1);
    y = Vc_symm(:, :, 2);
    z = Vc_symm(:, :, 3);

    % Define angles

    % phi : rotation around 001 axis,
    %       from 100 axis to Vc vector,
    %       projected on (100,010) plane
    Vc_phi = rad2deg(atan2(y, x));

    % chi : rotation around 010 axis,
    %       from 001 axis to Vc vector,
    %       projected on (100,001) plane
    Vc_chi = rad2deg(atan2(x, z));

    % psi : angle from 001 axis to Vc vector
    Vc_psi = acosd(z);

    switch crystal_system
    case {'cubic', 'cub'}
        anglesR   = 45 - Vc_chi; % red color proportional to (45 - chi)
        minAngleR =  0;
        maxAngleR = 45;
        anglesB   = Vc_phi; % blue color proportional to phi
        minAngleB =  0;
        maxAngleB = 45;
    case {'hexagonal', 'hex'}
        anglesR   = 90 - Vc_psi; % red color proportional to (90 - psi)
        minAngleR =  0;
        maxAngleR = 90;
        anglesB   = Vc_phi; % blue color proportional to phi
        minAngleB =  0;
        maxAngleB = 30;
    otherwise
        gtError('gtCrystVector2SST:wrong_crystal_system', ...
                'Unsupported crystal system!');
    end

    % Pick equiv that is in our triangle
    iSST = anglesR>=minAngleR & anglesR<maxAngleR & anglesB>=minAngleB & anglesB<maxAngleB;

    ok = squeeze(sum(iSST, 1));
    ok_no = (ok == 0);
    ok_too = (ok > 1);

    if any(ok_no)
        warning('gtCrystVector2SST:no_indices_found', ...
            'Problem when moving to SST triangle, no indices found for vectors:');
        disp('  indices:');
        disp(find(ok_no));
        disp('  vectors:');
        disp(Vc(ok_no, :));
        iSST(1, ok_no, :) = 1;
        anglesR(1, ok_no, :) = 0;
        anglesB(1, ok_no, :) = 0;
    end
    if any(ok_too)
        warning('gtCrystVector2SST:too_many_indices_found', ...
            'Might be a problem when moving to SST triangle, too many indices found!');
        disp('  indices:');
        disp(find(ok_too));
        disp('  vectors:');
        disp(Vc(ok_too, :));
        i_too = find(iSST(:, ok_too, :));
        iSST(i_too(2:end), ok_too, :) = 0;
    end

    angleR = squeeze(anglesR(iSST));
    angleB = squeeze(anglesB(iSST));

    if nargout > 1
        Vsst = reshape(Vc_symm(iSST(:, :, [1 1 1])), [], 3);
        if nargout > 2
            denom = 1 + Vsst(:, 3);
            Vsst_p = Vsst(:, 1:2) ./ denom(:, [1 1]);
        end
    end

    % Color map
    red   = angleR/maxAngleR;
    green = (maxAngleR - angleR)/maxAngleR .* (maxAngleB - angleB)/maxAngleB;
    blue  = (maxAngleR - angleR)/maxAngleR .* angleB/maxAngleB;

    rgb = [red green blue];
    if saturate
        rgb = gtCrystSaturateSSTColor(rgb);
    end

end % end of function