diff --git a/3_pairmatchingGUI/gtMatchSaveResults.m b/3_pairmatchingGUI/gtMatchSaveResults.m
index 2799ff015e8a8e8480938634e119f8827d53d061..9b0818c950ce65f0bb96663f2371ddacfe43c95d 100644
--- a/3_pairmatchingGUI/gtMatchSaveResults.m
+++ b/3_pairmatchingGUI/gtMatchSaveResults.m
@@ -108,12 +108,12 @@ scattveclab = gtGeoScattVecFromDiffVec(handles.pairsfilt.diffveclab, ...
                                        handles.labgeo.beamdir);
 
 handles.pairsfilt.planenorm = gtGeoLab2Sam(scattveclab,...
-    handles.pairsfilt.omegaA, handles.labgeo, samgeo, 1, 1);
+    handles.pairsfilt.omegaA, handles.labgeo, samgeo, true, 'only_rot', true);
 
 
 % Diffraction vector
 handles.pairsfilt.diffvecsam = gtGeoLab2Sam(handles.pairsfilt.diffveclab,...
-    handles.pairsfilt.omegaA, handles.labgeo, samgeo, 1, 1);
+    handles.pairsfilt.omegaA, handles.labgeo, samgeo, true, 'only_rot', true);
 
 
 % Difspot A in Sample reference
@@ -123,7 +123,7 @@ centAlab = gtGeoDet2Lab([handles.pairsfilt.centAU,handles.pairsfilt.centAV],...
            handles.labgeo, 0) - handles.pairsfilt.shiftA;
 
 handles.pairsfilt.centsamA = gtGeoLab2Sam(centAlab, handles.pairsfilt.omegaA,...
-                                          handles.labgeo, samgeo, 0, 0);
+                                          handles.labgeo, samgeo, false);
 
 % Spot B data is not important.
 %handles.pairsfilt.centsamB = zeros(length(handles.pairsfilt.indA),3);
