gtCalculateDist.m

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

% GTCALCULATEDIST  Calculates d-spacing and twice the Bragg angle for a
%     given hkl. Also gives crystal and lattice system of the spacegroup. 
%     [dsp,twotheta,crystal_system] = 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
%       twotheta       - twice the Bragg angle of the specified hkl family
%       crystal_system - crystallographic crystal system of the specified
%                        spacegroup
%     
%     Version 004 14-03-2012 by PReischig
%       Addition of lattice system as output; 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
%
%
%     FUNCTIONS CALLED:
%[use]- gtConvEnergyToWavelength

if ~exist('hkl','var') || isempty(hkl)
    hkl = input('Insert the Miller indeces 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 = '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 = '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 = '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 = '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 = 'trigonal';
    
    if ~isempty(strfind(hermann,'P')) % hexagonal lattice 
        % 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
        % 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 % snow case
    crystal = '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 = '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;

% 2theta
twotheta = 2*asind(lambda/(2*dsp));

crystal_system = crystal;

end % end of function