function xyz = gtCrystHKL2Cartesian(hklsp, Bmat)
% xyz = gtCrystHKL2Cartesian(hklsp, Bmat)
%
% Given the signed hkl indices of lattice planes in an arbitrary lattice 
% and unit cell, it returns the (X,Y,Z) coordinates of the plane normals in
% a Cartesian basis.
%
% The Cartesian (orthogonal) basis (X,Y,Z) is defined such that:
%   - lattice vectors of the unit cell are a1, a2, a3 (not orthogonal)
%   - X is parallel with a1
%   - Y lies in the a1-a2 plane
%   - Z is cross(X,Y)
% 
% The definition corresponds to the one in Poulsen's 3DRXRD book, Chapter
% 3. This function essentially computes the equation Gc = B*Ghkl  or 
% xyz = Bmat*hklsp  (for column vectors !)
%
% INPUT  hklsp - set of (signed) Miller indices (3,n) or (4,n)
%        Bmat  - transformation matrix from Miller indices into Cartesian
%                coordinates (3,3)
%
% OUTPUT xyz   - normalised Cartesian coordinates of the plane normals
%

    % If hexagonal material, transform the four indices into three
    if (size(hklsp, 1) == 4)
        %hklsp = hklsp([1 2 4], :);
        hklsp = gtCrystFourIndexes2Miller(hklsp', 'plane')';
    end

    % Reciprocal lattice point in the Cartesian basis (a vector normal to hkl
    % plane):
    % The [b1 b2 b3] matrix is in principal equal to the B matrix as defined in
    % Poulsen's 3DXRD book, Chapter 3, page 26.
    xyz = Bmat * double(hklsp);

    % Normalise xyz
    xyzl = sqrt(sum(xyz .* xyz, 1));
    xyz  = xyz ./ xyzl([1 1 1], :);
end