diff --git a/7_fed2/gtFedFwdProjExact.m b/7_fed2/gtFedFwdProjExact.m
index 066872fdc0164f11a7dbbc94f0475768e7dc7296..4ae88f7cb90dcc59943659638e42901acbb39442 100755
--- a/7_fed2/gtFedFwdProjExact.m
+++ b/7_fed2/gtFedFwdProjExact.m
@@ -31,7 +31,7 @@ function [gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars, parameters, det_num
     % beam direction and rotation axis are defined.
     % Use the Sample -> Lab orientation transformation assuming omega=0;
     % (vector length preserved for free vectors)
-    plorig = gtGeoSam2Lab(plorig, eye(3), labgeo, samgeo, true, false);
+    plorig = gtGeoSam2Lab(plorig, eye(3), labgeo, samgeo, true);
 
     if (isfield(fedpars, 'detector') ...
             && isfield(fedpars.detector, 'psf') ...
@@ -254,7 +254,7 @@ function [uvw_bl, rot_s2l] = predict_uvw(plsam, gvcs, defT, bl, labgeo, samgeo,
         sinth, labgeo.beamdir', labgeo.rotdir', labgeo.rotcomp, bl.omind);
 
     % Delete those where no reflection occurs
-    if any(isnan(om));
+    if (any(isnan(om)))
         inds_bad = find(isnan(om));
         gvd(:, inds_bad(1:10))
         error('gtFedFwdProjExact:bad_R_vectors', ...
diff --git a/zUtil_Deformation/gtDefFwdProjGvdm2NW.m b/zUtil_Deformation/gtDefFwdProjGvdm2NW.m
index 1b6735b0a343af63c439fb1a632174bcb01f7661..8e55d979a40f25fe9618524eca2753903c877317 100644
--- a/zUtil_Deformation/gtDefFwdProjGvdm2NW.m
+++ b/zUtil_Deformation/gtDefFwdProjGvdm2NW.m
@@ -67,7 +67,7 @@ function bl = gtDefFwdProjGvdm2NW(grain, ref_sel, gv, fedpars, parameters, eta_s
     % beam direction and rotation axis are defined.
     % Use the Sample -> Lab orientation transformation assuming omega=0;
     % (vector length preserved for free vectors)
-    pls_orig = gtGeoSam2Lab(pls_orig, eye(3), labgeo, samgeo, true, false);
+    pls_orig = gtGeoSam2Lab(pls_orig, eye(3), labgeo, samgeo, true);
 
     linear_interp = ~(isfield(fedpars, 'projector') ...
         && strcmpi(fedpars.projector, 'nearest'));
diff --git a/zUtil_Deformation/gtDefFwdProjGvdm2UVW.m b/zUtil_Deformation/gtDefFwdProjGvdm2UVW.m
index 826e1d1ed60d98730d8f1d76faf079919ed7edd7..81f867aa424cdf437dfb1ce35893843f3f077d25 100755
--- a/zUtil_Deformation/gtDefFwdProjGvdm2UVW.m
+++ b/zUtil_Deformation/gtDefFwdProjGvdm2UVW.m
@@ -43,7 +43,7 @@ function bl = gtDefFwdProjGvdm2UVW(grain, ref_sel, gv, fedpars, parameters, det_
     % beam direction and rotation axis are defined.
     % Use the Sample -> Lab orientation transformation assuming omega=0;
     % (vector length preserved for free vectors)
-    pls_orig = gtGeoSam2Lab(pls_orig, eye(3), labgeo, samgeo, true, false);
+    pls_orig = gtGeoSam2Lab(pls_orig, eye(3), labgeo, samgeo, true);
 
     if (isfield(fedpars, 'detector') ...
             && isfield(fedpars.detector, 'psf') ...
@@ -137,7 +137,7 @@ function bl = gtDefFwdProjGvdm2UVW(grain, ref_sel, gv, fedpars, parameters, det_
             dvec_lab = gtFedPredictDiffVecMultiple(pllab, labgeo.beamdir');
 
             rot_s2l = permute(rot_l2s, [2 1 3]);
-            gvcs_lab = gtGeoSam2Lab(gvpcs', rot_s2l, labgeo, samgeo, false, [], element_wise_gcs);
+            gvcs_lab = gtGeoSam2Lab(gvpcs', rot_s2l, labgeo, samgeo, false, 'element_wise', element_wise_gcs);
 
             uvw_bl = gtFedPredictUVWMultiple([], dvec_lab, gvcs_lab', ...
                 detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
diff --git a/zUtil_Deformation/gtDefFwdProjGvdm2W.m b/zUtil_Deformation/gtDefFwdProjGvdm2W.m
index fbe3a62957f0d6730be25f63a2ed7fa7b4a9d427..86e763c64ed54d2a748938bd6a89ccfe3453077b 100644
--- a/zUtil_Deformation/gtDefFwdProjGvdm2W.m
+++ b/zUtil_Deformation/gtDefFwdProjGvdm2W.m
@@ -60,7 +60,7 @@ function bl = gtDefFwdProjGvdm2W(grain, ref_sel, gv, fedpars, parameters, det_in
     % beam direction and rotation axis are defined.
     % Use the Sample -> Lab orientation transformation assuming omega=0;
     % (vector length preserved for free vectors)
-    pls_orig = gtGeoSam2Lab(pls_orig, eye(3), labgeo, samgeo, true, false);
+    pls_orig = gtGeoSam2Lab(pls_orig, eye(3), labgeo, samgeo, true);
 
     linear_interp = ~(isfield(fedpars, 'projector') ...
         && strcmpi(fedpars.projector, 'nearest'));
diff --git a/zUtil_Deformation/gtDefShapeFunctionsCreateNW.m b/zUtil_Deformation/gtDefShapeFunctionsCreateNW.m
index cda6a88b957ef7251f24becaf297132eaf00ab97..87a8ec31f213233869c9e8b4fc1c5faaa28cb957 100644
--- a/zUtil_Deformation/gtDefShapeFunctionsCreateNW.m
+++ b/zUtil_Deformation/gtDefShapeFunctionsCreateNW.m
@@ -45,7 +45,7 @@ function sf = gtDefShapeFunctionsCreateNW(sampler, eta_resoltion, varargin)
     g = gtMathsRod2OriMat(ref_gr.R_vector');
     pl_cry = gtVectorLab2Cryst(pl_samd, g)';
 
-%             pl_labd = gtGeoSam2Lab(pl_samd, eye(3), self.parameters.labgeo, self.parameters.samgeo, true, false)';
+%     pl_labd = gtGeoSam2Lab(pl_samd, eye(3), self.parameters.labgeo, self.parameters.samgeo, true)';
 
     fprintf('Computing bounds of NW shape functions: ')
     c = tic();
diff --git a/zUtil_Deformation/gtDefShapeFunctionsCreateW.m b/zUtil_Deformation/gtDefShapeFunctionsCreateW.m
index dad67c497d54511f0f1ce747641e66e6a09f3690..0521e7118057476c1eb331711d8d2f3ba294f9ce 100644
--- a/zUtil_Deformation/gtDefShapeFunctionsCreateW.m
+++ b/zUtil_Deformation/gtDefShapeFunctionsCreateW.m
@@ -20,7 +20,7 @@ function sf = gtDefShapeFunctionsCreateW(sampler)
         g = gtMathsRod2OriMat(ref_gr.R_vector');
         pl_cry = gtVectorLab2Cryst(pl_samd, g)';
     end
-%     pl_labd = gtGeoSam2Lab(pl_samd, eye(3), self.parameters.labgeo, self.parameters.samgeo, true, false)';
+%     pl_labd = gtGeoSam2Lab(pl_samd, eye(3), self.parameters.labgeo, self.parameters.samgeo, true)';
 
     sub_space_size = sampler.stats.sampling.gaps;
     ospace_samp_size = sampler.get_orientation_sampling_size();
diff --git a/zUtil_Deformation/gtDefShapeFunctionsNW2UVW.m b/zUtil_Deformation/gtDefShapeFunctionsNW2UVW.m
index bc94aa3ffe81bb99a563928546ca4a8a080f018b..fef1813f68b6b94029e96bd380f56bca0f934cd9 100644
--- a/zUtil_Deformation/gtDefShapeFunctionsNW2UVW.m
+++ b/zUtil_Deformation/gtDefShapeFunctionsNW2UVW.m
@@ -97,7 +97,7 @@ function sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, varargin)
         % Calculating all the sample center positions at each omega inside
         % the blobs
         rot_s2l_all = gtMathsRotationTensor(all_ws, rotcomp_w);
-        gclab_v_all = gtGeoSam2Lab(centers, rot_s2l_all, labgeo, samgeo, false, [], false);
+        gclab_v_all = gtGeoSam2Lab(centers, rot_s2l_all, labgeo, samgeo, false, 'element_wise', false);
         gclab_v_all = reshape(gclab_v_all', 3, 1, tot_ws, []);
 
         %%% Here we put the scattering vectors and center positions
diff --git a/zUtil_Fit/gtFitFwdSimGrains.m b/zUtil_Fit/gtFitFwdSimGrains.m
index 72181795cd2eff3822d985f4f3e36fb68ea72d88..95db81c3c266a61db00686a50ddb6175a95e5360 100644
--- a/zUtil_Fit/gtFitFwdSimGrains.m
+++ b/zUtil_Fit/gtFitFwdSimGrains.m
@@ -213,10 +213,10 @@ for ii = 1:length(grain)
         [plsamd4, gcsamd4, sinthd4] = gtFitApplyDrifts(plsamd, gcsam, dspd, om4, drift, energy);
         
         % Plane normals in Lab at omega positions after drifts
-        pllab1 = gtGeoSam2Lab(plsamd1', om1, labgeo, samgeo, 1, 1);
-        pllab2 = gtGeoSam2Lab(plsamd2', om2, labgeo, samgeo, 1, 1);
-        pllab3 = gtGeoSam2Lab(plsamd3', om3, labgeo, samgeo, 1, 1);
-        pllab4 = gtGeoSam2Lab(plsamd4', om4, labgeo, samgeo, 1, 1);
+        pllab1 = gtGeoSam2Lab(plsamd1', om1, labgeo, samgeo, true, 'only_rot', true);
+        pllab2 = gtGeoSam2Lab(plsamd2', om2, labgeo, samgeo, true, 'only_rot', true);
+        pllab3 = gtGeoSam2Lab(plsamd3', om3, labgeo, samgeo, true, 'only_rot', true);
+        pllab4 = gtGeoSam2Lab(plsamd4', om4, labgeo, samgeo, true, 'only_rot', true);
         pllab  = [pllab1; pllab2; pllab3; pllab4]';
     
         %pllab01 = sam2lab0 * plsamd1;
diff --git a/zUtil_Fit/gtFitSpotPosDeviations.m b/zUtil_Fit/gtFitSpotPosDeviations.m
index 25cd1617e3e1e609b354cd9dd546e8ec33d1a407..b16ccbd45a02c58aec06c21a022369597eb4db04 100644
--- a/zUtil_Fit/gtFitSpotPosDeviations.m
+++ b/zUtil_Fit/gtFitSpotPosDeviations.m
@@ -194,7 +194,7 @@ end
 
 
 % Plane normal in Lab reference at omega
-sp.pl = gtGeoSam2Lab(sp.pl', sp.om, labgeo, par.samgeo, 1, 1)';
+sp.pl = gtGeoSam2Lab(sp.pl', sp.om, labgeo, par.samgeo, ture, 'only_rot', true)';
 
 % From degrees to images:
 sp.om = sp.om/par.omstep;
@@ -276,8 +276,8 @@ error('Function is under development...' )
     % The plane normals and grain centers need to be brought in the Lab 
     % reference where the beam direction and rotation axis are defined.
     %(Should be done the other way, in the Sample ref., as that's faster)
-    pllab = gtGeoSam2Lab(pld', om, labgeo, samgeo, 1, 1)';
-    gclab = gtGeoSam2Lab(gcd', om, labgeo, samgeo, 0, 0)';
+    pllab = gtGeoSam2Lab(pld', om, labgeo, samgeo, true, 'only_rot', true)';
+    gclab = gtGeoSam2Lab(gcd', om, labgeo, samgeo, false)';
    
     % Check if all reflections occur
     %if any(isnan(omp))
@@ -286,7 +286,7 @@ error('Function is under development...' )
         
     % Diffraction vector
     dveclab = gtFedPredictDiffVecMultiple(pllab, labgeo.beamdir');
-    %dvecsam = gtGeoLab2Sam(dveclab', omp, labgeo, samgeo, 1, 1)';
+    %dvecsam = gtGeoLab2Sam(dveclab', omp, labgeo, samgeo, true, 'only_rot', true)';
 
     % Absolute detector coordinates U,V in pixels 
     uv = gtFedPredictUVMultiple([], dveclab, gclab, labgeo.detrefpos', ...
diff --git a/zUtil_ForwardSim/gtForwardSimulate_v2.m b/zUtil_ForwardSim/gtForwardSimulate_v2.m
index 0de4c59831d1c25007268f6dadf6cbad86f1266d..b15023b66c195c4f9968a8c0d3c36b7a721391f1 100644
--- a/zUtil_ForwardSim/gtForwardSimulate_v2.m
+++ b/zUtil_ForwardSim/gtForwardSimulate_v2.m
@@ -851,7 +851,7 @@ function proj_bls = predict_spot_masks_from_indexed(gr, grain_INDX, ...
         segmentedSpots(det_index).CentroidX(difspotID_INDXTR), ...
         segmentedSpots(det_index).CentroidY(difspotID_INDXTR)];
     projvecs = gtGeoDet2Lab(projvecs, detgeo, false);
-    projvecs = gtGeoLab2Sam(projvecs, Ws, labgeo, samgeo, 0, 0) - gr_center;
+    projvecs = gtGeoLab2Sam(projvecs, Ws, labgeo, samgeo, false) - gr_center;
 
     verts_rec = gtFwdSimComputeCircumscribingPolyhedron(gr.center, ...
         projvecs, Ws, BBs, parameters, det_index, 'convhull');
diff --git a/zUtil_ForwardSim/gtFwdSimComputeCircumscribingPolyhedron.m b/zUtil_ForwardSim/gtFwdSimComputeCircumscribingPolyhedron.m
index e6945b228f3f27f341de42810a6523866a70a01a..5229bac3b028bd5a1114b3a790380d24a523b7c7 100644
--- a/zUtil_ForwardSim/gtFwdSimComputeCircumscribingPolyhedron.m
+++ b/zUtil_ForwardSim/gtFwdSimComputeCircumscribingPolyhedron.m
@@ -45,10 +45,10 @@ function verts = gtFwdSimComputeCircumscribingPolyhedron(gr_center, projvecs, om
     cbr_pos_lab = gtGeoDet2Lab([(bb(:, 1) + bb(:, 3)), (bb(:, 2) + bb(:, 4))], detgeo, false);
 
     % positions in the reconstruction coordinates
-    ctl_pos_rec = gtGeoLab2Sam(ctl_pos_lab, rot_s2l_w, labgeo, recgeo, 0, 0);
-    ctr_pos_rec = gtGeoLab2Sam(ctr_pos_lab, rot_s2l_w, labgeo, recgeo, 0, 0);
-    cbl_pos_rec = gtGeoLab2Sam(cbl_pos_lab, rot_s2l_w, labgeo, recgeo, 0, 0);
-    cbr_pos_rec = gtGeoLab2Sam(cbr_pos_lab, rot_s2l_w, labgeo, recgeo, 0, 0);
+    ctl_pos_rec = gtGeoLab2Sam(ctl_pos_lab, rot_s2l_w, labgeo, recgeo, false);
+    ctr_pos_rec = gtGeoLab2Sam(ctr_pos_lab, rot_s2l_w, labgeo, recgeo, false);
+    cbl_pos_rec = gtGeoLab2Sam(cbl_pos_lab, rot_s2l_w, labgeo, recgeo, false);
+    cbr_pos_rec = gtGeoLab2Sam(cbr_pos_lab, rot_s2l_w, labgeo, recgeo, false);
 
     projvecs_rec = gtGeoSam2Sam(projvecs, samgeo, recgeo, true);
     projvecs_rec = gtMathsNormalizeVectorsList(projvecs_rec);
diff --git a/zUtil_Geo/gtGeoDiffVecInSample.m b/zUtil_Geo/gtGeoDiffVecInSample.m
index da1bdd4dc55e34941347c4e9777ca5d8a4c18e0c..ddf96b67a7f77ab6867240db3cc73f508db508bb 100644
--- a/zUtil_Geo/gtGeoDiffVecInSample.m
+++ b/zUtil_Geo/gtGeoDiffVecInSample.m
@@ -48,7 +48,7 @@ function [diffvec, samAXYZ, theta] = gtGeoDiffVecInSample(detAUV, omega, samBXYZ
     labA = gtGeoDet2Lab(detAUV, detgeo, 0);
 
     % Coordinate transformation from Lab to Sample
-    samAXYZ = gtGeoLab2Sam(labA, omega, labgeo, samgeo, 0, 0);
+    samAXYZ = gtGeoLab2Sam(labA, omega, labgeo, samgeo, false);
 
     % Diffraction vectors from points samXYZ to point samA:
     diffvec = samAXYZ - samBXYZ(ones(size(samAXYZ, 1), 1), :);
@@ -65,7 +65,7 @@ function [diffvec, samAXYZ, theta] = gtGeoDiffVecInSample(detAUV, omega, samBXYZ
         % Coordinate transformation from Lab to Sample for beam direction
         % (free vector)
         beamdirlab = labgeo.beamdir(ones(length(omega), 1), :);
-        beamdirsam = gtGeoLab2Sam(beamdirlab, omega, labgeo, samgeo, 1, 1);
+        beamdirsam = gtGeoLab2Sam(beamdirlab, omega, labgeo, samgeo, true, 'only_rot', true);
 
         % Use the dot products of beamdir and diffvec to get theta
         theta = 0.5 * acosd( sum(beamdirsam.*diffvec, 2) );
diff --git a/zUtil_Geo/gtGeoDiffVecInSamplePair.m b/zUtil_Geo/gtGeoDiffVecInSamplePair.m
index 7e242210af6691c146ebed83376e4f900ec1bdd6..a367e5ae6b7c457263ac313fbd2ea9b0e82d1e27 100644
--- a/zUtil_Geo/gtGeoDiffVecInSamplePair.m
+++ b/zUtil_Geo/gtGeoDiffVecInSamplePair.m
@@ -73,11 +73,11 @@ function [diffvec, samA, theta] = gtGeoDiffVecInSamplePair(detUVA, detUVB, ...
 
     % Coordinate transformation from Lab to Sample (free vector)
     omega   = (omegaA + omegaB) / 2;
-    diffvec = gtGeoLab2Sam(dveclab, omega, labgeo, samgeo, 1, 1);
+    diffvec = gtGeoLab2Sam(dveclab, omega, labgeo, samgeo, true, 'only_rot', true);
 
     if (nargout > 1)
         % Coordinate transformation from Lab to Sample
-        samA = gtGeoLab2Sam(labA, omegaA, labgeo, samgeo, 0, 0);
+        samA = gtGeoLab2Sam(labA, omegaA, labgeo, samgeo, false);
 
         if (nargout == 3)
             % Use the dot products of beamdir and diffvec to get theta
diff --git a/zUtil_Geo/gtGeoDiffractometerAlignPllab.m b/zUtil_Geo/gtGeoDiffractometerAlignPllab.m
new file mode 100644
index 0000000000000000000000000000000000000000..67eb4f377d0eb0c14672dbc2dbdb9074caf72c1a
--- /dev/null
+++ b/zUtil_Geo/gtGeoDiffractometerAlignPllab.m
@@ -0,0 +1,24 @@
+function [tilts, feasible] = gtGeoDiffractometerAlignPllab(diff, gr, pl_ind, lims)
+    num_sample_tilts = size(diff.axes_sam_tilt, 1);
+    if (~exist('pl_ind', 'var') || isempty(pl_ind))
+        pl_ind = 1:numel(gr.allblobs.omega);
+    end
+    pl_sam = permute(gr.allblobs.pl(pl_ind, :), [2 3 1]);
+    pl_sam(4, 1, :) = 1;
+
+    if (~exist('lims', 'var') || isempty(lims))
+        lims = [10, 20];
+    end
+
+    rot_ax = find(diff.axes_rotation);
+    other_axes = gtMathsCross(diff.axes_sam_tilt, diff.axes_rotation);
+    [ind_axes, ~] = find(other_axes);
+    for ii = num_sample_tilts:-1:1
+        tilts(ii, :) = atan2d(pl_sam(ind_axes(ii), 1, :), pl_sam(rot_ax, 1, :)); %#ok<FNDSB>
+        t = gtMathsRotationTranslationTensor(diff.axes_sam_tilt(ii, :), tilts(ii, :), diff.origin_sam_tilt, true);
+        pl_sam = gtMathsMatrixProduct(t, pl_sam);
+    end
+
+    feasible = all(bsxfun(@le, abs(tilts), lims'), 1);
+end
+
diff --git a/zUtil_Geo/gtGeoDiffractometerDefinition.m b/zUtil_Geo/gtGeoDiffractometerDefinition.m
new file mode 100644
index 0000000000000000000000000000000000000000..01a72ded64bf4106cfd65029b838e5721200910e
--- /dev/null
+++ b/zUtil_Geo/gtGeoDiffractometerDefinition.m
@@ -0,0 +1,27 @@
+function diff = gtGeoDiffractometerDefinition(types)
+    diff = struct( ...
+        'axes_basetilt', {zeros(1, 3)}, ...
+        'axes_rotation', {zeros(1, 3)}, ...
+        'axes_sam_tilt', {zeros(0, 3)}, ... % They are applied in reverse order (last->first) for SAM->LAB transformations in gtGeoDiffractometerTensor
+        'origin_basetilt', {zeros(1, 3)}, ...
+        'origin_rotation', {zeros(1, 3)}, ...
+        'origin_sam_tilt', {zeros(0, 3)} );
+
+    if (iscell(types))
+        diff(2:numel(types)) = diff;
+        for ii = 1:numel(types)
+            diff(ii) = gtAcqDiffractometerDefinition(types{ii});
+        end
+    else
+        switch (lower(types))
+            case 'esrf_id11'
+                diff.axes_basetilt = [0, 1, 0];
+                diff.axes_rotation = [0, 0, 1];
+                diff.axes_sam_tilt = [0, 1, 0; 1, 0, 0];
+                diff.origin_sam_tilt = [0, 0, 0; 0, 0, 0];
+            otherwise
+                error('gtAcqDiffractometerDefinition:wrong_argument', ...
+                    'Unknown diffractometer: %s', types)
+        end
+    end
+end
diff --git a/zUtil_Geo/gtGeoDiffractometerTensor.m b/zUtil_Geo/gtGeoDiffractometerTensor.m
new file mode 100644
index 0000000000000000000000000000000000000000..729dfc468752e3e94548d2ca287c0b8bb7935144
--- /dev/null
+++ b/zUtil_Geo/gtGeoDiffractometerTensor.m
@@ -0,0 +1,69 @@
+function [t, t_sam, t_rot, t_base] = gtGeoDiffractometerTensor(diff, type, freevec, varargin)
+    conf = struct( ...
+        'angles_basetilt', 0, ...
+        'angles_rotation', 0, ...
+        'angles_sam_tilt', zeros(0, 1) );
+    conf = parse_pv_pairs(conf, varargin);
+
+    num_rotations = numel(conf.angles_rotation);
+    num_basetilts = numel(conf.angles_basetilt);
+    if (num_basetilts > 1 && num_rotations > 1 && num_basetilts ~= num_rotations)
+        error('gtGeoDiffractometerTensor:wrong_argument', ...
+            'Only one direction allowed for scanning at a time')
+    end
+
+    % Whether we do a SAM->LAB or a LAB->SAM tensor
+    is_sam2lab = strcmpi(type, 'sam2lab');
+
+    % Building sample tilt components
+    if (iscell(conf.angles_sam_tilt))
+        error('gtGeoDiffractometerTensor:not_implemented', 'multiple_sample tilts not implemented yet')
+    elseif (isempty(conf.angles_sam_tilt) || isempty(diff.axes_sam_tilt))
+        t_sam = eye(4);
+    else
+        t_sam = eye(4);
+        if (is_sam2lab)
+            axes_order = size(diff.axes_sam_tilt, 1):-1:1;
+        else
+            axes_order = 1:size(diff.axes_sam_tilt, 1);
+        end
+
+        for ii = axes_order
+            sam_tilt_ii = build_tensor(diff.axes_sam_tilt(ii, :), ...
+                conf.angles_sam_tilt(ii), diff.origin_sam_tilt(ii, :), ...
+                is_sam2lab, freevec);
+            t_sam = sam_tilt_ii * t_sam;
+        end
+    end
+
+    % Building rotation components
+    if (isempty(conf.angles_rotation))
+        t_rot = eye(4);
+    else
+        t_rot = build_tensor(diff.axes_rotation, conf.angles_rotation, ...
+            diff.origin_rotation, is_sam2lab, freevec);
+    end
+
+    % Building base tilt components
+    if (isempty(conf.angles_basetilt))
+        t_base = eye(4);
+    else
+        t_base = build_tensor(diff.axes_basetilt, conf.angles_basetilt, ...
+            diff.origin_basetilt, is_sam2lab, freevec);
+    end
+
+    if (is_sam2lab)
+        t = gtMathsMatrixProduct(t_base, gtMathsMatrixProduct(t_rot, t_sam));
+    else
+        t = gtMathsMatrixProduct(t_sam, gtMathsMatrixProduct(t_rot, t_base));
+    end
+end
+
+function t = build_tensor(t_axis, t_angles, t_orig, is_sam2lab, freevec)
+    if (is_sam2lab)
+        t = gtMathsRotationTranslationTensor(t_axis, t_angles, t_orig, freevec);
+    else
+        t = gtMathsRotationTranslationTensor(t_axis, -t_angles, t_orig, freevec);
+    end
+end
+
diff --git a/zUtil_Geo/gtGeoLab2Sam.m b/zUtil_Geo/gtGeoLab2Sam.m
index 3b82fa6fd16b2629a413b1980586f4631b31f669..e238754ff3198e413ca2d7ad2e6c3ee1d133d6f6 100644
--- a/zUtil_Geo/gtGeoLab2Sam.m
+++ b/zUtil_Geo/gtGeoLab2Sam.m
@@ -1,9 +1,10 @@
-function [samXYZ, lab2sam] = gtGeoLab2Sam(labXYZ, omega, labgeo, samgeo, ...
-                             freevec, rot_flag, element_wise)
+function [samXYZ, lab2sam, rottensors] = gtGeoLab2Sam(labXYZ, omega, labgeo, samgeo, freevec, varargin)
 % GTGEOLAB2SAM  LAB -> SAMPLE reference coordinate transformation
 %
-% [samXYZ, lab2sam] = gtGeoLab2Sam(labXYZ, omega, labgeo, ...
-%                                  samgeo, freevec, rot_flag)
+% [samXYZ, lab2sam] = gtGeoLab2Sam(labXYZ, omega, labgeo, samgeo, varargin)
+%
+% conf = struct('only_rot', false, 'element_wise', true, ...
+%     'diffractometer', [], 'phi', 0, 'omega', 0, 'alpha', 0, 'beta', 0);
 %
 % ----------------------------------------------------------------------
 %
@@ -21,7 +22,8 @@ function [samXYZ, lab2sam] = gtGeoLab2Sam(labXYZ, omega, labgeo, samgeo, ...
 %   samgeo.voxsize - [x,y,z] voxel sizes in lab units
 %   freevec        - 0 for fixed vectors (e.g. location);
 %                    1 for free vectors (e.g. orientation/direction)
-%   rot_flag       - if true, only applies rotation so the norm of the 
+% (Optional)
+%   only_rot       - if true, only applies rotation so the norm of the 
 %                    input vectors will be unchanged; 
 %                    i.e. no scaling applied according to the 
 %                    voxel sizes (only works if 'samgeo.lab2sam' is undefined)
@@ -36,13 +38,16 @@ function [samXYZ, lab2sam] = gtGeoLab2Sam(labXYZ, omega, labgeo, samgeo, ...
 %                    lab vectors.
 %
 
-if (~exist('rot_flag', 'var') || isempty(rot_flag))
-    rot_flag = false;
-end
+conf = struct('only_rot', false, 'element_wise', true, ...
+    'diffractometer', [], 'phi', 0, 'omega', 0, 'alpha', 0, 'beta', 0);
+conf = parse_pv_pairs(conf, varargin);
 
-if (~exist('element_wise', 'var') || isempty(element_wise))
-    element_wise = true;
-end
+[size_omega(1), size_omega(2), size_omega(3)] = size(omega);
+
+use_instr_tensor_4x4 = all(size_omega(1:2) == 4);
+use_instr_tensor_3x3 = all(size_omega(1:2) == 3);
+
+use_diffractometer = ~isempty(conf.diffractometer) || use_instr_tensor_4x4;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%% Change of reference from Lab to Sample
@@ -53,7 +58,7 @@ end
 %   Therefore needs calculation of the inverse, and not only transpose.
 if (~isfield(samgeo,'lab2sam') || (isfield(samgeo,'lab2sam') && isempty(samgeo.lab2sam)))
     
-    if (rot_flag)
+    if (conf.only_rot)
         % Keep norm of input vectors -> keep unit lengths in the rows of 
         % the transformation matrix 
         sam2lab = [samgeo.dirx; samgeo.diry; samgeo.dirz];
@@ -66,8 +71,8 @@ if (~isfield(samgeo,'lab2sam') || (isfield(samgeo,'lab2sam') && isempty(samgeo.l
 
     lab2sam = inv(sam2lab);
 
-elseif (rot_flag)
-    error('Transformation matrix already defined. Set ''rot_flag'' to false.')
+elseif (conf.only_rot)
+    error('Transformation matrix already defined. Set ''only_rot'' to false.')
 else
     lab2sam = samgeo.lab2sam;   
 end
@@ -76,47 +81,74 @@ end
 %%% Rotation in Lab
 %%%%%%%%%%%%%%%%%%%%%
 
-ones_in = ones(size(labXYZ, 1), 1);
+if (use_diffractometer)
+    if (use_instr_tensor_4x4)
+        rottensors = omega;
+    else
+        rottensors = gtGeoDiffractometerTensor(conf.diffractometer, 'sam2lab', freevec, ...
+            'angles_basetilt', conf.basetilt, ...
+            'angles_rotation', conf.omega, ...
+            'angles_sam_tilt', conf.sample_tilts );
+    end
+    if (conf.element_wise)
+        labXYZ = permute(labXYZ, [2 3 1]);
+        labXYZ(4, 1, :) = 1;
+    else
+        labXYZ = labXYZ';
+        labXYZ(4, :) = 1;
+    end
 
-[size_omega(1), size_omega(2), size_omega(3)] = size(omega);
-% Tensor which rotates with omega from Lab back to Sample:
-if (all(size_omega(1:2) == [3 3]))
-    % This means that the rotation tensors were passed already!
-    if (size_omega(3) == 1)
-        rottensors = omega(:, :, ones_in);
+    labXYZ_rot = gtMathsMatrixProduct(rottensors, labXYZ);
+
+    if (conf.element_wise)
+        labXYZ_rot = labXYZ_rot(1:3, 1, :);
+        labXYZ_rot = permute(labXYZ_rot, [3 1 2]);
     else
-        rottensors = omega;
+        labXYZ_rot = permute(labXYZ_rot, [2 3 1]);
+        labXYZ_rot = labXYZ_rot(:, 1:3);
     end
 else
-    rotcomp    = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
-    rottensors = gtMathsRotationTensor(omega, rotcomp);
-end
+    ones_in = ones(size(labXYZ, 1), 1);
+
+    % Tensor which rotates with omega from Lab back to Sample:
+    if (use_instr_tensor_3x3)
+        % This means that the rotation tensors were passed already!
+        if (size_omega(3) == 1)
+            rottensors = omega(:, :, ones_in);
+        else
+            rottensors = omega;
+        end
+    else
+        rotcomp    = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
+        rottensors = gtMathsRotationTensor(omega, rotcomp);
+    end
 
-if (~freevec)
-    % Vectors relative to rotation axis position
-    labXYZ = labXYZ - labgeo.rotpos(ones_in, :);
-end
+    if (~freevec)
+        % Vectors relative to rotation axis position
+        labXYZ = labXYZ - labgeo.rotpos(ones_in, :);
+    end
 
-rottensors_t = reshape(rottensors, 3, []);
+    rottensors_t = reshape(rottensors, 3, []);
 
-if (element_wise)
-    % Multiply element-wise to avoid loop (applying rottensors to labXYZ)
-    labXYZ_t = [labXYZ, labXYZ, labXYZ]';
-    labXYZ_t = reshape(labXYZ_t, 3, []);
-    pro = labXYZ_t .* rottensors_t;
+    if (conf.element_wise)
+        % Multiply element-wise to avoid loop (applying rottensors to labXYZ)
+        labXYZ_t = [labXYZ, labXYZ, labXYZ]';
+        labXYZ_t = reshape(labXYZ_t, 3, []);
+        pro = labXYZ_t .* rottensors_t;
 
-    % Sum, reshape, transpose to get rotated vector
-    labXYZ_rot = reshape(sum(pro, 1), 3, [])';
-else
-    % output wll be groupped by (all rotations per center) x number centers
-    labXYZ_rot = labXYZ * rottensors_t;
-    labXYZ_rot = reshape(labXYZ_rot', 3, [])';
-end
+        % Sum, reshape, transpose to get rotated vector
+        labXYZ_rot = reshape(sum(pro, 1), 3, [])';
+    else
+        % output wll be groupped by (all rotations per center) x number centers
+        labXYZ_rot = labXYZ * rottensors_t;
+        labXYZ_rot = reshape(labXYZ_rot', 3, [])';
+    end
 
-if (~freevec)
-    ones_out = ones(size(labXYZ_rot, 1), 1);
-    % Set offset back to get rotated vector
-    labXYZ_rot = labXYZ_rot + labgeo.rotpos(ones_out, :) - samgeo.orig(ones_out, :);
+    if (~freevec)
+        ones_out = ones(size(labXYZ_rot, 1), 1);
+        % Set offset back to get rotated vector
+        labXYZ_rot = labXYZ_rot + labgeo.rotpos(ones_out, :) - samgeo.orig(ones_out, :);
+    end
 end
 
 samXYZ = labXYZ_rot * lab2sam;
diff --git a/zUtil_Geo/gtGeoProjForReconstruction.m b/zUtil_Geo/gtGeoProjForReconstruction.m
index 7584c2140e53bb3b0f05efcbf6363d3519625a42..a5d3853ef2445e97b47b049d1fe223328a6887e2 100644
--- a/zUtil_Geo/gtGeoProjForReconstruction.m
+++ b/zUtil_Geo/gtGeoProjForReconstruction.m
@@ -79,18 +79,18 @@ function proj_geom = gtGeoProjForReconstruction(projvec_sam, omega, vol_center,.
 
     % Bounding box center X,Y,Z coordinates in LAB and REC reference
     bbcent_lab  = gtGeoDet2Lab(bbcent_det, detgeo, 0);
-    bbcent_rec  = gtGeoLab2Sam(bbcent_lab, rot_omega, labgeo, recgeo, 0, 0);
+    bbcent_rec  = gtGeoLab2Sam(bbcent_lab, rot_omega, labgeo, recgeo, false);
 
     % Detector orientation vectors in REC reference
     detdiru = detgeo.detdiru ./ sqrt(sum(detgeo.detdiru .^ 2)) * detgeo.pixelsizeu;
     detdirv = detgeo.detdirv ./ sqrt(sum(detgeo.detdirv .^ 2)) * detgeo.pixelsizev;
 
-    detdiru_rec = gtGeoLab2Sam(detdiru(ones_omega, :), rot_omega, labgeo, recgeo, true, false);
-    detdirv_rec = gtGeoLab2Sam(detdirv(ones_omega, :), rot_omega, labgeo, recgeo, true, false);
+    detdiru_rec = gtGeoLab2Sam(detdiru(ones_omega, :), rot_omega, labgeo, recgeo, true);
+    detdirv_rec = gtGeoLab2Sam(detdirv(ones_omega, :), rot_omega, labgeo, recgeo, true);
 
     if (strcmpi(rectype, 'ASTRA_absorption'))
         projvec_lab = labgeo.beamdir(ones_omega, :);
-        projvec_rec = gtGeoLab2Sam(projvec_lab, rot_omega, labgeo, recgeo, true, false);
+        projvec_rec = gtGeoLab2Sam(projvec_lab, rot_omega, labgeo, recgeo, true);
     else
         % Normalised projection vector in REC reference
         projvec_rec = gtGeoSam2Sam(projvec_sam, samgeo, recgeo, true, false);
diff --git a/zUtil_Geo/gtGeoSam2Lab.m b/zUtil_Geo/gtGeoSam2Lab.m
index 6605fea299c5f1dd2c800ca2e87c12cabc675854..b6da8499aede0191c128755a89d287d5e7aaf5b3 100644
--- a/zUtil_Geo/gtGeoSam2Lab.m
+++ b/zUtil_Geo/gtGeoSam2Lab.m
@@ -1,9 +1,10 @@
-function [labXYZ, sam2lab] = gtGeoSam2Lab(samXYZ, omega, labgeo, samgeo, ...
-                             freevec, rot_flag, element_wise)
+function [labXYZ, sam2lab, rottensors] = gtGeoSam2Lab(samXYZ, omega, labgeo, samgeo, freevec, varargin)
 % GTGEOSAM2LAB  SAMPLE -> LAB reference coordinate transformation.
 %
-% [labXYZ, sam2lab] = gtGeoSam2Lab(samXYZ, om, labgeo, samgeo, ...
-%                              freevec, rot_flag)
+% [labXYZ, sam2lab] = gtGeoSam2Lab(samXYZ, om, labgeo, samgeo, freevec, varargin)
+%
+% conf = struct('only_rot', false, 'element_wise', true, ...
+%     'diffractometer', [], 'phi', 0, 'omega', 0, 'alpha', 0, 'beta', 0);
 %
 % ----------------------------------------------------------------------- 
 %
@@ -21,7 +22,8 @@ function [labXYZ, sam2lab] = gtGeoSam2Lab(samXYZ, omega, labgeo, samgeo, ...
 %   samgeo.voxsize - [x,y,z] voxel sizes in lab units
 %   freevec        - 0 for fixed vectors (e.g. location);
 %                    1 for free vectors (e.g. orientation)
-%   rot_flag       - if true, only applies rotation and the norm of the 
+% (optional)
+%   only_rot       - if true, only applies rotation and the norm of the 
 %                    input vectors will be unchanged; i.e. no scaling 
 %                    applied according to the voxel sizes (only works if 
 %                    'samgeo.sam2lab' is undefined)
@@ -36,13 +38,16 @@ function [labXYZ, sam2lab] = gtGeoSam2Lab(samXYZ, omega, labgeo, samgeo, ...
 %                    lab vectors.
 %
 
-if (~exist('rot_flag', 'var') || isempty(rot_flag))
-    rot_flag = false;
-end
+conf = struct('only_rot', false, 'element_wise', true, ...
+    'diffractometer', [], 'basetilt', [], 'omega', [], 'sample_tilts', []);
+conf = parse_pv_pairs(conf, varargin);
 
-if (~exist('element_wise', 'var') || isempty(element_wise))
-    element_wise = true;
-end
+[size_omega(1), size_omega(2), size_omega(3)] = size(omega);
+
+use_instr_tensor_4x4 = all(size_omega(1:2) == 4);
+use_instr_tensor_3x3 = all(size_omega(1:2) == 3);
+
+use_diffractometer = ~isempty(conf.diffractometer) || use_instr_tensor_4x4;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%% Change of reference from Sample to Lab
@@ -52,7 +57,7 @@ end
 %   (Lab is an orthonormal basis by definition.)
 if (~isfield(samgeo,'sam2lab') || (isfield(samgeo,'sam2lab') && isempty(samgeo.sam2lab)))
     
-    if (rot_flag)
+    if (conf.only_rot)
         % Keep norm of input vectors -> keep unit lengths in the rows of
         % the transformation matrix
         sam2lab = [samgeo.dirx; samgeo.diry; samgeo.dirz];
@@ -63,8 +68,8 @@ if (~isfield(samgeo,'sam2lab') || (isfield(samgeo,'sam2lab') && isempty(samgeo.s
             samgeo.dirz * samgeo.voxsize(3) ];
     end
     
-elseif (rot_flag)
-    error('Transformation matrix already defined. Set ''rot_flag'' to false.')
+elseif (conf.only_rot)
+    error('Transformation matrix already defined. Set ''only_rot'' to false.')
 else
     sam2lab = samgeo.sam2lab;
 end
@@ -73,51 +78,79 @@ end
 %%% Rotation in Lab
 %%%%%%%%%%%%%%%%%%%%%
 
-ones_in = ones(size(samXYZ, 1), 1);
+% Vectors in Lab reference at the origin before rotation
+labXYZprerot = samXYZ * sam2lab;
 
-[size_omega(1), size_omega(2), size_omega(3)] = size(omega);
-% Tensor which rotates with omega from Sample to Lab
-if (all(size_omega(1:2) == [3 3]))
-    % This means that the rotation tensors were passed already!
-    if (size_omega(3) == 1)
-        rottensors = omega(:, :, ones_in);
-    else
+if (use_diffractometer)
+    if (use_instr_tensor_4x4)
         rottensors = omega;
+    else
+        rottensors = gtGeoDiffractometerTensor(conf.diffractometer, 'sam2lab', freevec, ...
+            'angles_basetilt', conf.basetilt, ...
+            'angles_rotation', conf.omega, ...
+            'angles_sam_tilt', conf.sample_tilts );
+    end
+    if (conf.element_wise)
+        labXYZprerot = permute(labXYZprerot, [2 3 1]);
+        labXYZprerot(4, 1, :) = 1;
+    else
+        labXYZprerot = labXYZprerot';
+        labXYZprerot(4, :) = 1;
+    end
+
+    labXYZ = gtMathsMatrixProduct(rottensors, labXYZprerot);
+
+    if (conf.element_wise)
+        labXYZ = labXYZ(1:3, 1, :);
+        labXYZ = permute(labXYZ, [3 1 2]);
+    else
+        labXYZ = labXYZ(1:3, :, :);
+        labXYZ = reshape(labXYZ, 3, []);
+        labXYZ = permute(labXYZ, [2 3 1]);
     end
 else
-    rotcomp    = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
-    rottensors = gtMathsRotationTensor(-omega, rotcomp);
-end
+    ones_in = ones(size(samXYZ, 1), 1);
 
-% Vectors in Lab reference at the origin before rotation
-labXYZprerot = samXYZ * sam2lab;
+    % Tensor which rotates with omega from Sample to Lab
+    if (use_instr_tensor_3x3)
+        % This means that the rotation tensors were passed already!
+        if (size_omega(3) == 1)
+            rottensors = omega(:, :, ones_in);
+        else
+            rottensors = omega;
+        end
+    else
+        rotcomp    = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
+        rottensors = gtMathsRotationTensor(-omega, rotcomp);
+    end
 
-if (~freevec)
-    % Vectors relative to rot. axis position in Lab reference
-    vv           = samgeo.orig - labgeo.rotpos;
-    labXYZprerot = labXYZprerot + vv(ones_in,:);
-end
+    if (~freevec)
+        % Vectors relative to rot. axis position in Lab reference
+        vv           = samgeo.orig - labgeo.rotpos;
+        labXYZprerot = labXYZprerot + vv(ones_in,:);
+    end
 
-rottensors_t = reshape(rottensors, 3, []);
+    rottensors_t = reshape(rottensors, 3, []);
 
-if (element_wise)
-    % Multiply element-wise to avoid loop
-    labXYZprerot_t = [labXYZprerot, labXYZprerot, labXYZprerot]';
-    labXYZprerot_t = reshape(labXYZprerot_t, 3, []);
-    pro = labXYZprerot_t .* rottensors_t;
+    if (conf.element_wise)
+        % Multiply element-wise to avoid loop
+        labXYZprerot_t = [labXYZprerot, labXYZprerot, labXYZprerot]';
+        labXYZprerot_t = reshape(labXYZprerot_t, 3, []);
+        pro = labXYZprerot_t .* rottensors_t;
 
-    % Sum, reshape, transpose
-    labXYZ = reshape(sum(pro, 1), 3, [])';
-else
-    % output wll be groupped by (all rotations per center) x number centers
-    labXYZ = labXYZprerot * rottensors_t;
-    labXYZ = reshape(labXYZ', 3, [])';
-end
+        % Sum, reshape, transpose
+        labXYZ = reshape(sum(pro, 1), 3, [])';
+    else
+        % output wll be groupped by (all rotations per center) x number centers
+        labXYZ = labXYZprerot * rottensors_t;
+        labXYZ = reshape(labXYZ', 3, [])';
+    end
 
-if (~freevec)
-    ones_out = ones(size(labXYZ, 1), 1);
-    % Set offset back
-    labXYZ = labXYZ + labgeo.rotpos(ones_out, :);
+    if (~freevec)
+        ones_out = ones(size(labXYZ, 1), 1);
+        % Set offset back
+        labXYZ = labXYZ + labgeo.rotpos(ones_out, :);
+    end
 end
 
 end % of function
diff --git a/zUtil_Geo/gtGeoSamEnvInSampleRef.m b/zUtil_Geo/gtGeoSamEnvInSampleRef.m
index afa6ca8e65e5d6bb252a15a13636b9b9c5e77d92..0031c8e8c96f380d6f30edfce1ce7da0390a88f1 100644
--- a/zUtil_Geo/gtGeoSamEnvInSampleRef.m
+++ b/zUtil_Geo/gtGeoSamEnvInSampleRef.m
@@ -1,5 +1,5 @@
-function samenv = gtGeoSamEnvInSampleRef(labgeo,samgeo)
-%     samenv = gtSampleVolInSampleRef(labgeo,samgeo)
+function samenv = gtGeoSamEnvInSampleRef(labgeo, samgeo)
+%     samenv = gtSampleVolInSampleRef(labgeo, samgeo)
 %     ----------------------------------------------
 %     Returns the coordinates of the irradiated sample volume in the Sample 
 %     reference.
@@ -17,19 +17,18 @@ function samenv = gtGeoSamEnvInSampleRef(labgeo,samgeo)
 %                                  size of the three)
 %       samenv.axis = <double>     rotation axis direction
 
-% Envelope top and bottom point in Lab
-toppoint = labgeo.rotpos + labgeo.rotdir*labgeo.samenvtop;
-botpoint = labgeo.rotpos + labgeo.rotdir*labgeo.samenvbot;
+    % Envelope top and bottom point in Lab
+    toppoint = labgeo.rotpos + labgeo.rotdir * labgeo.samenvtop;
+    botpoint = labgeo.rotpos + labgeo.rotdir * labgeo.samenvbot;
 
-% Envelope top and bottom point in Sample
-samenv.top  = gtGeoLab2Sam(toppoint,0,labgeo,samgeo,0,0);
-samenv.bot  = gtGeoLab2Sam(botpoint,0,labgeo,samgeo,0,0);
+    % Envelope top and bottom point in Sample
+    samenv.top  = gtGeoLab2Sam(toppoint, 0, labgeo, samgeo, false);
+    samenv.bot  = gtGeoLab2Sam(botpoint, 0, labgeo, samgeo, false);
 
-% Radius
-samenv.rad  = labgeo.samenvrad/max(samgeo.voxsize);
+    % Radius
+    samenv.rad  = labgeo.samenvrad / max(samgeo.voxsize);
 
-% Normalised rotation axis direction
-samenv.axis = gtGeoLab2Sam(labgeo.rotdir,0,labgeo,samgeo,1,1);
-samenv.axis = samenv.axis/sqrt(sum(samenv.axis.*samenv.axis,2));
-
-end % end of function
+    % Normalised rotation axis direction
+    samenv.axis = gtGeoLab2Sam(labgeo.rotdir, 0, labgeo, samgeo, true, 'only_rot', true);
+    samenv.axis = samenv.axis / sqrt(sum(samenv.axis .* samenv.axis, 2));
+end
diff --git a/zUtil_Indexter/gtIndexDifSpots_2.m b/zUtil_Indexter/gtIndexDifSpots_2.m
index 05c942214b8af85560b16ebf484b59158052a06c..9b4349fbd57afc3ba9a4983b56b848d1c8cfe525 100644
--- a/zUtil_Indexter/gtIndexDifSpots_2.m
+++ b/zUtil_Indexter/gtIndexDifSpots_2.m
@@ -86,7 +86,7 @@ function [gid, conf, sp] = gtIndexDifSpots_2(grain, sp, parameters, fsimID, lim)
         xyzB  = [sp.samcentXB sp.samcentYB sp.samcentZB];
         xyz   = [xyzA; xyzB];
 
-        uvwLab = gtGeoSam2Lab(xyz, sp.om, labgeo, parameters.samgeo, 0, false);
+        uvwLab = gtGeoSam2Lab(xyz, sp.om, labgeo, parameters.samgeo, false);
         sp.uvw(:, 1:2) = gtGeoLab2Det(uvwLab, detgeo, 0);
 
         pairsinfo = true;
diff --git a/zUtil_Indexter/gtIndexFwdSim.m b/zUtil_Indexter/gtIndexFwdSim.m
index 4a35a068b586706cdf7f633a22791dff7dbb9375..cf807a93e3191fde7c70ba8bfa735fcb444a6e51 100644
--- a/zUtil_Indexter/gtIndexFwdSim.m
+++ b/zUtil_Indexter/gtIndexFwdSim.m
@@ -57,7 +57,7 @@ fwdSim.rot(:, :, ind4) = rot4(:, :, ind4, 4);
 
 % Diffraction vector
 fwdSim.dveclab = gtFedPredictDiffVecMultiple(fwdSim.pllab, labgeo.beamdir');
-fwdSim.dvecsam = gtGeoLab2Sam(fwdSim.dveclab', fwdSim.om, labgeo, samgeo, 1, 1)';
+fwdSim.dvecsam = gtGeoLab2Sam(fwdSim.dveclab', fwdSim.om, labgeo, samgeo, true, 'only_rot', true)';
 
 % Make sure grain center is extended
 gcsam = gcsam(:, ones(size(fwdSim.dveclab, 2), 1));
diff --git a/zUtil_Indexter/gtIndexFwdSimGrains.m b/zUtil_Indexter/gtIndexFwdSimGrains.m
index 6c58738539bf96bf3cd7a617354989e3638e4c4e..88a2d48fe847bdc0710aaf337b3a3b9197cc95e7 100644
--- a/zUtil_Indexter/gtIndexFwdSimGrains.m
+++ b/zUtil_Indexter/gtIndexFwdSimGrains.m
@@ -215,7 +215,7 @@ for ii = 1:length(grain)
 
     % Diffraction vectors in Lab and Sample
     dveclab = gtFedPredictDiffVecMultiple(pllab4, labgeo.beamdir');
-    dvecsam = gtGeoLab2Sam(dveclab', om, labgeo, samgeo, 1, 1)';
+    dvecsam = gtGeoLab2Sam(dveclab', om, labgeo, samgeo, true, 'only_rot', true)';
 
     % Extend grain center
     gc_sam = grain{ii}.center';
diff --git a/zUtil_Maths/gtMathsRotationTranslationTensor.m b/zUtil_Maths/gtMathsRotationTranslationTensor.m
new file mode 100644
index 0000000000000000000000000000000000000000..fd424e937af551a350fee18e528df64f28677f9f
--- /dev/null
+++ b/zUtil_Maths/gtMathsRotationTranslationTensor.m
@@ -0,0 +1,22 @@
+function t = gtMathsRotationTranslationTensor(t_axis, t_angles, t_orig, freevec)
+% function t = gtMathsRotationTranslationTensor(t_axis, t_angles, t_orig, freevec)
+%
+% It returns a <4 x 4 x n> stack of rotation + translation tensors
+
+    rot_comp = gtMathsRotationMatrixComp(t_axis', 'col');
+    t = eye(4);
+    t = repmat(t, 1, 1, numel(t_angles));
+    t(1:3, 1:3, :) = gtMathsRotationTensor(t_angles, rot_comp);
+
+    t_tr_pre = eye(4);
+    if (~freevec)
+        t_tr_pre(:, 4) = [-t_orig'; 1];
+    end
+
+    t_tr_post = eye(4);
+    if (~freevec)
+        t_tr_post(:, 4) = [t_orig'; 1];
+    end
+
+    t = gtMathsMatrixProduct(t_tr_post, gtMathsMatrixProduct(t, t_tr_pre));
+end