gtCalculateDist.m

function [dsp,theta] = gtCalculateDist(hkl,cryst,energy)

% GTCALCULATEDIST  Calculates d-spacing and the Bragg angle for a given hkl
%     [dsp,theta] = gtCalculateDist(hkl,cryst,energy)
%     -----------------------------------------------
%    
%     INPUT:
%       hkl                   = row vector of Miller indices;
%                               size 1x3 or 1x4 for hexagonal
%       cryst.latticepar      = lattice parameters
%       cryst.spacegroup      = spacegroup
%       cryst.hermann_mauguin = Hermann-Mauguin symbol
%       energy                = beam energy in keV 
% 
%     OUTPUT:
%       dsp   = d-spacing of the specified hkl family
%       theta = the Bragg angle of the specified hkl family
%
%
%     Version 006 10-12-2012 by LNervo
%       Changed output from twotheta to theta
%
%     Version 005 14-05-2012 by LNervo
%       Remove crystal_system calculation
%
%     Version 004 14-03-2012 by PReischig
%       formatting.
%
%     Version 003 08-03-2012 by PReischig
%       Seperated 'parameters' input into 'cryst' and 'energy'.
%       Bug fix: definition of variable 'hermann'.
% 
%     Version 002 29-11-2011 by LNervo
%       hkl <double 1x3> row vector!
%       parameters = parameters file
%
%

if ~exist('hkl','var') || isempty(hkl)
    hkl = input('Insert the Miller indexes as row vector: ');
end
if ~exist('energy','var') || isempty(energy)
    energy = input('Insert the beam energy in keV: ');
end

latticepar = cryst.latticepar;
spacegroup = cryst.spacegroup;
hermann    = cryst.hermann_mauguin;
lambda     = gtConvEnergyToWavelength(energy);

hkl = double(hkl);

% Lattice parameters
a     = latticepar(1);
b     = latticepar(2);
c     = latticepar(3);
alpha = latticepar(4);
beta  = latticepar(5);
gamma = latticepar(6);

if size(hkl,2)==4  % hexagonal materials
    h = hkl(1);
    k = hkl(2);
    l = hkl(4);
else
    h = hkl(1);
    k = hkl(2);
    l = hkl(3);
end


% Spacegroups

if between(spacegroup,1,2) % crystal_system = 'triclinic'; % a ~= b ~= c;  alpha ~= beta ~= gamma ~= 90
    
    Y = h^2/a^2*(sind(alpha))^2 + k^2/b^2*(sind(beta))^2 + l^2/c^2*(sind(gamma))^2;
    Z = 2*k*l/b/c*(cosd(beta) *cosd(gamma) - cosd(alpha)) + ...
        2*l*h/c/a*(cosd(gamma)*cosd(alpha) - cosd(beta))  + ...
        2*h*k/a/b*(cosd(alpha)*cosd(beta)  - cosd(gamma)) ;
    
    invdsp2 = 1/(1 - (cosd(alpha))^2 - (cosd(beta))^2 - (cosd(gamma))^2 + ...
              2*cosd(alpha)*cosd(beta)*cosd(gamma))*(Y+Z);
end

if between(spacegroup,3,15) % crystal_system = 'monoclinic'; % a ~= b ~= c;  alpha = gamma = 90 ~= beta
    
    invdsp2 = h^2/(a*sind(beta))^2 + k^2/b^2 + l^2/(c*sind(beta))^2 - ...
              (2*h*l*cosd(beta)) / (a*c*(sind(beta))^2);
end

if between(spacegroup,16,74) % crystal_system = 'orthorhombic'; % a ~= b ~= c;  alpha = beta = gamma = 90
    
    invdsp2 = h^2/a^2 + k^2/b^2 + l^2/c^2;
end

if between(spacegroup,75,142) % crystal_system = 'tetragonal'; % a = b ~= c;  alpha = beta = gamma = 90
 
    invdsp2=(h^2+k^2)/a^2+l^2/c^2;
end

if between(spacegroup,143,167) % crystal_system = 'trigonal';
    
    if ~isempty(strfind(hermann,'P')) % hexagonal lattice_system
        % a1 = a2 = a3 ~= c;  alpha = beta = 90; gamma = 120
        
        invdsp2 = 4/(3*a^2) * (h^2 + k^2 + h*k) + l^2/c^2;
    
    elseif ~isempty(strfind(hermann,'R')) % rhombohedral lattice_system
        % a = b = c;  alpha = beta = gamma ~= 90
        
        Y1 = h^2 + k^2 + h*k;
        Y2 = 2*(h*k + h*l + k*l);
        Z1 = (sind(alpha))^2;
        Z2 = (cosd(alpha))^2 - cosd(alpha);
        W  = 1 + 2*(cosd(alpha))^3 - 3*(cosd(alpha))^2;
        
        invdsp2 = 1/a^2*(Y1*Z1 + Y2*Z2)/W;
    end
end

if between(spacegroup,168,194) || spacegroup==663 % crystal_system = 'hexagonal'; % a1 = a2 = a3 ~= c;  alpha = beta = 90; gamma = 120
    
    invdsp2 = 4/(3*a^2)*(h^2 + k^2 + h*k) + l^2/c^2;
end

if between(spacegroup,195,230) % crystal_system = 'cubic'; % a = b = c;  alpha = beta = gamma = 90
    
    invdsp2 = (h^2+k^2+l^2)/a^2;
end


% d-spacing
dsp = 1/(invdsp2)^0.5;

% theta from the Bragg's law
theta = asind(lambda/(2*dsp));

end % end of function