function dsp = gtCrystDSpacing(hkl, Bmat)
%     dsp = gtCrystDSpacing(hkl, Bmat)
%     --------------------------------
%     Calculates d-spacing for given hkl-s and lattice parameters.
% 
%     INPUT:
%       hkl  = <double>     set of Miller indices (3,n) or (4,n)
%       Bmat = <double>     'B matrix' - (hkl) to Cartesian matrix (3,3)
%
%     OUTPUT:
%       dsp  = <single>     row vector of d-spacings (1,n)


% Lattice vectors in real space, in a Cartesian basis

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

% Reciprocal lattice point in the Cartesian basis (a vector normal to hkl plane)
ghkl = Bmat*hkl;

% D-spacings
dsp = 1./sqrt(sum(ghkl.*ghkl,1));

% TODO: work around the single problem
%dsp = single(dsp);

end % of function


% % To be checked against:
% % Lattice base vectors in rows:
% av = [a1; a2; a3]
% % Reciprocal space vectors in coloumns:
% bv = inv(av)
% ghkl = bv*hkl
% dsp = 1/norm(ghkl)