diff --git a/zUtil_Indexter/gtIndexFwdSimGrains.m b/zUtil_Indexter/gtIndexFwdSimGrains.m
new file mode 100644
index 0000000000000000000000000000000000000000..5e5c6e44a00c70c4fcfbab0744802e9b95cef290
--- /dev/null
+++ b/zUtil_Indexter/gtIndexFwdSimGrains.m
@@ -0,0 +1,167 @@
+function grain = gtIndexFwdSimGrains(grain, parameters)
+% GTINDEXFWDSIMGRAINS Forward simulates the spot positions for all
+% reflections in all grains.
+%
+% grain = gtIndexFwdSimGrains(grain, parameters)
+%
+% ---------------------------------------------------
+%
+% INPUT
+% grain - cell array containing grains as in index.mat
+% parameters - parameters as in the parameters file
+%
+% OUTPUT
+% grain.fsim. - added field containing fwd simulation data:
+% .cu - U coordinate of center
+% .cv - U coordinate of center
+% .cw - W (image no.) coordinate of center
+% .eta - eta diffraction angle in degrees
+% .sinth - sin(theta)
+% .plref - plane normal in reference crystal (unstrained theoretical)
+% .refind - linear index of reflection in grain
+% .thetatype - thetatype
+
+%tic
+
+if ~exist('parameters','var')
+ parameters = load('parameters.mat');
+ parameters = parameters.parameters;
+end
+
+
+labgeo = parameters.labgeo;
+samgeo = parameters.samgeo;
+cryst = parameters.cryst;
+energy = parameters.acq.energy;
+omstep = 180/parameters.acq.nproj;
+
+% Scaling for Lab -> Detector transformation (mm -> pixels)
+detscaleu = 1/labgeo.pixelsizeu ;
+detscalev = 1/labgeo.pixelsizev ;
+
+% Rotation tensor components
+rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir','col');
+
+% Detector projection matrix
+Qdet = gtFedDetectorProjectionTensor(labgeo.detdiru', labgeo.detdirv', ...
+ detscaleu, detscalev);
+
+% B matrix for transformation from (hkl) to Cartesian Crystal reference
+Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
+
+% Coordinates of the plane normals in the Cartesian Crystal system
+pl_ref = gtCrystHKL2Cartesian(double(cryst.hklsp), Bmat);
+refind0 = 1:size(pl_ref,2);
+
+% D-spacings in the reference crystal for the given reflections
+dsp_ref = gtCrystDSpacing(double(cryst.hkl(:,cryst.thetatypesp)), Bmat);
+
+
+% Rotation tensor describing crystal orientation; calculate all grains
+% at once for speed-up.
+% (gtCrystRodriguesVector and gtFedRod2RotTensor (or Rod2g) used together
+% gives the transformation tensor from CRYSTAL TO SAMPLE coordinates)
+Rvec = gtIndexAllGrainValues(grain,'R_vector',[],1,1:3);
+cry2sam = gtFedRod2RotTensor(Rvec');
+
+
+% Parallel loop through grains
+parfor ii = 1:length(grain)
+
+ % Plane normals in the Sample reference
+ pl_sam = cry2sam(:,:,ii) * pl_ref;
+
+ % Impose strain on grain
+ if ~isfield(grain{ii},'strain')
+ grain{ii}.strain.strainT = NaN(3,3);
+ end
+
+ if any(isnan(grain{ii}.strain.strainT(:)))
+ %fprintf('No strain info in grain #%d. Using as measured plane normals.\n',ii)
+ pl_samd = pl_sam;
+ drel = ones(1,size(pl_sam,2));
+ else
+ % New deformed plane normals and relative elongations
+ [pl_samd, drel] = gtStrainPlaneNormals(pl_sam, grain{ii}.strain.strainT);
+ end
+
+ % d-spacing after applying the strain tensor
+ dsp_d = drel .* dsp_ref;
+
+ % Bragg angles theta in the deformed state
+ sinth = 0.5*gtConvEnergyToWavelength(energy)./dsp_d;
+
+ % The plane normals need to be brought in the Lab reference where the
+ % beam direction and rotation axis are defined.
+ % Use the Sample -> Lab orientation transformation assuming omega=0;
+ % (vector length preserved for free vectors)
+ sam2lab = [samgeo.dirx; samgeo.diry; samgeo.dirz]';
+ pl_labd = sam2lab*pl_samd;
+
+
+ % Predict omega angles: 4 for each (hkl) plane normal (the Friedel pairs are
+ % om1-om3 and om2-om4); -(hkl) plane normals are assumed not to be in the list
+ [om4, pllab4, ~, rot4, omind4] = gtFedPredictOmegaMultiple(...
+ pl_labd, sinth, labgeo.beamdir', labgeo.rotdir', rotcomp, []);
+
+ % warning('Check function here! Possible bug!')
+
+ % Delete those where no reflection occurs
+ todel = (isnan(om4(1,:)) | isnan(om4(3,:)));
+ om4(:,todel) = [];
+ omind4(:,todel) = [];
+ sinth(todel) = [];
+ pllab4(:,todel,:) = [];
+ rot4(:,:,todel,:) = [];
+
+ plref = pl_ref(:,~todel);
+ refind = refind0(~todel);
+ thtype = cryst.thetatypesp(~todel);
+ hklsp = cryst.hklsp(:,~todel);
+
+ om = reshape(om4',[],1)';
+ plref4 = [plref, plref, plref, plref];
+ sinth4 = [sinth, sinth, sinth, sinth];
+ thtype4 = [thtype, thtype, thtype, thtype];
+ refind4 = [refind, refind, refind, refind];
+
+ % Diffraction vectors in Lab and Sample
+ dveclab = gtFedPredictDiffVecMultiple(reshape(pllab4,3,[]), labgeo.beamdir');
+ %dvecsam = gtGeoLab2Sam(dveclab', om, labgeo,...
+ % parameters.samgeo, 1, 1)';
+
+ % Etend grain center
+ gc_sam = grain{ii}.center';
+ gc_sam = gc_sam(:,ones(size(dveclab,2),1));
+
+ % Absolute detector coordinates U,V in pixels
+ uv = gtFedPredictUVMultiple(rot4, dveclab, gc_sam, labgeo.detrefpos', ...
+ labgeo.detnorm', Qdet, [labgeo.detrefu; labgeo.detrefv]);
+
+ % Keep only those spots that fall on the detector
+ keep = (uv(1,:) >= 1) & (labgeo.detsizeu >= uv(1,:)) & ...
+ (uv(2,:) >= 1) & (labgeo.detsizev >= uv(2,:));
+
+ % Absolute coordinates U,V,W (pixel,pixel,image); to account for
+ % measurement errors, do not correct for the image number
+ % to be mod(om,360)/omstep
+ grain{ii}.fsim.cu = uv(1,keep);
+ grain{ii}.fsim.cv = uv(2,keep);
+ grain{ii}.fsim.cw = om(keep)/omstep;
+
+ grain{ii}.fsim.eta = gtGeoEtaFromDiffVec(dveclab(:,keep)', labgeo)';
+ grain{ii}.fsim.sinth = sinth4(keep);
+ grain{ii}.fsim.plref = plref4(:,keep);
+ grain{ii}.fsim.refind = refind4(:,keep);
+ grain{ii}.fsim.thetatype = thtype4(keep);
+ %grain{ii}.fsim.dveclab = reshape(dveclab,3,[],4);
+ %grain{ii}.fsim.dvecsam = reshape(dvecsam,3,[],4);
+
+end
+
+%toc
+
+end % of main function
+
+
+
diff --git a/zUtil_Indexter/gtIndexFwdSimPairs.m b/zUtil_Indexter/gtIndexFwdSimPairs.m
new file mode 100644
index 0000000000000000000000000000000000000000..432f53225b222d81d5d31eef56ace968a56fe073
--- /dev/null
+++ b/zUtil_Indexter/gtIndexFwdSimPairs.m
@@ -0,0 +1,272 @@
+function grain = gtIndexFwdSimPairs(grain, labgeo, samgeo, drift, ...
+ latticepar, omstep, energy)
+% GTINDEXFWDSIMPAIRS Forward simulates spot positions on the detector
+% for all spots indexed as Friedel pairs.
+%
+% grain = gtIndexFwdSimPairs(grain, labgeo, samgeo, drift, ...
+% latticepar, omstep, energy)
+%
+% ---------------------------------------------------
+%
+% INPUT
+% grain - cell array containing grains as in index.mat
+% labgeo - as in parameters.labgeo
+% samgeo - as in parameters.samgeo
+% drift - sample drifts; if empty, no drifts are considered
+% latticepar - lattice parameters (1x6)
+% omstep - omega rotational stepping (degrees/image)
+% energy - beam energy in keV
+%
+% OUTPUT
+% grain.fsimA - added field containing fwd simulation data for spots A
+% grain.fsimB - added field containing fwd simulation data for spots B
+
+
+%%%%%%%%%%%%%%%%%%%
+%% Prepare input
+%%%%%%%%%%%%%%%%%%%
+
+% B matrix for transformation from (hkl) to Cartesian Crystal reference
+Bmat = gtCrystHKL2CartesianMatrix(latticepar);
+
+% Scaling for Lab -> Detector transformation (mm -> pixels)
+detscaleu = 1/labgeo.pixelsizeu ;
+detscalev = 1/labgeo.pixelsizev ;
+
+% Rotation tensor components
+rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir','col');
+
+% Detector projection matrix
+Qdet = gtFedDetectorProjectionTensor(labgeo.detdiru', labgeo.detdirv', ...
+ detscaleu, detscalev);
+
+
+% Rotation tensor describing crystal orientation; calculate all grains
+% at once for speed-up.
+% (gtCrystRodriguesVector and gtFedRod2RotTensor (or Rod2g) used together
+% gives the transformation tensor from CRYSTAL TO SAMPLE coordinates)
+Rvec = gtIndexAllGrainValues(grain,'R_vector',[],1,1:3);
+cry2sam = gtFedRod2RotTensor(Rvec');
+
+
+
+for ii = 1:length(grain)
+
+ if ~isfield(grain{ii},'strain')
+ grain{ii}.strain.strainT = NaN(3,3);
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %% Theoretical plane normals in deformed state
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ % Coordinates of the plane normals in the Cartesian Crystal system
+ pl_cry = gtCrystHKL2Cartesian(grain{ii}.hklsp, Bmat);
+ grain{ii}.plref = pl_cry;
+
+ % Plane normals in the undrifted SAMPLE reference
+ pl_sam = cry2sam(:,:,ii) * pl_cry;
+
+ % Impose strain on grain
+ if any(isnan(grain{ii}.strain.strainT(:)))
+ %fprintf('No strain info in grain #%d. Using as measured plane normals.\n',ii)
+ pl_samd = pl_sam;
+ drel = ones(1,size(pl_sam,2));
+ else
+ % New deformed plane normals and relative elongations
+ [pl_samd, drel] = gtStrainPlaneNormals(pl_sam, grain{ii}.strain.strainT);
+ end
+
+ % D-spacings in the reference crystal for the given reflections
+ dsp_ref = gtCrystDSpacing(grain{ii}.hkl, Bmat);
+
+ % d-spacing after applying the strain tensor
+ dsp_d = drel .* dsp_ref;
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %% Apply rotational drift
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ if ~isempty(drift)
+ % In case of drifts, need to calculate the omega angles for A and
+ % B reflections separately.
+
+ grain{ii}.fsimA = sfApplyDrifts(pl_samd, dsp_d, grain{ii}.center',...
+ grain{ii}.omegaA, drift, labgeo, samgeo, rotcomp, Qdet, omstep, energy);
+
+ grain{ii}.fsimB = sfApplyDrifts(pl_samd, dsp_d, grain{ii}.center',...
+ grain{ii}.omegaB, drift, labgeo, samgeo, rotcomp, Qdet, omstep, energy);
+
+ else
+ % In case of no drifts, the omega angles for reflections A and B
+ % can be calculated together.
+
+ % Bragg angles theta in the deformed state
+ sinth_d = 0.5*gtConvEnergyToWavelength(energy)./dsp_d;
+
+ % The plane normals need to be brought in the Lab reference where the
+ % beam direction and rotation axis are defined.
+ % Use the Sample -> Lab orientation transformation assuming omega=0;
+ % (vector length preserved for free vectors)
+ sam2lab = [samgeo.dirx; samgeo.diry; samgeo.dirz]';
+ pld = sam2lab*pl_samd;
+
+ % Predict omega angles: 4 for each plane normal (the Friedel pairs are
+ % om1-om3 and om2-om4)
+ % Do it only once, for A and B at the same time.
+ [om4, pl_lab4, ~, rot4, omind4] = gtFedPredictOmegaMultiple(...
+ pld, sinth_d, labgeo.beamdir', labgeo.rotdir', rotcomp, []);
+
+ % Check if all reflections occur
+ noref = any(isnan(om4), 1);
+ if any(noref)
+ msg = sprintf('No reflection was predicted for some pairs in grain #%d', ii);
+ error(msg)
+ end
+
+
+ % Get fwd simulated data for reflections A and B.
+ % Find the correct omega index for the pairs where the difference is
+ % the smallest (shouldn't matter if grain.omega is the measured one
+ % or the one corrected for the pair).
+ grain{ii}.fsimA = sfFwdSim(grain{ii}.omegaA, grain{ii}.center', ...
+ sinth_d, om4, pl_lab4, rot4, omind4, labgeo, samgeo, Qdet, omstep);
+
+ grain{ii}.fsimB = sfFwdSim(grain{ii}.omegaB, grain{ii}.center', ...
+ sinth_d, om4, pl_lab4, rot4, omind4, labgeo, samgeo, Qdet, omstep);
+
+ end
+
+ % Extra info
+ if 1
+ grain{ii}.fsimA.drel = drel;
+ grain{ii}.fsimB.drel = drel;
+ grain{ii}.fsimA.plsam = pl_samd;
+ grain{ii}.fsimB.plsam = pl_samd;
+ end
+end
+
+
+end % of main function
+
+
+%%%%%%%%%%%%%%%%%%%%%%
+%% SUB-FUNCTIONS
+%%%%%%%%%%%%%%%%%%%%%%
+
+function fsim = sfApplyDrifts(pl, dsp, gc, om, drift, labgeo, ...
+ samgeo, rotcomp, Qdet, omstep, energy)
+% pl - plane normals before drifts
+% sinth - sin(theta) in deformed state
+% om - omega-s to interpolate the rotational drift
+% gc - grain center in Sample reference before drift
+%
+% Use the measured omega position as an approximation to calculate
+% the drifts, because the fwd simulated is not yet known. As the
+% drifts are assumed smooth and slow, this is a good approximation.
+
+ % Plane normals, grain centers and sin(theta) after drifts
+ [pld, gcd, sinthd] = gtFitApplyDrifts(pl, gc, dsp, om, drift, energy);
+
+ % The plane normals need to be brought in the Lab reference where the
+ % beam direction and rotation axis are defined.
+ % Use the Sample -> Lab rotation assuming omega=0;
+ % (samgeo.dir norm is always 1, so sam2lab will be preserve vector norms)
+ sam2lab = [samgeo.dirx; samgeo.diry; samgeo.dirz]';
+ pld = sam2lab*pld;
+
+
+ % Predict omega angles: 4 for each plane normal (the Friedel pairs are
+ % om1-om3 and om2-om4)
+ [om4, pllab4, ~, rot4, omind4] = gtFedPredictOmegaMultiple(...
+ pld, sinthd, labgeo.beamdir', labgeo.rotdir', rotcomp, []);
+
+ % Check if all reflections occur
+ noref = any(isnan(om4), 1);
+ if any(noref)
+ error('No reflection was predicted for some pairs in grain.');
+ end
+
+ % Find the correct omega index for the pairs where the difference is
+ % the smallest (shouldn't matter if grain.omega is the measured one
+ % or the one corrected for the pair).
+
+ % Get fwd simulated data for the observed reflections
+ fsim = sfFwdSim(om, gcd, ...
+ sinthd, om4, pllab4, rot4, omind4, labgeo, samgeo, Qdet, omstep);
+
+end
+
+
+function fsim = sfFwdSim(om_mes, gcsam, sinth, om4, pllab4, rot4, ...
+ omind4, labgeo, samgeo, Qdet, omstep)
+ % Provides the fwd simulation data for a given omega index set.
+
+ % Absolute differences
+ dom = om4 - om_mes([1 1 1 1],:);
+
+ % Find the indices with the smallest absolute difference; consider
+ % errors and periodicity of 360deg in omega
+ [~, ind] = min( abs( [dom; dom - 360; dom + 360] ) , [], 1);
+
+ % Predicted omega is 360deg more
+ omhigh = ( (5 <= ind) & (ind <= 8) );
+
+ % Predicted omega is 360deg less
+ omlow = ( (9 <= ind) & (ind <= 12) );
+
+ % Correct indices
+ ind(omhigh) = ind(omhigh) - 4;
+ ind(omlow) = ind(omlow) - 8;
+
+ % Get final omega angles (can be <0) using the corresponding linear
+ % indices of omind4
+ linind = ind + (0:length(ind)-1)*4;
+ fsim.om = om4(linind);
+ fsim.om(omhigh) = fsim.om(omhigh) - 360;
+ fsim.om(omlow) = fsim.om(omlow) + 360;
+
+ % Omega index
+ fsim.omind = omind4(linind);
+
+ % Omegas in a column vector
+ ind1 = (ind==1);
+ ind2 = (ind==2);
+ ind3 = (ind==3);
+ ind4 = (ind==4);
+
+ % Extract the correct pllab vectors according to omind
+ fsim.pllab = zeros(3, length(ind));
+ fsim.pllab(:,ind1) = pllab4(:,ind1,1);
+ fsim.pllab(:,ind2) = pllab4(:,ind2,2);
+ fsim.pllab(:,ind3) = pllab4(:,ind3,3);
+ fsim.pllab(:,ind4) = pllab4(:,ind4,4);
+
+ % Extract the correct rotation matrices according to omind
+ fsim.rot = zeros(3, 3, length(ind));
+ fsim.rot(:,:,ind1) = rot4(:,:,ind1,1);
+ fsim.rot(:,:,ind2) = rot4(:,:,ind2,2);
+ fsim.rot(:,:,ind3) = rot4(:,:,ind3,3);
+ fsim.rot(:,:,ind4) = rot4(:,:,ind4,4);
+
+
+ % Diffraction vector
+ fsim.dveclab = gtFedPredictDiffVecMultiple(fsim.pllab, labgeo.beamdir');
+ fsim.dvecsam = gtGeoLab2Sam(fsim.dveclab', fsim.om, labgeo, samgeo, 1, 1)';
+
+ % Make sure grain center is extended
+ gcsam = gcsam(:,ones(size(fsim.dveclab,2),1));
+
+ % Absolute detector coordinates U,V in pixels
+ uv = gtFedPredictUVMultiple(fsim.rot, fsim.dveclab, gcsam, labgeo.detrefpos', ...
+ labgeo.detnorm', Qdet, [labgeo.detrefu; labgeo.detrefv]);
+
+
+ % Absolute coordinates U,V,W (pixel,pixel,image); to account for
+ % measurement errors, do not correct for the image number
+ % to be mod(om,360)/omstep
+ fsim.uvw = [uv; fsim.om/omstep];
+
+end
+
+