function results = gtSymmetricReflections(cryst, energy)
% GTSYMMETRICREFLECTIONS  Gets all the hkl signed reflections and calculate
%                         properly d-spacing and twotheta-type
%     results = gtSymmetricReflections(cryst, energy)
%     -----------------------------------------------
%
%     INPUT:
%       cryst.hkl             = column vector of Miller indices
%                               size 3xn or 4xn for hexagonal
%       cryst.crystal_system  = crystal system of lattice
%       cryst.latticepar      = lattice parameters
%       cryst.spacegroup      = spacegroup
%       cryst.hermann_mauguin = Hermann-Mauguin symbol
%       energy                = beam energy in keV 
%
%     OUTPUT:
%       results = structure with fields:
%              .hklsp       
%              .thetasp     
%              .dspacingsp  
%              .thetatypesp 
%              .intsp       
%              .mult        
%
%     Version 005 01-08-2013 by LNervo
%       Replaced gtCalculateDist with gtCrystDSpacing and gtCrystTheta
%
%     Version 004 10-07-2013 by LNervo
%       Keep the full list, with also {-h-k-l} reflections
%       Added check of rows/columns for hkl list
%
%     Version 003 14-03-2012 by PReischig
%       Adapted to gtCalculateDist.
%
%     Version 002 08-03-2012 by PReischig
%       Changed input format: separated 'parameters' into 'cryst' and
%       'energy'.
%
%     Version 001 27-04-2012 by AKing 
%       recalculate multiplicities here because possible errors arising
%       from how xop dat file is processed


% Get symmetry operators
symm = gtCrystGetSymmetryOperators(cryst.crystal_system, cryst.spacegroup);

hkltypes_used = double(cryst.hkl);
if isrow(hkltypes_used) % we want columns to give to gtCrystSignedHKLs
    hkltypes_used = hkltypes_used';
end
% 4 columns needed for hexagonal
if strcmpi(cryst.crystal_system,'hexagonal') && size(hkltypes_used,1)==3
    hkltypes_used(4,:) =  hkltypes_used(3,:);
    hkltypes_used(3,:) = -hkltypes_used(1,:)-hkltypes_used(2,:);
end


[~, ~, mult, hklsp, thtype]  = gtCrystSignedHKLs(hkltypes_used, symm);

intensity = [];
for ii = 1:length(hkltypes_used)
    nhkls = mult(ii)/2;
    intensity = [intensity; repmat(cryst.int(ii),nhkls,1)];
end

Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
dspacingsp = gtCrystDSpacing(hklsp', Bmat);
thetasp    = gtCrystTheta(dspacingsp, energy);


% simpler output - rows output
results.hklsp       = hklsp';
results.thetasp     = thetasp;
results.dspacingsp  = dspacingsp;
results.thetatypesp = thtype';
results.intsp       = intensity';
results.mult        = mult';

end % end of function
