gtCalculateDist.m

function [dsp,crystal_system,tt]=gtCalculateDist(hkl,parameters,savepar)
% [dsp,crystal_system]=gtCalculateDist(hkl,parameters,[savepar])
%      hkl <double 1x3>
%      parameters = parameters file
%      savepar = logical to save

if ~exist('savepar','var') || isempty(savepar)
    savepar=false;
end
if ~exist('parameters','var') || isempty(parameters)
    load parameters
end
if hkl==[0 0 0]
    savepar=false;
end
if isfield(parameters.acq,'latticepar')
    latticepar=parameters.acq.latticepar;
else
    latticepar=[1 1 1 1 1 1];
    savepar=false;
end
if isfield(parameters.acq,'spacegroup')
    spacegroup=parameters.acq.spacegroup;
else 
    disp('spacegroup is not saved in parameters.mat')
    dsp=[];
    return
end

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
        
if ~isfield(parameters.acq,'hermann_mauguin')
   hermann=gtReadSpaceGroup(spacegroup); 
end

% calculate twotheta (twoth) for different spacegroups
if between(spacegroup,1,2)
    crystal='triclinic';
    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';
    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';
    invdsp2=h^2/a^2+k^2/b^2+l^2/c^2;
end
if between(spacegroup,75,142)
    crystal='tetragonal';
    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 system
        invdsp2=4/(3*a^2) * (h^2 + k^2 + h*k) + l^2/c^2;
    elseif ~isempty(strfind(hermann,'R')) % rhombohedral system
        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';
    invdsp2=4/(3*a^2) * (h^2 + k^2 + h*k) + l^2/c^2;
end
if between(spacegroup,195,230)
    crystal='cubic';
    invdsp2=(h^2+k^2+l^2)/a^2;
end

dsp=1/(invdsp2)^0.5;
if isfield(parameters.acq,'energy')
   lambda=gtConvEnergyToWavelength(parameters.acq.energy);
else
   lambda=0; 
   savepar=false;
end
twoth=2*asind(lambda/(2*dsp));

if nargout == 2
    crystal_system=crystal;
end

if nargout == 3
    crystal_system=crystal;
    tt=twoth;
end

% saving parameters.mat
if savepar
    
    if ~isfield(parameters.acq,'crystal_system') || ...
            isempty(parameters.acq.crystal_system) || ...
            ~strcmpi(parameters.acq.crystal_system,crystal)

        parameters.acq.crystal_system=crystal;
        save parameters parameters
    end


    if ~isfield(parameters.acq,'hermann_mauguin') || ...
            isempty(parameters.acq.hermann_mauguin) || ...
            ~strcmpi(parameters.acq.hermann_mauguin,hermann)

        parameters.acq.hermann_mauguin=hermann;
        save parameters parameters
    end

end


end