function [rgb, Vsst, Vsst_p] = gtCrystVector2SST(Vc, crystal_system, symm, saturate)
% GTCRYSTVECTOR2SST  Moves crystal vectors to SST zone.
%
%    [rgb, Vsst, Vsst_p] = gtCrystVector2SST(Vc, crystal_system, symm[, saturate])
%    -------------------------------------------------------------------------------
%
%    INPUT:
%      Vc             = <double>     Row vectors expressed in crystal frame
%      crystal_system = <string>     Crystal system type ('cubic'|'hexagonal')
%      symm           = <cell_array> Symmetry operators
%
%    OPTINAL INPUT:
%      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
%
%    Version 001 15-10-2012 by YGuilhem

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

    % Normalize vectors?
    %Vc_norm = sqrt(sum(Vc.^2, 2));
    %Vc = Vc ./ Vc_norm(:, [1 1 1]);

    % All equivalents by symm
    M = size(Vc, 1);
    N = length(symm);
    for isymm=1:N
        symm{N+isymm} = -symm{isymm};
    end
    Vc_symm = gtVectorLab2Cryst(Vc, symm);

    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 = anglesR(iSST);
    angleB = anglesB(iSST);
    if nargout > 1
        Vsst = reshape(Vc_symm(iSST(:, [1 1 1], :)), [3 M]).';
        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 = sqrt(rgb);
        sat = max(rgb, [], 2);
        rgb = rgb./sat(:, [1 1 1]);
    end

end % end of function