function [dsp,crystal_system,tt]=gtCalculateDist(hkl,parameters)
% [dsp,crystal_system]=gtCalculateDist(hkl,parameters)
%      hkl <double 1x3>

if nargin < 2
   disp('Usage: [dsp,crystal]=gtCalculateDist(hkl,parameters)')
   return
end

latticepar=parameters.acq.latticepar;
spacegroup=parameters.acq.spacegroup;

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=readSpaceGroup(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;
lambda=gtConvEnergyToWavelength(parameters.acq.energy);
twoth=2*asind(lambda/(2*dsp));

if nargout == 2
    crystal_system=crystal;
end

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

if ~isfield(parameters.acq,'crystal_system')
   parameters.acq.crystal_system=crystal;
   save parameters parameters
end


end