diff --git a/zUtil_Indexter/gtINDEXTER.m b/zUtil_Indexter/gtINDEXTER.m
index a1ddd3ec9179edbad55d27948a393f58a1ce765f..c3a608534ed5d6282234722fe87cafbecbf26f2d 100644
--- a/zUtil_Indexter/gtINDEXTER.m
+++ b/zUtil_Indexter/gtINDEXTER.m
@@ -1,5 +1,11 @@
-function gtINDEXTER(phaseID, flag_update, flag_inputfile, flag_log, forcemerge, parameters)
-% FUNCTION gtINDEXTER(phaseID, flag_update, flag_inputfile, flag_log, forcemerge, parameters)
+function gtINDEXTER(phaseID, flag_update, flag_inputfile, flag_log, ...
+                    forcemerge, parameters)
+% GTINDEXTER Grain indexing algorithm Indexter main function.
+%
+% gtINDEXTER(phaseID, flag_update, flag_inputfile, flag_log, ...
+%            forcemerge, parameters)
+%
+% --------------------------------------------------------------
 %
 % Main user function for the Friedel pair indexing algorithm gtIndexter.
 % It loads the Friedel pair data from the Spotpair table in the 
@@ -388,10 +394,4 @@ fprintf('Total elapsed time: %0.0f sec\n',tEnd)
 diary off
 
 
-end % of function
-
-
-
-
-
-
+end % of function
\ No newline at end of file
diff --git a/zUtil_Indexter/gtIndexAngularCons.m b/zUtil_Indexter/gtIndexAngularCons.m
index c9341e2dc1c56049f2032834088e7d34d9cd232e..99b87bd14ac988ea61bf3fcf85976abc22178094 100644
--- a/zUtil_Indexter/gtIndexAngularCons.m
+++ b/zUtil_Indexter/gtIndexAngularCons.m
@@ -1,6 +1,6 @@
 function uniangles = gtIndexAngularCons(f1,f2,cryst)
 %
-% FUNCTION uniangles = gtIndexAngularCons2(f1,f2,cryst)
+% FUNCTION uniangles = gtIndexAngularCons(f1,f2,cryst)
 %
 % Given two lattice plane families, returns a unique list of
 % theoretical angles between plane normals of the two types in real space.
@@ -15,11 +15,11 @@ function uniangles = gtIndexAngularCons(f1,f2,cryst)
 %   cryst.Bmat
 %
 % OUTPUT
-%   uniangles - unique list of angles in degrees
+%   uniangles - unique list of angles in radians
 %
 
 % Tolerance for uniqueness in degrees
-TOL = 1e-10; 
+TOL = 1e-12; 
 
 % Thetatypes (not assuming that they go from 1..n)
 th1 = cryst.thetatype(f1);
@@ -42,7 +42,7 @@ n2 = size(v2,1);
 angles = zeros(n1*n2,1);
 
 for ii = 1:size(v1,1)
-	angles((ii-1)*n2+1:ii*n2) = gtMathsVectorsAngles(v1(ii,:),v2) ;
+	angles((ii-1)*n2+1:ii*n2) = gtMathsVectorsAnglesRad(v1(ii,:), v2, 0) ;
 end
 
 % Unique angles with a tolerance
diff --git a/zUtil_Indexter/gtIndexAngularConsMatrix.m b/zUtil_Indexter/gtIndexAngularConsMatrix.m
index 02c051712cccca44498f570af7400c73c108978e..8f010c69331870b3ac6bd166c3c4083af5b355e8 100644
--- a/zUtil_Indexter/gtIndexAngularConsMatrix.m
+++ b/zUtil_Indexter/gtIndexAngularConsMatrix.m
@@ -10,7 +10,7 @@ function ACM = gtIndexAngularConsMatrix(cryst)
 %   cryst     - as in parameters.cryst; 
 %
 % OUTPUT
-%   ACM{i,j}  - theoretical angles in degrees for any i,j family indices 
+%   ACM{i,j}  - theoretical angles in radians for any i,j family indices 
 %               in cryst; symmetric cell array of size nxn where 
 %               n = length(cryst.thetatype)
 %
diff --git a/zUtil_Indexter/gtIndexBuildGrains.m b/zUtil_Indexter/gtIndexBuildGrains.m
index d519eee402fd2c89fa2ac760e4d42c0e46317b32..d22e11d18501711a3f4c1e94839e6044dae7f07a 100644
--- a/zUtil_Indexter/gtIndexBuildGrains.m
+++ b/zUtil_Indexter/gtIndexBuildGrains.m
@@ -1,18 +1,17 @@
-function [grain, unindexed] = gtIndexBuildGrains(remaining_in, tot, ...
-                            tol_build, cryst, samenv)
-%
-% FUNCTION [grain, unindexed] = gtIndexBuildGrains(remaining_in, tot, ...
-%                               tol_build, cryst, samenv)
+function [grain, tot] = gtIndexBuildGrains(tot, tol_build, cryst, samenv)
+% GTINDEXBUILDGRAINS Finds new grains (one of the main stages of Indexter).
+% 
+% [grain, tot] = gtIndexBuildGrains(tot, tol_build, cryst, samenv)
 % 
-% 1st/3 Step of the iterative indexing strategy. Finds grains among the 
-% reflections, returns the grains with their orientation, position and
+% ----------------------------------------------------------------
+%
+% 1st step out of 3 of the iterative indexing strategy. Finds grains among
+% the reflections, returns the grains with their orientation, position and
 % statistical data on the reflections.
-% The signed hkl indices of the reflections are determined.
+% The signed hkl indices of the reflections are also determined.
 % Unindexed reflection id-s are returned in 'unindexed'.  
 %
 % INPUT
