diff --git a/4_grains/gtPredictPositions_3D.m b/4_grains/gtPredictPositions_3D.m
index 58300afb79789c698fc377d559fbf5078445e7d0..103c2533b4bb5b770feaa0b86b430117d24d4544 100644
--- a/4_grains/gtPredictPositions_3D.m
+++ b/4_grains/gtPredictPositions_3D.m
@@ -1,6 +1,6 @@
-function [images,u_direct,v_direct,u_diffr,v_diffr,Omega,Eta,Theta] = gtPredictPositions_3D(normal, hkl, grainCentroid, plot_flag, parameters, codeflag)
+function [images,u_direct,v_direct,u_diffr,v_diffr,Omega,Eta,Theta] = gtPredictPositions_3D(phaseid, normal, hkl, grainCentroid, plot_flag, parameters, codeflag)
 % GTPREDICTPOSITIONS_3D  predict positions of spots in images
-%     [images,u_direct,v_direct,u_diffr,v_diffr,Omega,Eta,Theta] = gtPredictPositions_3D(normal, hkl, grainCentroid, plot_flag, parameters, codeflag)
+%     [images,u_direct,v_direct,u_diffr,v_diffr,Omega,Eta,Theta] = gtPredictPositions_3D(phaseid, normal, hkl, grainCentroid, plot_flag, parameters, codeflag)
 %     -----------------------------------------------------------------------------------------------------------------------------------------------
 %
 %     if parameters.geo includes the neccessary fields for the arbitary
@@ -12,6 +12,9 @@ function [images,u_direct,v_direct,u_diffr,v_diffr,Omega,Eta,Theta] = gtPredictP
 %     positions of the ext and dif spots.
 %     plot_flag=1 plots images
 %
+%     Version 005 01-08-2013 by LNervo
+%       Replaced gtCalculateDist with gtCrystDSpacing and gtCrystTheta
+%
 %     Modified on March 08 and March 14, 2012 by PReischig
 %       Adapted input for gtCalculateDist.
 %    
@@ -27,7 +30,7 @@ function [images,u_direct,v_direct,u_diffr,v_diffr,Omega,Eta,Theta] = gtPredictP
 %       Works with old coordinate system 
 
 acq=parameters.acq;
-geo=parameters.geo;
+geo=parameters.geo;%obsolete...it does not exist anymore
 labgeo = parameters.labgeo;
 
 % do we have all the required fields?
@@ -55,7 +58,7 @@ if all(isfield(geo, FEDgeofields)) && codeflag
     omstep    = 180/acq.nproj;
 
     % given a hkl plane, calculate the interplanar distance.
-    [~,Theta] = gtCalculateDist(hkl,parameters.cryst,parameters.acq.energy);
+    Theta = gtCrystTheta( gtCrystDSpacing(hkl, gtCrystHKL2CartesianMatrix(parameters.cryst(phaseid).latticepar) ), parameters.acq.energy);
     sinth = sind(Theta);
 
     % calc Omega angles
diff --git a/zUtil_Cryst/gtCrystTheta.m b/zUtil_Cryst/gtCrystTheta.m
new file mode 100644
index 0000000000000000000000000000000000000000..6aa4c89011b986ca40b4cb91915b6416e0e11757
--- /dev/null
+++ b/zUtil_Cryst/gtCrystTheta.m
@@ -0,0 +1,28 @@
+function theta = gtCrystTheta(dsp,energy)
+% GTCRYSTTHETA  Calculates theta angle given a list of d-spacing
+%     theta = gtCrystTheta(dsp,energy)
+%     --------------------------------
+%    
+%     INPUT:
+%       dsp    = d-spacing (Angstron)
+%       energy = beam energy in keV 
+% 
+%     OUTPUT:
+%       theta = the Bragg angle of the specified hkl family by given the
+%               d-spacing list
+%
+%
+%     Version 001 10-07-2013 by LNervo
+
+
+if ~exist('energy','var') || isempty(energy)
+    load parameters
+    energy = parameters.acq.energy;
+end
+
+lambda = gtConvEnergyToWavelength(energy);
+
+% theta from the Bragg's law
+theta = asind(lambda./(2*dsp));
+
+end % end of function
diff --git a/zUtil_Cryst/gtDetectorTwotheta.m b/zUtil_Cryst/gtDetectorTwotheta.m
index 0acfd076a426d5ac9290e8652ed5327789cca75e..4bb697bf10a9222090ba32c0c768da0f0718b775 100644
--- a/zUtil_Cryst/gtDetectorTwotheta.m
+++ b/zUtil_Cryst/gtDetectorTwotheta.m
@@ -1,5 +1,5 @@
 function results=gtDetectorTwotheta(parameters, phaseID)
-% GTCALCULATETWOTHETA  Returns reflections and theta values on detector
+% GTDETECTORTWOTHETA  Returns reflections and theta values on detector
 %     results=gtDetectorTwotheta(parameters, phaseID)
 %     ----------------------------------------------------------
 %     INPUT:
@@ -56,13 +56,10 @@ disp(['hkl total families:       ' num2str(size(hkl, 1))])
 
 % ensure that the angles we have are correct for the energy and lattice parameters 
 % in the parameters file
-for j = 1:size(hkl, 1)
-    % calculate d-spacing and twotheta for each reflection
-    [d0(j), tt(j)] = gtCalculateDist(hkl(j, :), parameters.cryst(phaseID), parameters.acq.energy);
-end
+Bmat = gtCrystHKL2CartesianMatrix(parameters.cryst(phaseID).latticepar);
+d0 = gtCrystDSpacing(hkl', Bmat); % hkl as column vector
+tt = gtCrystTheta(d0, parameters.acq.energy);
 
-d0 = d0'; % fix bug 22/08/2012 LNervo
-tt = tt';
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Keep only reflections on the detector
@@ -71,8 +68,8 @@ minangle = parameters.labgeo.detanglemin;
 maxangle = parameters.labgeo.detanglemax;
 
 % remove reflections out of detector
-ind = find(tt > maxangle);   % ind contains reflections whichdo not impinge on the detector
-ind = [ind; find(tt < minangle)];
+ind = find(tt > maxangle/2);   % ind contains reflections whichdo not impinge on the detector
+ind = [ind; find(tt < minangle/2)];
 % ...remove them
 tt(ind)    = [];            
 hkl(ind, :) = [];
diff --git a/zUtil_Cryst/gtSymmetricReflections.m b/zUtil_Cryst/gtSymmetricReflections.m
index b61f334fe9d6c4677aead561e70fae8ceb528ad8..f54b305cd40f8064399b240bc40a22ec7dbc2508 100644
--- a/zUtil_Cryst/gtSymmetricReflections.m
+++ b/zUtil_Cryst/gtSymmetricReflections.m
@@ -22,6 +22,9 @@ function results = gtSymmetricReflections(cryst,energy)
 %              .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
@@ -67,10 +70,10 @@ for ii = 1:size(hkltypes_used,1)
     mult      = [mult; nhkls];
 end
 
-for jj = 1:size(hklsp,1)
-    % calculate d-spacing and twotheta for each reflection
-    [dsp(jj), theta(jj)] = gtCalculateDist(hklsp(jj,:), cryst, energy);
-end
+Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
+dsp = gtCrystDSpacing(hklsp', Bmat);
+theta = gtCrystTheta(dsp, energy);
+
 
 % simpler output (should ditch the hashtable)
 results.hklsp       = int16(hklsp);