function symm = gtCrystGetSymmetryOperators(crystal_system, spacegroup)
% GTCRYSTGETSYMMETRYOPERATORS  Gives the list of symmetry operators.
%
%     symm = gtCrystGetSymmetryOperators(crystal_system[, spacegroup])
%     ----------------------------------------------------------------
%     This function returns the crystal symmetry operators for different
%     crystal systems/symmetry:
%       - cubic     (sg = 225)
%       - hexagonal (sg = 194, 663)
%       - trigonal  (sg = 167)
%     With no input, everything is read from the parameters.mat file.
%
%     OPTIONAL INPUT:
%       crystal_system = <string>  Name of the crystal system
%                                  ('cubic', 'hexagonal', 'trigonal')
%       spacegroup     = <int>     Space group number (not really need by now
%
%     Note:
%       All the symmetry operators here are tabulated as in the appendix from
%       Randle & Engler, 2000
%
%     Version 002 13-11-2012 by YGuilhem yoann.guilhem@esrf.fr
%       Bugfix and minor changes
%
%     Version 001 05-10-2012 by YGuilhem yoann.guilhem@esrf.fr

if (~exist('crystal_system', 'var') || isempty(crystal_system))
    if (exist('spacegroup', 'var') && ~isempty(spacegroup))
        if between(spacegroup, 195, 230)
            crystal_system = 'cubic';
        elseif (between(spacegroup, 168, 194) || spacegroup == 663)
            crystal_system = 'hexagonal';
        elseif between(spacegroup, 143, 167)
            crystal_system = 'trigonal';
        else
            gtError('gtCrystGetSymmetryOperators:wrong_input', [ ...
                'Unimplemented spacegroup (' num2str(spacegroup) ')!']);
        end
    elseif (exist('parameters.mat', 'file'))
        parameters = gtLoadParameters();
        crystal_system = parameters.cryst.crystal_system;
    else
        gtError('gtCrystGetSymmetryOperators:wrong_input', [ ...
           'crystal_system is missing and parameters file cannot be ' ...
           'found in the current directory!']);
    end
elseif (~exist('spacegroup', 'var') || isempty(spacegroup))
    if (exist('parameters.mat', 'file'))
        parameters = gtLoadParameters();
        if isfield(parameters, 'cryst')
            spacegroup = parameters.cryst.spacegroup;
        end
    end
end

switch (crystal_system)
    case 'cubic'
        if exist('spacegroup', 'var') && ~isempty(spacegroup)
            if ~between(spacegroup, 195, 230)
                gtError('gtCrystGetSymmetryOperators:wrong_input', [ ...
                    'spacegroup (' num2str(spacegroup) ') does not ' ...
                    'match with the crystal system (' crystal_system ')!']);
            end
        end

        % Randle & Engler 2000
        symm(1).g3  = [ 1  0  0 ;...
            0  1  0 ;...
            0  0  1 ];

        symm(2).g3  = [ 0  0 -1 ;...
            0 -1  0 ;...
            -1  0  0 ];

        symm(3).g3  = [ 0  0 -1 ;...
            0  1  0 ;...
            1  0  0 ];

        symm(4).g3  = [-1  0  0 ;...
            0  1  0 ;...
            0  0 -1 ];

        symm(5).g3  = [ 0  0  1 ;...
            0  1  0 ;...
            -1  0  0 ];

        symm(6).g3  = [ 1  0  0 ;...
            0  0 -1 ;...
            0  1  0 ];

        symm(7).g3  = [ 1  0  0 ;...
            0 -1  0 ;...
            0  0 -1 ];

        symm(8).g3  = [ 1  0  0 ;...
            0  0  1 ;...
            0 -1  0 ];

        symm(9).g3  = [ 0 -1  0 ;...
            1  0  0 ;...
            0  0  1 ];

        symm(10).g3 = [-1  0  0 ;...
            0 -1  0 ;...
            0  0  1 ];

        symm(11).g3 = [ 0  1  0 ;...
            -1  0  0 ;...
            0  0  1 ];

        symm(12).g3 = [ 0  0  1 ;...
            1  0  0 ;...
            0  1  0 ];

        symm(13).g3 = [ 0  1  0 ;...
            0  0  1 ;...
            1  0  0 ];

        symm(14).g3 = [ 0  0 -1 ;...
            -1  0  0 ;...
            0  1  0 ];

        symm(15).g3 = [ 0 -1  0 ;...
            0  0  1 ;...
            -1  0  0 ];

        symm(16).g3 = [ 0  1  0 ;...
            0  0 -1 ;...
            -1  0  0 ];

        symm(17).g3 = [ 0  0 -1 ;...
            1  0  0 ;...
            0 -1  0 ];

        symm(18).g3 = [ 0  0  1 ;...
            -1  0  0 ;...
            0 -1  0 ];

        symm(19).g3 = [ 0 -1  0 ;...
            0  0 -1 ;...
            1  0  0 ];

        symm(20).g3 = [ 0  1  0 ;...
            1  0  0 ;...
            0  0 -1 ];

        symm(21).g3 = [-1  0  0 ;...
            0  0  1 ;...
            0  1  0 ];

        symm(22).g3 = [ 0  0  1 ;...
            0 -1  0 ;...
            1  0  0 ];

        symm(23).g3 = [ 0 -1  0 ;...
            -1  0  0 ;...
            0  0 -1 ];

        symm(24).g3 = [-1  0  0 ;...
            0  0 -1 ;...
            0 -1  0 ];

    case 'hexagonal'
        if exist('spacegroup', 'var') && ~isempty(spacegroup)
            if ~between(spacegroup, 168, 194) && spacegroup ~= 663
                gtError('gtCrystGetSymmetryOperators:wrong_input', [ ...
                    'spacegroup (' num2str(spacegroup) ') does not ' ...
                    'match with the crystal system (' crystal_system ')!']);
            end
        else
            spacegroup = 194;
        end
        switch spacegroup
            %case {663, 194, 225} % Heritage from Sabine function..
            % Why space group 225 was in hexagonal crystals?
            % What is 663 space group? (seems to be snow)
            case {194, 663}

                symm(1).g  = [ 1  0  0  0 ;...
                    0  1  0  0 ;...
                    0  0  1  0 ;...
                    0  0  0  1 ];

                symm(2).g  = [ 0  0  1  0 ;...
                    1  0  0  0 ;...
                    0  1  0  0 ;...
                    0  0  0  1 ];

                symm(3).g  = [ 0  1  0  0 ;...
                    0  0  1  0 ;...
                    1  0  0  0 ;...
                    0  0  0  1 ];

                symm(4).g  = [ 1  0  0  0 ;...
                    0  1  0  0 ;...
                    0  0  1  0 ;...
                    0  0  0 -1 ];

                symm(5).g  = [ 0  0  1  0 ;...
                    1  0  0  0 ;...
                    0  1  0  0 ;...
                    0  0  0 -1 ];

                symm(6).g  = [ 0  1  0  0 ;...
                    0  0  1  0 ;...
                    1  0  0  0 ;...
                    0  0  0 -1 ];

                symm(7).g =  [-1  0  0  0 ;...
                    0 -1  0  0 ;...
                    0  0 -1  0 ;...
                    0  0  0  1 ];

                symm(8).g  = [ 0  0 -1  0 ;...
                    -1  0  0  0 ;...
                    0 -1  0  0 ;...
                    0  0  0  1 ];

                symm(9).g  = [ 0 -1  0  0 ;...
                    0  0 -1  0 ;...
                    -1  0  0  0 ;...
                    0  0  0  1 ];

                symm(10).g = [-1  0  0  0 ;...
                    0 -1  0  0 ;...
                    0  0 -1  0 ;...
                    0  0  0 -1 ];

                symm(11).g = [ 0  0 -1  0 ;...
                    -1  0  0  0 ;...
                    0 -1  0  0 ;...
                    0  0  0 -1 ];

                symm(12).g = [ 0 -1  0  0 ;...
                    0  0 -1  0 ;...
                    -1  0  0  0 ;...
                    0  0  0 -1 ];

                % 3 components version, should not match with g order, i.e g(i) ~= g3(i)
                % From "Introduction to Texture Analysis", by Randle & Engler, 2000

                % Rotation operator of 180 around x to emulate mirror symmetry without turning to left-handed CS
                mirror = [ 1  0  0 ;...
                    0 -1  0 ;...
                    0  0 -1 ];
                for ii=1:6
                    a = (ii-1)*60;
                    ca = cosd(a);
                    sa = sind(a);
                    symm(ii).g3 = [ ca sa  0 ;...
                        -sa ca  0 ;...
                        0  0  1 ];
                    symm(ii+6).g3 = mirror * symm(ii).g3;
                end
            otherwise
                gtError('gtCrystGetSymmetryOperators:wrong_input', ...
                    ['Unsupported spacegroup: ', num2str(spacegroup)]);
        end

    case 'trigonal'
        if (exist('spacegroup', 'var') && ~isempty(spacegroup))
            if (~between(spacegroup, 143, 167))
                gtError('gtCrystGetSymmetryOperators:wrong_input', [ ...
                    'spacegroup (' num2str(spacegroup) ') does not ' ...
                    'match with the crystal system (' crystal_system ')!']);
            end
        else
            spacegroup = 167;
        end
        switch (spacegroup)
            case {167}
                % 6 operators, also return in 3 index format
                % however, this doesn't seem right when looking at the alumina
                % datasets... more work to do.

                % Reference
                symm(1).g = [ 1  0  0  0 ;...
                    0  1  0  0 ;...
                    0  0  1  0 ;...
                    0  0  0  1 ];

                % Rotate 120 around z
                symm(2).g = [ 0  1  0  0 ;...
                    0  0  1  0 ;...
                    1  0  0  0 ;...
                    0  0  0  1 ];

                % Rotate 240 around z
                symm(3).g = [ 0  0  1  0 ;...
                    1  0  0  0 ;...
                    0  1  0  0 ;...
                    0  0  0  1 ];

                % Rotate 180 around x
                symm(4).g = [ 1  0  0  0 ;...
                    0  0  1  0 ;...
                    0  1  0  0 ;...
                    0  0  0 -1 ];

                % Rotate 180 around y
                symm(5).g = [ 0  0  1  0 ;...
                    0  1  0  0 ;...
                    1  0  0  0 ;...
                    0  0  0 -1 ];

                % Rotate 180 around dummy
                symm(6).g = [ 0  1  0  0 ;...
                    1  0  0  0 ;...
                    0  0  1  0 ;...
                    0  0  0 -1 ];

                % 3 axis
                % Reference
                c = cosd(120);
                s = sind(120);

                symm(1).g3 = [ 1  0  0 ;...
                    0  1  0 ;...
                    0  0  1 ];

                % hkil -> ihkl : Rotate 120 around z
                symm(2).g3 = [ c  s  0 ;...
                    -s  c  0 ;...
                    0  0  1 ];

                % hkil -> kihl : Rotate 240 around z
                symm(3).g3 = [ c -s  0 ;...
                    s  c  0 ;...
                    0  0  1 ];

                % Rotate 180 around hex axis 1 (// y in cart)
                symm(4).g3 = [-1  0  0 ;...
                    0  1  0 ;...
                    0  0 -1 ];

                % Rotate 180 around hex axis 2
                symm(5).g3 = symm(2).g3 * symm(4).g3 * symm(3).g3;

                % Rotate 180 around hex axis 3
                symm(6).g3 = symm(3).g3 * symm(4).g3 * symm(2).g3;

            otherwise
                gtError('gtCrystGetSymmetryOperators:wrong_spacegroup', ...
                    ['Unsupported spacegroup: ', num2str(spacegroup)]);
        end
end % end switch crystal_system

end % end of function
