gtCrystCalculateReflections.m

function list = gtCrystCalculateReflections(cryst, detgeo, energy, filename)
% GTCRYSTCALCULATEREFLECTIONS  Calculate the reflections list using python 
%                              from xfab library in fable
%     list = gtCrystCalculateReflections(cryst, detgeo, energy, [filename])
%     -------------------------------------------------------------------
%     List of theoretical reflections for a given crystal system and spacegroup
%
%     INPUT: 
%       cryst    = <struct>     parameters.cryst(phaseid)
%       labgeo   = <struct>     parameters.labgeo
%       energy   = <double>     parameters.acq.energy (keV)
%       filename = <string>     filename for reflections list {'reflections_*.txt'}
%
%     OUTPUT:
%       list = structure with fields (the same with 'sp' for the full lists)
%           .hkl       = reflection list <double>
%           .theta     = theta value for each reflection <double>
%           .thetatype = theta type <double>
%           .sinthl    = sind(theta)/lambda <double>
%           .dspacing  = d-spacing <double>
%           .mult      = multiplicity of each reflections <double>
%
%
%     Version 002 12-12-2012 by LNervo
%       Compact output as a structure
%
%     Version 001 03-12-2012 by LNervo


latticepar = cryst.latticepar;
sg = cryst.spacegroup;

tmp = detgeo.detanglemin/2; % theta
minsintl = min(0.0001, sind(tmp)/gtConvEnergyToWavelength(energy));
tmp = detgeo.detanglemax/2;
maxsintl = sind(tmp)/gtConvEnergyToWavelength(energy);

if ~exist('filename','var') || isempty(filename)
    filename = 'reflections_';
end


latticepar_str = sprintf('%f ',latticepar);
latticepar_str = strtrim(latticepar_str);

% python script to be used
global GT_MATLAB_HOME;
script_file = fullfile(GT_MATLAB_HOME, 'zUtil_Python', 'reflections_list.py');

% run the command for all the reflections
cmd = sprintf('%s %s %f %f %d %s %s', ...
              script_file,latticepar_str,minsintl,maxsintl,sg,'all',filename);
[~, msg] = gtPythonCommand(cmd, true); disp(msg);

% run the command for only the reflections unique
cmd = sprintf('%s %s %f %f %d %s %s', ...
              script_file,latticepar_str,minsintl,maxsintl,sg,'unique',filename);
[~, msg] = gtPythonCommand(cmd, true); disp(msg);

% read the full list of reflections
[~,Cmat] = gtReadTextFile([filename 'all.txt'],'%f %f %f %f',[1 3],true,'Delimiter','\n','CommentStyle','#');

% extracting information
hklsp = Cmat(:,1:3);
sinthlsp = Cmat(:,4);

% read the full list of reflections
[~,Cmat] = gtReadTextFile([filename 'unique.txt'],'%f %f %f %f',[1 3],true,'Delimiter','\n','CommentStyle','#');

% extracting information
hkl = Cmat(:,1:3);
sinthl = Cmat(:,4);


if strcmp(cryst.crystal_system, 'hexagonal')
    % getting 4-indexes notation
    hklsp(:,4) = hklsp(:,3);
    hklsp(:,3) = -hklsp(:,1)-hklsp(:,2);
    hkl(:,4)   = hkl(:,3);
    hkl(:,3)   = -hkl(:,1)-hkl(:,2);
    disp('  h   k   i   l     theta   thetatype  mult    d-spacing');
    disp('----------------------------------------------------------');
else
    disp('  h   k   l     theta   thetatype  mult    d-spacing');
    disp('------------------------------------------------------');
end

% extracting theta from sinthl = sin(theta)/lambda
thetasp = asind(sinthlsp*gtConvEnergyToWavelength(energy)); % fixed BUG !!! 
theta = asind(sinthl*gtConvEnergyToWavelength(energy)); % fixed BUG !!!

[~, ind, thetatypesp] = unique(sinthlsp, 'stable');
thetatype = thetatypesp(ind);
mult = [ind(2:end)',length(thetatypesp)+1]-ind';

Bmat = gtCrystHKL2CartesianMatrix(latticepar);
dspacingsp = gtCrystDSpacing(hklsp', Bmat)';
dspacing = gtCrystDSpacing(hkl', Bmat)';


list.hklsp       = hklsp';
list.hkl         = hkl';
list.thetasp     = thetasp';
list.theta       = theta';
list.sinthlsp    = sinthlsp';
list.sinthl      = sinthl';
list.dspacingsp  = dspacingsp';
list.dspacing    = dspacing';
list.thetatypesp = thetatypesp';
list.thetatype   = thetatype';
list.mult        = mult;
list.int         = ones(1,length(thetatype)); %dummy values
list.intsp       = ones(1,length(thetatypesp)); %dummy values
list.indfam      = ind';

for ii=1:length(list.thetatype)
    if ~strcmp(cryst.crystal_system, 'hexagonal')
        fprintf( '%3d %3d %3d %9.3f %5d %10d %10.3f\n',list.hkl(:,ii),list.theta(ii),list.thetatype(ii),list.mult(ii),list.dspacing(ii) );
    else
        fprintf( '%3d %3d %3d %3d %9.3f %5d %10d %10.3f\n',list.hkl(:,ii),list.theta(ii),list.thetatype(ii),list.mult(ii),list.dspacing(ii) );
    end
end

end % end of function