function Vc = gtVectorLab2Cryst(Vs, g, mode)
% GTVECTORLAB2CRYST  Express laboratory row vector(s) in crystal CS.
%
%     Vc = gtVectorLab2Cryst(Vs, g[, mode])
%     -----------------------------------
%
%     INPUT:
%       Vs   = <double>  Vector(s) expressed in the lab CS.
%       g    = <double>  Array containing crystal orientation matrix(ces).
%       mode = <string>  Can be either 'diag' or a 3 coord vector {[3 1 2]},
%              <int>     This affect the output shape when multiple input 
%                        vectors and matrices, i.e. M>1 and N>1
%                        See further explanations about the output.
%
%     OUTPUT:
%       Vc   = <double>  Input vector(s) Vs expressed in the crystal CS.
%                        Shape of output depends of input size (see below).
%
%     Note:
%       This function can be used in different modes:
%       - Vs = 1x3 array (1 row vector)  , g = 3x3xN array
%         output: N vectors expressed in the crystal CS defined by g matrices
%         output size: [N 3]
%
%       - Vs = Mx3 array (M row vectors) , g = 3x3 array
%         output: M vectors expressed in the crystal CS defined by g matrix
%         output size: [M 3]
%
%       - Vs = Mx3 array (M row vectors) , g = 3x3xN array
%         output: MxN vectors expressed in the crystals CS defined by g matrices
%         output size: [M 3 N] by default, can be modified using numeric mode
%           1: input vector coordinate dimension  (3)
%           2: input orientation matrix dimension (N)
%           3: input vector index dimension       (M)
%
%       - diag mode (only if M=N, need 'diag' mode, otherwise previous case)
%         Vs = Mx3 array (M row vectors) , g = 3x3xN array (M=N)
%         output: M vectors expressed in the crystal CS defined by g matrices
%         output size: [M 3 M]
%
%     Version 001 17-10-2012 by YGuilhem

% Get number of vectors (M) and orientation matrices (N)
M = size(Vs, 1);
N = size(g, 3);

diag = false; 
if ~exist('mode', 'var') || isempty(mode)
    mode = [3 1 2];
elseif isnumeric(mode)
    mode = mode(:).';
    if size(mode, 2) ~= 3 || ~all(between(mode, ones(size(mode)), 3*ones(size(mode))))
        gtError('gtVectorLab2Cryst:wrong_mode', ...
            'Numeric mode should have 3 coordinates in [1 3]!');
    end
elseif ischar(mode) && strcmp(mode, 'diag')
    if M ~= N
        gtError('gtVectorLab2Cryst:wrong_mode', ...
            ['''diag'' mode needs as many vectors as matrices!']);
    end
    diag = true;
else
    gtError('gtVectorLab2Cryst:wrong_mode', ...
        ['Unknown mode ''' mode '''!']);
end

% Check input orientation matrices g size
if (~isreal(g))
    gtError('gtVectorLab2Cryst:wrong_input_type', ...
        'Wrong input g type, wait a 3x3xN real array!');
elseif (size(g, 1) ~= 3 || size(g, 2) ~= 3)
    gtError('gtVectorLab2Cryst:wrong_input_size', ...
        'Input array g should be sized 3x3xN!');
end

% Check input vectors size
if (size(Vs, 2) ~= 3)
    gtError('gtVectorLab2Cryst:wrong_input_size', ...
        'Input array Vs should be sized Mx3!');
end

if (M == 1 && N > 1)
    % Transpose and reshape Vs
    Vs = reshape(Vs.', [1 3 M]);

    % Multiply and sum along 2nd axis to reproduce Vc = g . Vs
    %Vc = squeeze(sum(bsxfun(@times, Vs, g), 2));
    Vc = squeeze(sum(Vs([1 1 1], :, ones(N, 1)).*g, 2)).';

elseif (M >= 1 && N == 1)
%    % Transpose g to get gT, reshape it and expand it along 1st axis
%    g = reshape(g, [1 3 3]);
%
%    % Reshape Vs
%    Vs = reshape(Vs, [M 1 3]);
%
%    % Multiply and sum along 3rd axis to reproduce Vc = g . Vs
%    Vc = sum(bsxfun(@times, g, Vs), 3);

    % Can be simply done with this formulae 
    Vc = (g * Vs.').';

elseif (M == N && diag)
    % Transpose and reshape Vc
    Vs = reshape(Vs.', [1 3 M]);

    % Multiply and sum along 2nd axis to reproduce Vc = g . Vs
    %Vc = squeeze(sum(bsxfun(@times, g, Vs), 2)).';
    Vc = squeeze(sum(g.*Vs([1 1 1], :, :), 2)).';

elseif (M > 1 && N > 1)
%     % Transpose and reshape Vs
%     Vs = reshape(Vs.', [1 3 1 M]);
%
%     % Multiply and sum along 2nd axis to reproduce Vc = g . Vs
%     Vc = permute(squeeze(sum(bsxfun(@times, g, Vs), 2)), mode);

    % Faster Version of the code in the comment
    g_temp = reshape(permute(g, [2 1 3]), 3, [])';
    Vc = g_temp * (Vs.');
    Vc = reshape(Vc, 3, size(g, 3), []);
    Vc = permute(Vc, mode);

else
    gtError('gtVectorLab2Cryst:wrong_input_size', ...
            ['Cannot handle the case where size(Vs)=' num2str(size(Vs)) ...
            ' and length(g)=' num2str(N)]);
end

end % end of function