function Vc = gtVectorLab2Cryst(Vs, g, diag)
% GTVECTORLAB2CRYST  Express lab row vector(s) to crystal CS
%
%     Vc = gtVectorLab2Cryst(Vs, g, diag)
%     -----------------------------------
%
%     INPUT:
%       Vs   = <double>   Vector(s) in the lab CS.
%       g    = <cell>     Cell containing crystal orientation matrix(ces).
%       diag = <logical>  
%
%     OUTPUT:
%       Vc  = <double>  Vector Vs expressed in the crystal CS.
%
%     Note:
%     This function can be used in 4 different modes:
%       - Vs = 1x3 matrix (1 row vector)  , g = {N}x3x3 matrices
%         output N vectors expressed in the crystal CS defined by g matrices
%
%       - Vs = Mx3 matrix (M row vectors) , g = 3x3 matrix
%         output M vectors expressed in the crystal CS defined by g matrix
%
%       - Vs = Mx3 matrix (M row vectors) , g = {N}x3x3 matrices
%         output Nx3xM matrix containing NxM row vectors, which are all the Vs
%         vectors expressed in the crystals CS defined by g matrices
%
%       - diag mode (need diag argument, otherwise previous case)
%         Vs = Mx3 matrix (M row vectors) , g = {N}3x3 matrices (M=N)
%         output M (equal to N) vectors in the crystal CS defined by g matrices
%
%     Version 001 17-10-2012 by YGuilhem

if ~exist('diag', 'var') || isempty(diag)
    diag = false;
else
    diag = logical(diag);
end

% Get number of vectors and orientation matrices
M = size(Vs, 1);
if iscell(g)
    N = length(g);
else
    N = 1;
end

if M==1 && N >= 1
    % Convert cell to array of matrices
    sam2crys = cat(3, g{:});

    % Expand Vs in 3rd dimensiodifferent n
    Vext = Vs([1 1 1], :, ones(N, 1));

    % Compute Vc and permute to get in in row vectors
    Vc = permute(sum(sam2crys.*Vext, 2), [3 1 2]);

elseif M > 1 && N == 1
    % Convert cell to array of matrices
    sam2crys = g(:, :, ones(M, 1));

    % Expand Vs in 3rd dimension
    Vext = reshape(Vs.', [1 3 M]);

    % Compute Vc and permute to get in in row vectors
    Vc = permute(sum(sam2crys.*Vext([1 1 1], :, :), 2), [3 1 2]);

elseif M == N && diag
    % Convert cell to array of matrices
    sam2crys = cat(3, g{:});

    % Expand Vs in 3rd dimension
    Vext = reshape(Vs.', [1 3 M]);

    % Compute Vc and permute to get in in row vectors
    Vc = permute(sum(sam2crys.*Vext([1 1 1], :, :), 2), [3 1 2]);

elseif M>0 && N>0
    % Convert cell to array of matrices
    sam2crys = cat(3, g{:});
    sam2crys2 = sam2crys(:, :, :, ones(M, 1));

    % Expand Vs in 3rd dimension
    Vext = reshape(Vs.', [1 3 1 M]);
    VV = Vext([1 1 1], :, ones(N, 1), :);

    tmp = sum(sam2crys2.*VV, 2);
    Vc = permute(squeeze(tmp), [2 1 3]);
    %Vc = permute(squeeze(sum(sam2crys2.*VV, 2)), [2 1 3]);

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