-%
-%  remaining_in - yet unindexed reflection id-s (as in tot.id)
 %  tot.         - all reflection data (including indexed and unindexed)
 %  tol_build.   - tolerances for building grains;
 %                 (as in paramers.index.strategy.b.beg or b.end)
@@ -21,31 +20,29 @@ function [grain, unindexed] = gtIndexBuildGrains(remaining_in, tot, ...
 %
 % OUTPUT
 %  grain        - grain data (cell vector)
-%  unindexed    - remaining unindexed reflection id-s
+%  tot.indexed  - updated; true for the indexed reflection ID-s
 %
 
-
-                        
+                      
 fprintf('\n BUILDING GRAINS...\n\n')
 disp(tol_build)
 disp(' ')
 tic
 
-remaining  = remaining_in;
 nof_grains = 0;
-unindexed  = [];
 grain      = [];
 
 
-% Loop until there are reflections left to check
+% Loop through all reflections
 
-while length(remaining) >= 1
+for ii = 1:length(tot.indexed)
     
-    % Current reflection to be tested
-    act = gtIndexSelectRefs(remaining(1), tot);     
+    if tot.indexed(ii)
+        continue
+    end
     
-    % Remaining candidate reflections
-    cand_rem = gtIndexSelectRefs(remaining(2:end), tot); 
+    % Current reflection to be tested
+    act = gtIndexSelectRefsReduced(ii, tot);     
     
     % Tolerance for distance between diffraction paths is set equal to the 
     % approximate grain size times a constant, but is limited to a max.
@@ -53,18 +50,22 @@ while length(remaining) >= 1
                      tol_build.distf ;
     tol_build.dist = min(tol_build.dist, tol_build.distmax) ;
     
+    % Indices to check
+    inind       = ~tot.indexed;
+    inind(1:ii) = false;
+    
     % Check pair consistency of reflections (all candidates against
     % actual).
     % 'paircons' is a vector of the linear indices of consistent 
-    % reflections in cand_rem.
-    [paircons, isecs] = gtIndexCheckPairCons(act, cand_rem, ...
-                        tol_build, samenv, cryst.ACM, tol_build.ming-1);
+    % reflections in tot.
+    [paircons, isecs] = gtIndexCheckPairCons(act, tot, ...
+                        tol_build, samenv, cryst.ACM, tol_build.ming-1, inind);
     
     % If there are enough pair consistent candidates
     if length(paircons) >= tol_build.ming-1
         
         % Load reflection data for candidates
-        cand_paircons = gtIndexSelectRefs([act.id; cand_rem.id(paircons)], tot);
+        cand_paircons = gtIndexSelectRefs([act.id; tot.id(paircons)], tot);
         
         % Find a consistent group among those candidates (check group 
         % consistency)
@@ -81,24 +82,11 @@ while length(remaining) >= 1
             fprintf('Found new grain #%d:   %s\n', nof_grains, ...
                     num2str(grain{nof_grains}.pairid))
             
-            todel = cand_paircons.id(pgroup);
-            
-        % No consistent group was found for the actual reflection. 
-        % Put it in unindexed set.
-        else
-            unindexed = [unindexed; act.id];
-            todel = act.id;
+            tot.indexed(cand_paircons.id(pgroup)) = true;           
         end
         
-    % If there are not enough consistent candidates, put actual reflection 
-    % in unindexed set
-    else
-        unindexed = [unindexed; act.id];
-        todel = act.id;
     end
-    
-    % Delete reflections which have been dealt with from remaining
-    remaining(ismember(remaining,todel)) = [];
+   
     
 end
 
@@ -108,3 +96,4 @@ toc
 
 
 end % of function
+
diff --git a/zUtil_Indexter/gtIndexCheckAllPairCons.m b/zUtil_Indexter/gtIndexCheckAllPairCons.m
index 640509607a04b592b2d3d80815c3780c1f7ca172..2a65ff003c215b14d427e1177d87e78da6cbc519 100644
--- a/zUtil_Indexter/gtIndexCheckAllPairCons.m
+++ b/zUtil_Indexter/gtIndexCheckAllPairCons.m
@@ -1,7 +1,11 @@
 function ok = gtIndexCheckAllPairCons(act, cand, tol, samenv, ACM)
-% FUNCTION  ok = gtIndexCheckAllPairCons(act, cand, tol, samenv, ACM)
+% GTINDEXCHECKALLPAIRCONS All positive pair consistency check.
 %
-% Checks pair consistency of all reflections in 'cand' against 'act'.
+% ok = gtIndexCheckAllPairCons(act, cand, tol, samenv, ACM)
+%
+% ---------------------------------------------------------
+%
+% Checks whether all reflections in 'cand' are pair consistent with 'act'.
 %
 % INPUT
 %   act.    - the actual reflection data to be tested against
@@ -60,6 +64,7 @@ end
 
 pofinters = gtMathsLinePairsIntersections([act.ca act.dir],[cand.ca cand.dir]);
 cons      = gtIndexIsPointInSample(pofinters,samenv);
+
 if ~all(cons)
 	return
 end
@@ -68,33 +73,33 @@ end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Check angle
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 % Angles between plane normals consistent?
-cons  = false(length(cand.id),1);
 
-candhklind = cand.hklind;
-candpl     = cand.pl;
-actpl      = act.pl;
+cons = false(size(cand.pl,1),1);
+
+% Angular tolerance in radians
+tol_ang_rad = tol.ang*pi/180;
 
 % Angular consistency matrix; hklind{} may have more than one element !
 ACM_reduced = ACM(act.hklind{1},:);
 
-% !!! This could be a parfor loop!
-for ii = 1:length(cand.id)
-  ACV = ACM_reduced(:,candhklind{ii});
+% Angle between actual and candidates
+ang = gtMathsVectorsAnglesRad(act.pl, cand.pl, 0);  
+
+% !!! This could be a parfor loop! But may be too short.
+for ii = 1:size(cand.pl,1)
+  ACV = ACM_reduced(:,cand.hklind{ii});
   ACV = vertcat(ACV{:});
- 
-  ang = gtMathsVectorsAngles(actpl, candpl(ii,:));  
-  cons(ii) = min(abs(ACV-ang)) < tol.ang;  
+  cons(ii) = min(abs(ACV - ang(ii))) < tol_ang_rad;  
 end
 
-
 if ~all(cons)
 	return
 end
 
-ok = true;
 
+% All pairs are consistent
+ok = true;
 
-end % of function
 
+end % of function
\ No newline at end of file
diff --git a/zUtil_Indexter/gtIndexCheckPairCons.m b/zUtil_Indexter/gtIndexCheckPairCons.m
index 622d56ac8c5cef3f58cc45245bd9b9edc0b7d509..f6cf064b41e98e5b56f6193fe8587affd8fea40f 100644
--- a/zUtil_Indexter/gtIndexCheckPairCons.m
+++ b/zUtil_Indexter/gtIndexCheckPairCons.m
@@ -1,8 +1,11 @@
 function [paircons, pofinters] = gtIndexCheckPairCons(act, cand, tol, ...
-                                                      samenv, ACM, minn)
+                                 samenv, ACM, minn, cons)
+% GTINDEXCHECKPAIRCONS Pair consistency check.
 %
-% FUNCTION  [paircons, pofinters] = gtIndexCheckPairCons(act, cand, ...
-%                                                        tol, samenv, ACM)
+% [paircons, pofinters] = gtIndexCheckPairCons(act, cand, tol, ...
+%                                              samenv, ACM, minn, cons)
+%
+% ---------------------------------------------------------------------
 %
 % Checks pair consistency of reflection 'act' against a set of others 
 % in 'cand'.
@@ -14,11 +17,12 @@ function [paircons, pofinters] = gtIndexCheckPairCons(act, cand, tol, ...
 %   samenv - sample envelope
 %   ACM    - angular consistency matrix
 %   minn   - minimum number of consistent reflections to find
+%   cons   - candidates to be checked (logical vector)
 %
 % OUTPUT
 %   paircons  - a vector of logicals; true where 'cand' is pair consistent
-%   pofinters - point of 'intersection' between diff. paths of 'act' and
-%               consistent 'cand'-s pair wise; empty if found less than 
+%   pofinters - pair-wise points of 'intersection' between diff. paths of 
+%               'act' and the consistent 'cand'-s; empty if found less than 
 %               'minn' consistent reflections
 %
 
@@ -27,47 +31,27 @@ function [paircons, pofinters] = gtIndexCheckPairCons(act, cand, tol, ...
 %% Check bbox sizes and intensity
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-% initial indices of consistent reflections 
-paircons = (1:length(cand.id))';  
-
-
 % bbox y size close enough?
-cons = (act.bbys/tol.bbys < cand.bbys) & (cand.bbys < act.bbys*tol.bbys);
-paircons(~cons) = []; 
-
-if length(paircons) < minn
-	pofinters = [];
-    return
-end
-
-
+consm(:,4) = (act.bbys/tol.bbys < cand.bbys) & (cand.bbys < act.bbys*tol.bbys);
 % bbox x size close enough?
-cons = (act.bbxs/tol.bbxs < cand.bbxs(paircons)) & (cand.bbxs(paircons) < act.bbxs*tol.bbxs);
-paircons(~cons) = []; 
-
-if length(paircons) < minn
-	pofinters = [];
-    return
-end
-
-
+consm(:,3) = (act.bbxs/tol.bbxs < cand.bbxs) & (cand.bbxs < act.bbxs*tol.bbxs);
 % integrated intensities close enough?
-cons = (act.int/tol.int < cand.int(paircons)) & (cand.int(paircons) < act.int*tol.int);
-paircons(~cons) = [];
+consm(:,2) = (act.int/tol.int < cand.int) & (cand.int < act.int*tol.int);
 
-if length(paircons) < minn
-	pofinters = [];
-    return
-end
+consm(:,1) = cons;
 
+cons = all(consm,2);
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Check distance and location
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-% spatial distance small enough?
-cons = gtMathsLinesDists([act.ca act.dir],[cand.ca(paircons,:) cand.dir(paircons,:)]) < tol.dist ;
-paircons(~cons) = []; 
+paircons = find(cons);
+
+% Spatial distance small enough?
+cons2 = gtMathsLinesDists([act.ca act.dir], ...
+                          [cand.ca(cons,:) cand.dir(cons,:)]) < tol.dist ;
+paircons(~cons2) = []; 
 
 if length(paircons) < minn
 	pofinters = [];
@@ -85,10 +69,8 @@ pofinters = gtMathsLinePairsIntersections([act.ca act.dir], ...
             [cand.ca(paircons,:) cand.dir(paircons,:)]);
 
 cons = gtIndexIsPointInSample(pofinters,samenv);
-
 pofinters(~cons,:) = [];
-paircons(~cons) = []; 
-
+paircons(~cons)    = []; 
 
 if length(paircons) < minn
 	pofinters = [];
@@ -103,24 +85,31 @@ cons  = false(length(paircons),1);
 
 candhklind  = cand.hklind(paircons);
 candpl      = cand.pl(paircons,:);
-actpl       = act.pl;
 
-% Angular consistency matrix; hklind{} may have more than one element 
-ACM_reduced = ACM(act.hklind{1},:);
+% Angular tolerance in radians 
+tol_ang_rad = tol.ang*pi/180;
+
+% Angular consistency matrix for 'act';
+% act.hklind{1} may have more than one element 
+ACMact = ACM(act.hklind{1},:);
+
+% Angles between 'act' and 'cand' in radians
+ang = gtMathsVectorsAnglesRad(act.pl, candpl, 0);
 
 % !!! This could be a parfor loop!
 for ii = 1:length(paircons)
-  ACV = ACM_reduced(:,candhklind{ii});
+  ACV = ACMact(:,candhklind{ii});
   ACV = vertcat(ACV{:});
-  
-  ang = gtMathsVectorsAngles(actpl, candpl(ii,:));
-  cons(ii)  = min(abs(ACV-ang)) < tol.ang;
+  cons(ii) = min(abs(ACV - ang(ii))) < tol_ang_rad;
 end
 
-
 pofinters(~cons,:) = [];
 paircons(~cons)    = []; 
 
+if length(paircons) < minn
+	pofinters = [];
+end
+
 
 end % of function
 
diff --git a/zUtil_Indexter/gtIndexExecuteStrategy.m b/zUtil_Indexter/gtIndexExecuteStrategy.m
index 7ad13b49be5926239826d0d3dba0a84d5124a16c..0b48a747c63181038ada49e1f66e3a1ae8b2c2b2 100644
--- a/zUtil_Indexter/gtIndexExecuteStrategy.m
+++ b/zUtil_Indexter/gtIndexExecuteStrategy.m
@@ -1,24 +1,25 @@
 function [grain, remaining] = gtIndexExecuteStrategy(strategy, ...
                               tot, samenv, cryst, grain)
+% GTINDEXEXECUTESTRATEGY Executes the iterative indexing strategy.
 %
-% FUNCTION [grain, remaining] = gtIndexExecuteStrategy(strategy, grain, ...                                tot, samenv, cryst)
-%                               tot, samenv, cryst)
+% [grain, remaining] = gtIndexExecuteStrategy(strategy, tot, samenv, ...
+%                                             cryst, grain)
+%
+% -----------------------------------------------------------------------
 %
 % Executes the iterative indexing strategy. Returns the grains that have 
 % been found and the unindexed reflections. 
 %
 % INPUT
-%
-%  strategy.    - indexing strategy (as in parameters.index.strategy)
-%  tot.         - all reflection data (including indexed and unindexed)
-%  tol_build.   - tolerances for building grains;
-%                 (as in paramers.index.strategy.b.beg or b.end)
-%  cryst.       - crystallography and other data
-%  samenv.      - sample envelope in the Sample reference
+%   strategy.    - indexing strategy (as in parameters.index.strategy)
+%   tot.         - all reflection data (including indexed and unindexed)
+%   samenv.      - sample envelope in the Sample reference
+%   cryst.       - crystallography and other data
+%   grain.       - initial grain data
 %
 % OUTPUT
-%  grain        - grain data (cell vector)
-%  remaining    - remaining unindexed reflection id-s
+%   grain.       - grain data (cell vector)
+%   remaining    - remaining unindexed reflection id-s
 %
 
 
@@ -57,13 +58,23 @@ else
     
     strfields = fieldnames(strategy.b.beg);
     for ii = 1:length(strfields)
-        dif = strategy.b.end.(strfields{ii}) - strategy.b.beg.(strfields{ii});
+        if ((strategy.b.beg.(strfields{ii}) == Inf) || ...
+            (strategy.b.end.(strfields{ii}) == Inf))
+            dif = Inf;
+        else
+            dif = strategy.b.end.(strfields{ii}) - strategy.b.beg.(strfields{ii});
+        end
         strategy.l.b.(strfields{ii}) = dif/(strategy.iter-1);
     end
     
     strfields = fieldnames(strategy.m.beg);
     for ii = 1:length(strfields)
-        dif = strategy.m.end.(strfields{ii}) - strategy.m.beg.(strfields{ii});
+        if ((strategy.m.beg.(strfields{ii}) == Inf) || ...
+            (strategy.m.end.(strfields{ii}) == Inf))
+            dif = Inf;
+        else
+            dif = strategy.m.end.(strfields{ii}) - strategy.m.beg.(strfields{ii});
+        end
         strategy.l.m.(strfields{ii}) = dif/(strategy.iter-1);
     end
     
@@ -104,12 +115,6 @@ end
 
 nof_refs = length(tot.pairid);
 
-%grain = [];
-
-% List of all unindexed reflections
-remaining = 1:length(tot.pairid);
-
-
 for ii = 1:length(strategy.proc)
         
     fprintf('\n------------  STEP %d/%d  ------------\n', ...
@@ -121,7 +126,7 @@ for ii = 1:length(strategy.proc)
             
             % Update the build and merge tolerances.
             
-            fprintf('\n\n LOOSENING TOLERANCES...\n')
+            fprintf('\n LOOSENING TOLERANCES...\n')
             
             strfields = fieldnames(strategy.b.beg);
             
@@ -143,11 +148,9 @@ for ii = 1:length(strategy.proc)
         case 'b'  % BUILD GRAINS
             
             % Find new grains among unindexed reflections.
-            %   grainB     - newly created grains
-            %   remaining  - remaining unindexed reflection id-s
-            
-            [grainB, remaining] = gtIndexBuildGrains(remaining, tot, ...
-                                  tol.b, cryst, samenv);
+            %   grainB      - newly created grains
+            %   tot.indexed - updated
+            [grainB, tot] = gtIndexBuildGrains(tot, tol.b, cryst, samenv);
             
             grain = [grain, grainB];
             
@@ -155,35 +158,27 @@ for ii = 1:length(strategy.proc)
                 grain{jj}.id = jj;
             end
          
-            cfCheckOutput(grain, remaining, nof_refs)
-            
+            cfCheckOutput(grain, sum(~tot.indexed), nof_refs)
             
             
         case 'm'  % MERGE GRAINS
             
             % Find and merge identical grains that have been split.
-            %   grain  - all grain data (will be updated) 
-            
-            if ~isempty(grain)
-                [grain, remaining] = gtIndexMergeGrains(grain, remaining, ... 
-                                     tot, tol.m, cryst);
-            end
+            %   grain       - all grain data (will be updated)
+            %   tot.indexed - updated
+            [grain, tot] = gtIndexMergeGrains(grain, tot, tol.m, cryst);
             
-            cfCheckOutput(grain, remaining, nof_refs)
+            cfCheckOutput(grain, sum(~tot.indexed), nof_refs)
           
             
         case 'r'  % REFINE GRAINS
             
             % Add unindexed reflections to the grains.
-            %   grain      - all grain data (will be updated)
-            %   remaining  - remaining unindexed reflection id-s
+            %   grain       - all grain data (will be updated)
+            %   tot.indexed - updated
+            [grain, tot] = gtIndexRefineGrains(grain, tot, tol.r, cryst);
             
-            if ~isempty(grain)
-                [grain, remaining] = gtIndexRefineGrains(grain, ...
-                                     remaining, tot, tol.r, cryst);
-            end
-            
-            cfCheckOutput(grain, remaining, nof_refs)
+            cfCheckOutput(grain, sum(~tot.indexed), nof_refs)
             
     end
     
@@ -193,6 +188,8 @@ end
 fprintf('\n------------  END OF ITERATION  ------------\n')
 
 
+remaining = find(~tot.indexed);
+
 
 end % of main function
 
@@ -204,7 +201,7 @@ end % of main function
 
 % Checks if the number of indexed and unindexed reflection is consistent. 
 
-function cfCheckOutput(grain, remaining, nref_orig)
+function cfCheckOutput(grain, nref_remain, nref_orig)
         
     nof_grains = length(grain);
     
@@ -217,7 +214,7 @@ function cfCheckOutput(grain, remaining, nref_orig)
     end
     
     nref_indexed = length(allgrains.pairid);
-    nref_remain  = length(remaining);
+    %nref_remain  = length(remaining);
     nref_sum     = nref_indexed + nref_remain;
     
     fprintf('\nNumber of ...')
diff --git a/zUtil_Indexter/gtIndexFindGroupInCand.m b/zUtil_Indexter/gtIndexFindGroupInCand.m
index 50c1afeff87bf47c99f68ce08a2900e3780b82a6..6f84bb449da882a098012d11eaa8534e05f8559f 100644
--- a/zUtil_Indexter/gtIndexFindGroupInCand.m
+++ b/zUtil_Indexter/gtIndexFindGroupInCand.m
@@ -1,28 +1,34 @@
 function [goods, grain] = gtIndexFindGroupInCand(cand, isecs, tol_build,...
                                                  samenv, cryst)    
-% This is the function that finds possible new grains in the 
-% indexing process.
-% It searches for a possible grain in the input set of candidate 
+% GTINDEXFINDGROUPINCAND Core function to find new grains among the
+% reflections
+%
+% [goods, grain] = gtIndexFindGroupInCand(cand, isecs, tol_build,...
+%                                         samenv, cryst)    
+%
+% ------------------------------------------------------------------
+%
+% This is the function that finds possible new grains in the
+% indexing process. It searches for a grain in the input set of candidate 
 % reflections 'cand'.
-% All reflections in 'cand' are pair-consistent with the first reflection
-% in 'cand'. By looping through the rest of the set and checking group-
-% consistency, a 2nd, 3rd, 4th, etc. valid reflection is potentially found.   
+% All reflections in 'cand' are assumed to be pair-consistent with the 
+% first reflection in 'cand'. By looping through the rest of the set 
+% and checking group consistency, a 2nd, 3rd, 4th, etc. valid reflection 
+% is potentially found.   
 % 
 % INPUT 
-% 
 %  cand      - reflection data; all the Friedel pairs that possibly belong 
 %              to the same grain as the first one in cand 
 %  isecs     - intersection points of the diffraction paths in 3D space;  
 %              size (n,6), where n = length(cand)-1
-%  tol_build - indexing tolerances for "building grains" 
+%  tol_build - indexing tolerances for finding new grains (Build grains)
 %  samenv    - coordinates of the sample envelope
 %  cryst     - crystallography and other data 
 %
 % OUTPUT
-%
 %  goods     - vector of logicals; true for the reflections that 
 %              constitute the new grain; size (n,1)
-%  grain     - grain data; empty if no valid grain is found
+%  grain     - grain data; empty if no valid grain was found
 %
 
 nof_cand    = length(cand.id); 
@@ -38,14 +44,14 @@ goods       = false(nof_cand,1);
 for c1 = 2:nof_cand 
 
 	% Get data for actual reflection and the following one
-	baserefs = gtIndexSelectRefs([1,c1],cand);
+	baserefs = gtIndexSelectRefsReduced([1,c1],cand);
 
 	% Loop for finding a 3rd reflection of the grain
 	for c2 = c1+1:nof_cand 
 	
 		% Check group consistency of this third reflection c2 against the
 		% other two
-		tryref = gtIndexSelectRefs(c2,cand);
+		tryref = gtIndexSelectRefsReduced(c2,cand);
 
 		ok = gtIndexCheckGroupCons(tryref, baserefs, isecs(c1-1,:), ...
              tol_build, samenv, cryst.ACM);
diff --git a/zUtil_Indexter/gtIndexGrainOutputBasic.m b/zUtil_Indexter/gtIndexGrainOutputBasic.m
index affa8374f284e65d4c091f891c35f51967d76a4e..bed94ec77477c3068fb39395cb422138b24bf1cb 100644
--- a/zUtil_Indexter/gtIndexGrainOutputBasic.m
+++ b/zUtil_Indexter/gtIndexGrainOutputBasic.m
@@ -1,4 +1,10 @@
 function [grain, goodind] = gtIndexGrainOutputBasic(grainrefs, cryst, tol_ang)
+% GTINDEXGRAINOUTPUTBASIC Computes basic grain data from a set of
+% reflections.
+%
+% [grain, goodind] = gtIndexGrainOutputBasic(grainrefs, cryst, tol_ang)
+%
+% ---------------------------------------------------------------------
 %
 % Creates the basic 'grain' data structure from reflections 'grainrefs'. 
 % Used when a new grain is created or an existing one is modififed.
@@ -11,7 +17,6 @@ function [grain, goodind] = gtIndexGrainOutputBasic(grainrefs, cryst, tol_ang)
 % computed. 
 %
 % INPUT
-%
 %  grainrefs.   - grain reflection data
 %  cryst.       - crystallography and other data
 %  tol_ang      - max. angular orientation deviation between two grains
@@ -75,6 +80,10 @@ end
     gtCrystRodriguesTestCore(rline, refind, cryst.Rfzone_acc, ...
     cryst.Rfzone_ext, tol_Rdist);
 
+if isempty(Rvec)
+    grain = [];
+    return 
+end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Check angular deviation in real space
@@ -122,21 +131,21 @@ dang = acosd(dotprod);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 grain.id         = NaN;                      
-grain.refid      = grainrefs.id(goodind)';       % reflection id-s used in indexing
-grain.pairid     = grainrefs.pairid(goodind)';   % Friedel pair IDs from matching
+grain.refid      = grainrefs.id(goodind)';     % reflection id-s used in indexing
+grain.pairid     = grainrefs.pairid(goodind)'; % Friedel pair IDs from matching
 
 grain.R_vector   = Rvec;                     % Rodrigues vector
 grain.R_onedge   = Ronedge;                  % Rodrigues vector is close to
                                              % edge of fundamental zone
 
-grain.nof_pairs  = length(goodind);              % no. of pairs in grain
+grain.nof_pairs  = length(goodind);          % no. of pairs in grain
 
-grain.pl         = grainrefs.pl(goodind,:);      % plane normals
-grain.hkl        = grainrefs.hkl(goodind,:);     % hkl families
-grain.theta      = grainrefs.theta(goodind);     % theta
-grain.int        = grainrefs.int(goodind);       % intensity
-grain.bbxs       = grainrefs.bbxs(goodind);      % bbox x size
-grain.bbys       = grainrefs.bbys(goodind);      % bbox y size
+grain.pl         = grainrefs.pl(goodind,:);  % plane normals
+grain.hkl        = grainrefs.hkl(goodind,:); % hkl families
+grain.theta      = grainrefs.theta(goodind); % theta
+grain.int        = grainrefs.int(goodind);   % intensity
+grain.bbxs       = grainrefs.bbxs(goodind);  % bbox x size
+grain.bbys       = grainrefs.bbys(goodind);  % bbox y size
 
 grain.hklind     = hklind;            % unsigned hkl family linear index
 grain.shklind    = shklind;           % signed hkl linear index
@@ -155,14 +164,8 @@ grain.dcom = gtMathsPointLinesDists(grain.center, ...
 
 
 grain.indST   = zeros(3);  % strain tensor in indexing
-%grain.pla     = NaN;       % theoretical strained plane normal           
-%grain.pld     = NaN;       % theoretical d-spacing
 grain.normint = zeros(size(grain.int));
 
-%grain.difspotidA = [];
-%grain.difspotidB = [];
-
-
 % List of shkl-s that have been indexed
 grain.shklfound = cryst.shklfound_template;  % blank template
 
@@ -202,9 +205,4 @@ grain.stat.intstd     = std(grain.int);  % std of intensities
 grain.stat.normintstd = std(grain.int);  % std of normalised intensities
 
 
-
-
-end % of function
-
-
-
+end % of function
\ No newline at end of file
diff --git a/zUtil_Indexter/gtIndexMergeGrains.m b/zUtil_Indexter/gtIndexMergeGrains.m
index d13ec6f35d6c15700d8b3121558c3d8cf30f39fe..ff6d91a8007729dc6c0de1671caad8f51706d05c 100644
--- a/zUtil_Indexter/gtIndexMergeGrains.m
+++ b/zUtil_Indexter/gtIndexMergeGrains.m
@@ -1,17 +1,19 @@
-function [newgrains, unindexed] = gtIndexMergeGrains(grains, ...
-                                  unindexed, tot, tol_merge, cryst)
+function [newgrains, tot] = gtIndexMergeGrains(grains, tot, tol_merge, cryst)
+% GTINDEXMERGEGRAINS Merges grains that seem identical (one of the main 
+% stages of Indexter).
 %
-% FUNCTION newgrains = gtIndexMergeGrains(grains, tot, tol_merge, cryst)
+% [newgrains, tot] = gtIndexMergeGrains(grains, tot, tol_merge, cryst)
 % 
-% 2nd/3 Step of the iterative indexing strategy. It looks for  
-% combinations of exsiting grains and merges them into one grain if they 
+% --------------------------------------------------------------------
+%
+% 2nd of the 3 steps of the iterative indexing strategy. It looks for  
+% combinations of existing grains and merges them into one grain if they 
 % seem to be identical - meaning that the differences in bounding box 
 % sizes, intensities, location and orientation are within the specified 
 % tolerances.  
-% Grain parameters and statistics for merged grains are updated.
+% Updates grain parameters and statistics for merged grains.
 %
 % INPUT
-%
 %  grain        - all grain data (cell vector)
 %  tot.         - all reflection data (including indexed and unindexed)
 %  tol_merge.   - tolerances for merging grains;
@@ -20,9 +22,10 @@ function [newgrains, unindexed] = gtIndexMergeGrains(grains, ...
 %
 % OUTPUT
 %  grain        - updated grain data including merged and unmerged grains
+%  tot.indexed  - updated; true for reflections indexed in newgrains, false
+%                 for the rest
 %
 
-
 %%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Initialize grain data
 %%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -30,11 +33,17 @@ function [newgrains, unindexed] = gtIndexMergeGrains(grains, ...
 fprintf('\n MERGING GRAINS...\n\n')
 disp(tol_merge)
 disp(' ')
-tic
 
 nof_newgrains = 0;
 nof_grains    = length(grains);
 
+if (nof_grains == 0)
+    newgrains = grains;
+    return
+end
+
+tic
+
 % Grains indices to be checked
 remaining = 1:nof_grains;
 
@@ -177,7 +186,7 @@ while length(remaining) >= 1
                         nof_newgrains, num2str(newgrefs.pairid(excind)'))
                 
                 % Add excluded ones to unindexed set
-                unindexed = [unindexed; newgrefs.id(excind)];
+                tot.indexed(newgrefs.id(excind)) = false;
 
             end
             
@@ -202,9 +211,6 @@ while length(remaining) >= 1
 end
 
 
-unindexed = sort(unindexed);
-
-
 disp(' ')
 toc
 
@@ -212,12 +218,10 @@ toc
 end % of main function
 
 
-
 %%%%%%%%%%%%%%%%%%%%%%%%%%
 %% SUB-FUNCTION
 %%%%%%%%%%%%%%%%%%%%%%%%%%
 
-
 function grsel = sfSelectGrains(keeplist,gr)
 
   grsel.id        = gr.id(keeplist);
@@ -233,4 +237,3 @@ function grsel = sfSelectGrains(keeplist,gr)
   grsel.dang      = gr.dang(keeplist);
 
 end
-
diff --git a/zUtil_Indexter/gtIndexRefineGrains.m b/zUtil_Indexter/gtIndexRefineGrains.m
index b3bf6fd952e568077006699ebe0288f7489b5640..e8b3f21856ccb550b7c64824e641968e33c1032b 100644
--- a/zUtil_Indexter/gtIndexRefineGrains.m
+++ b/zUtil_Indexter/gtIndexRefineGrains.m
@@ -1,33 +1,39 @@
-function [grain, remaining] = gtIndexRefineGrains(grain, tryrefids, tot,...
-                                                  tol_b, cryst)
+function [grain, tot] = gtIndexRefineGrains(grain, tot, tol_b, cryst)
+% GTINDEXREFINEGRAINS Adds unindexed reflections to grains (one of the 
+% main stages of Indexter).
 %
-% FUNCTION [grain, remaining] = gtIndexRefineGrains(grain, tryrefids, ...
-%                                                   tot, tol_b, cryst)
+% [grain, tot] = gtIndexRefineGrains(grain, tot, tol_b, cryst)
 % 
-% 3rd/3 Step of the iterative indexing strategy. It loops through all the 
-% grains and checks if additional unindexed reflections can be added to 
-% them based on tolerances specified for building grains.  
+% ------------------------------------------------------------
+%
+% 3rd out of 3 steps of the iterative indexing strategy. It loops through
+% all the grains and checks if additional unindexed reflections can be 
+% added to them based on tolerances specified for building grains.  
 % Grain parameters and statistics for extended grains are updated, 
 % reflections are checked for multiplicity and least fitting ones are 
 % rejected if multiplicity exceeds the maximum.
 %
 % INPUT
-%
 %  grain.       - all grain data (cell vector)
-%  tryrefids    - unindexed reflection id-s in 'tot' to be tested 
 %  tot.         - all reflection data (including indexed and unindexed)
-%  tol_merge.   - tolerances for building grains;
+%  tol_b.       - tolerances for building grains;
 %                 (as in paramers.index.strategy.b.beg or b.end)
 %  cryst.       - crystallography and other data
 %
 % OUTPUT
 %  grain        - updated grain data
-%  remaining    - remaining unindexed reflection id-s
+%  tot.indexed  - updated; true for all reflection ID-s indexed in grain,
+%                 false for the rest
 %
 
 fprintf('\n REFINING GRAINS...\n\n')
 disp(tol_b)
 disp(' ')
+
+if isempty(grain)
+    return
+end
+
 tic
 
 % Grains to delete
@@ -54,7 +60,7 @@ for ii = 1:length(grain);
     grain{ii}.shkldirs = indshkldirs_cry*cry2sam'; 
     
     % load data for all candidate reflections
-    tryrefs = gtIndexSelectRefs(tryrefids, tot);
+    tryrefs = gtIndexSelectRefs(~tot.indexed, tot);
     
     % check consistency of bbox sizes, intensity, path distance 
     % and orientation
@@ -66,7 +72,7 @@ for ii = 1:length(grain);
         grain{ii}.extended = true;
         
         fprintf('Pairs added to grain #%d      :   %s\n', ii,...
-                num2str(tot.pairid(tryrefids(goodfits))'))
+                num2str(tot.pairid(tryrefs.id(goodfits))'))
     end
    
     
@@ -74,7 +80,7 @@ for ii = 1:length(grain);
     % reflection indices (if grain is found invalid, excrefids contains all
     % its refid-s)
     [grain{ii}, excrefids] = gtIndexModifyGrain(grain{ii}, [], ...
-                             tryrefids(goodfits), tot, tol_b, cryst);
+                             tryrefs.id(goodfits), tot, tol_b, cryst);
    
     
     if isempty(grain{ii})
@@ -84,13 +90,11 @@ for ii = 1:length(grain);
     
     % Delete the reflections from list that were added to the grain and add
     % back the ones that were rejected:
-    tryrefids(goodfits) = [];
-    tryrefids = [tryrefids; excrefids];
+    tot.indexed(tryrefs.id(goodfits)) = true;
+    tot.indexed(excrefids)            = false;
     
 end
 
-% Remaining reflections
-remaining = tryrefids;
 
 % Delete bad grains
 grain(gtodel) = [];
@@ -99,6 +103,4 @@ grain(gtodel) = [];
 disp(' ')
 toc
 
-end
-																			
-																			
+end % of function
\ No newline at end of file
diff --git a/zUtil_Indexter/gtIndexSelectRefs.m b/zUtil_Indexter/gtIndexSelectRefs.m
index 4d95dfdc7db6fd8d1b84f51a0c68565104b2e306..8f4219c6a3d4cd6e61550f57f25ef3b823113f33 100644
--- a/zUtil_Indexter/gtIndexSelectRefs.m
+++ b/zUtil_Indexter/gtIndexSelectRefs.m
@@ -1,6 +1,20 @@
 function refssel = gtIndexSelectRefs(keeplist, refs)
-% A helper function to load a sub-set 'keeplist' of reflection data from 
+% GTINDEXSELECTREFS Helper to load a sub-set of reflections.
+%
+% refssel = gtIndexSelectRefs(keeplist, refs)
+%
+% --------------------------------------------------
+%
+% A helper function to load the sub-set 'keeplist' of reflection data from 
 % the complete set 'refs'.
+%
+% INPUT
+%   keeplist  - linear or logical indices of reflections to be selected
+%   refs      - reflection data
+%
+% OUTPUT
+%   refssel   - selected sub-set of reflection data
+%
 
 
 refssel.id     = refs.id(keeplist);
diff --git a/zUtil_Indexter/gtIndexSelectRefsReduced.m b/zUtil_Indexter/gtIndexSelectRefsReduced.m
new file mode 100644
index 0000000000000000000000000000000000000000..148b1731aa4075fc321d6705cb9b87a82d07e524
--- /dev/null
+++ b/zUtil_Indexter/gtIndexSelectRefsReduced.m
@@ -0,0 +1,32 @@
+function refssel = gtIndexSelectRefsReduced(keeplist, refs)
+% GTINDEXSELECTREFSREDUCED Helper to load a sub-set of reflections.
+%
+% refssel = gtIndexSelectRefsReduced(keeplist, refs)
+%
+% --------------------------------------------------
+%
+% A helper function to load the sub-set 'keeplist' of reflection data from 
+% the complete set 'refs'.
+%
+% INPUT
+%   keeplist  - linear or logical indices of reflections to be selected
+%   refs      - reflection data
+%
+% OUTPUT
+%   refssel   - selected sub-set of reflection data
+%
+
+
+refssel.id     = refs.id(keeplist);
+refssel.int    = refs.int(keeplist);
+refssel.bbxs   = refs.bbxs(keeplist);
+refssel.bbys   = refs.bbys(keeplist);
+
+refssel.ca     = refs.ca(keeplist,:);
+refssel.dir    = refs.dir(keeplist,:);
+refssel.pl     = refs.pl(keeplist,:);
+
+refssel.hklind = refs.hklind(keeplist);   % cell vector
+
+
+end
\ No newline at end of file
diff --git a/zUtil_Indexter/gtIndexStatisticalFit.m b/zUtil_Indexter/gtIndexStatisticalFit.m
index cd0ca5fa3770643ded8142fe1724e3a6756361aa..d0a54e5541cc28c30a799f05f0c9fd59c6238204 100644
--- a/zUtil_Indexter/gtIndexStatisticalFit.m
+++ b/zUtil_Indexter/gtIndexStatisticalFit.m
@@ -1,28 +1,31 @@
 function [grain, remaining] = gtIndexStatisticalFit(grain, restids, tot,...
                                                    strategy_s, cryst)
+% GTINDEXSTATISTICALFIT Adds unindexed reflections to grains based on 
+% statistics.
+% 
+% [grain, remaining] = gtIndexStatisticalFit(grain, restids, ...
+%                                            tot, strategy_s, cryst)
 %
-% FUNCTION [grain, remaining] = gtIndexStatisticalFit(grain, restids, ...
-%                               tot, strategy_s, cryst)
+% ------------------------------------------------------------------
 %
-% Adds unindexed reflections to exiting grains based on statistical 
-% fitting. those reflections of which intensity, bounding box sizes, 
+% Adds unindexed reflections to existing grains based on statistical 
+% fitting. Those reflections of which intensity, bounding box sizes, 
 % distance and angles do not deviate more than the specified tolerance 
 % value (in std) from the average grain value will be added to the grain. 
 % If a reflection fits to more than one grain, it will be added to the
 % grain it fits best. 
 %
 % INPUT
-%  grain      - all grain data
-%  restids    - unindexed reflection id-s
-%  tot        - all original reflection data
-%  strategy_s - parameters for statistical fit
-%               (as in parameters.index.strategy.s)
-%  cryst      - crystallography and other data
+%   grain      - all grain data
+%   restids    - unindexed reflection id-s
+%   tot        - all original reflection data
+%   strategy_s - parameters for statistical fit
+%                (as in parameters.index.strategy.s)
+%   cryst      - crystallography and other data
 %
 % OUTPUT 
-%  grain      - all grain data updated with added reflections
-%  remaining  - remaining unindedxed reflection id-s after the fit
-%
+%   grain      - all grain data updated with added reflections
+%   remaining  - remaining unindedxed reflection ID-s after the fit
 %
 
 
@@ -120,7 +123,8 @@ remaining(todel) = [];
 % Grains to delete
 gtodel = [];
 
-parfor ii = 1:nof_grains
+% !!! ? can be a parfor loop ?
+for ii = 1:nof_grains
 	if ~isempty(addtograins{ii})
         
 		fprintf('Pairs fitted to grain #%d    :   %s\n', ii, ...
@@ -152,7 +156,6 @@ toc
 end % of main function
 
 
-
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% SUB-FUNCTIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%