diff --git a/5_reconstruction/gtFindSimilarGrains.m b/5_reconstruction/gtFindSimilarGrains.m
index 003bd1d6e43bb42c872ed53f7c55ea85edd87b58..8693ab86a74abb7fedf41a7d1d5d52368dce6356 100644
--- a/5_reconstruction/gtFindSimilarGrains.m
+++ b/5_reconstruction/gtFindSimilarGrains.m
@@ -24,8 +24,8 @@ map2.pos=map2.pos-repmat(mean(map2.pos,1), length(map2.pos), 1);
 map2.pos=map2.pos*parameters2.acq.pixelsize/parameters1.acq.pixelsize;
 
 % Compute all orientation matrices g on each map
-all_g1 = gtMathsRod2OriMat(map1.R_vec);
-all_g2 = gtMathsRod2OriMat(map2.R_vec);
+all_g1 = gtMathsRod2OriMat(map1.R_vec.');
+all_g2 = gtMathsRod2OriMat(map2.R_vec.');
 
 % Go through grain map 1, looking for matches in grain map 2
 output = [];
diff --git a/6_rendering/gtCrystAxisPoles.m b/6_rendering/gtCrystAxisPoles.m
index b5ac40b6674fa4ff3439a3979abffc402967f0df..d8d193b17697fe5ab3a20f662fd88bdf6f9d5463 100644
--- a/6_rendering/gtCrystAxisPoles.m
+++ b/6_rendering/gtCrystAxisPoles.m
@@ -68,7 +68,7 @@ poles = zeros(0,4);
 valid_poles = zeros(0,4);
 
 % Compute all orientation matrices g
-all_g = gtMathsRod2OriMat(r_vectors);
+all_g = gtMathsRod2OriMat(r_vectors.');
 
 % Express normalised hkls in cartesian SAMPLE CS
 all_normals = gtVectorCryst2Lab(normalised_hkls, all_g);
diff --git a/6_rendering/gtIPFCmap.m b/6_rendering/gtIPFCmap.m
index 1a4f2e4c59f9f385fc5675c123b0292a4191c66b..89bc23ca878834a830c105d3b5aee40dbc7be4fc 100644
--- a/6_rendering/gtIPFCmap.m
+++ b/6_rendering/gtIPFCmap.m
@@ -65,16 +65,15 @@ elseif exist(fullfile(phaseDir, 'r_vectors.mat'), 'file')
 elseif exist(fullfile(phaseDir, 'index.mat'), 'file')
     grain = [];
     load(fullfile(phaseDir, 'index.mat'), 'grain');
-    r_vectors = gtIndexAllGrainValues(grain, 'R_vector', [], 1:3, []);
+    r_vectors = gtIndexAllGrainValues(grain, 'R_vector', [], 1, 1:3);
 else
     gtError('gtIPFCmap:missing_file', 'Could not find r_vectors.mat or index.mat file!');
 end
 
-% Resize r_vectors if needed and transpose to column vectors
-if size(r_vectors,2) == 4
+% Resize r_vectors if needed
+if size(r_vectors, 2) == 4
     r_vectors = r_vectors(:, 2:4);
 end
-r_vectors = r_vectors.';
 
 % Get the crystal system
 if isempty(par.crystal_system)
@@ -84,8 +83,8 @@ end
 
 % Compute all orientation matrices g
 %   Gcry = g Gsam
-%   g = U = gtMathsRod2OriMat(R_vector)
-all_g = gtMathsRod2OriMat(r_vectors);
+%   g = U = gtMathsRod2OriMat(r_vectors.')
+all_g = gtMathsRod2OriMat(r_vectors.');
 
 % Compute LDc = LD expressed in the crystal CS
 LDc = gtVectorLab2Cryst(LD, all_g);
diff --git a/6_rendering/gtPlotHexagon.m b/6_rendering/gtPlotHexagon.m
index 09b6505d7e7a09dbf53653ec4ea2843ee4242673..96638148ef2f48811a82117caabfc3954688ff54 100644
--- a/6_rendering/gtPlotHexagon.m
+++ b/6_rendering/gtPlotHexagon.m
@@ -90,7 +90,7 @@ hex_faces = data.faces_end;
 vertices = zeros(12,3,size(r_vectors,1));
 
 % Compute all orientation matrices g
-all_g = gtMathsRod2OriMat(r_vectors(:, 2:4));
+all_g = gtMathsRod2OriMat(r_vectors(:, 2:4).');
 
 % Express corners3 of each grain in cartesian SAMPLE CS
 all_corners3a = gtVectorCryst2Lab(corners3, all_g);
diff --git a/6_rendering/gtRvectorCmap.m b/6_rendering/gtRvectorCmap.m
index 722312c2abc6dd8b2582d65679c1a7706ec20e73..f19868d628ff501ef9506b566078afdcb71aa07b 100644
--- a/6_rendering/gtRvectorCmap.m
+++ b/6_rendering/gtRvectorCmap.m
@@ -15,8 +15,8 @@ else
         'index.mat');
     try
     load(indexFile);
-    r_vectors=gtIndexAllGrainValues(grain, 'R_vector', '', 1:3, '');
-    r_vectors=[zeros(size(r_vectors, 1), 1) r_vectors];
+    r_vectors = gtIndexAllGrainValues(grain, 'R_vector', [], 1, 1:3);
+    r_vectors = [zeros(size(r_vectors, 1), 1) r_vectors];
     catch mexc
         gtPrintException('neither r_vectors.mat nor index.mat was found :-(');
         rmap=[]
diff --git a/6_rendering/gtShowSampleSurface.m b/6_rendering/gtShowSampleSurface.m
index 5c8e0ddd05f1849b2b89dde82137d156f99b47ec..6e97392810221fca002bb68c26794b469ad86562 100644
--- a/6_rendering/gtShowSampleSurface.m
+++ b/6_rendering/gtShowSampleSurface.m
@@ -119,7 +119,7 @@ disp('note that the section plane that we look at is y (--> left to right) z (--
 disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~');
 
 % Compute all orientation matrices g
-all_g = gtMathsRod2OriMat(r_vectors(:, 2:4));
+all_g = gtMathsRod2OriMat(r_vectors(:, 2:4).');
 
 % Express every c axis in cartesian SAMPLE CS
 c_axes = gtVectorCryst2Lab([0 0 1], all_g);
diff --git a/zUtil_Boundaries/gtAnalyseBoundaryIndexPlane.m b/zUtil_Boundaries/gtAnalyseBoundaryIndexPlane.m
index 02f1d34de3bd54c942652bc199b19cba90c18be3..92e79983c40a3ee8f3ea281759da43350b9c8ade 100644
--- a/zUtil_Boundaries/gtAnalyseBoundaryIndexPlane.m
+++ b/zUtil_Boundaries/gtAnalyseBoundaryIndexPlane.m
@@ -46,8 +46,13 @@ end
 tmp = sqrt(sum((all_hkls.*all_hkls),2));
 normalised_hkls = all_hkls./(repmat(tmp,1,3));
 
+% Extract all the Rodrigues vectors into column vectors
+%R_vectors1 = gtIndexAllGrainValues(list, 'grain1_R_vector', [], 1, 1:3); % not used
+R_vectors2 = gtIndexAllGrainValues(list, 'grain2_R_vector', [], 1, 1:3);
+
 % Compute all orientation matrices g
-all_g = gtMathsRod2OriMat(R);
+%all_g1 = gtMathsRod2OriMat(R_vectors1.'); % not used
+all_g2 = gtMathsRod2OriMat(R_vectors2.');
 
 for ii = 1:length(list) % list(i) is the id of the boundary in the volume
   
@@ -55,7 +60,6 @@ for ii = 1:length(list) % list(i) is the id of the boundary in the volume
 
     % Get boundaries_structure data
     grainids=[boundaries_structure(list(ii)).grain1 boundaries_structure(list(ii)).grain2];
-    R=[boundaries_structure(list(ii)).grain1_R_vector; boundaries_structure(list(ii)).grain2_R_vector];
   
     % Skip external boundaries / small boundaries
     if any(grainids==0) || boundaries_structure(list(ii)).count<4
@@ -72,7 +76,7 @@ for ii = 1:length(list) % list(i) is the id of the boundary in the volume
         % find(mod(x,2)==0 & mod(y,2)==0 & mod(z,2)==0)
         % or by using a grid of the boxsize (/2?) to create subregions
 
-        g = all_g(:, :, ii);
+        g = all_g2(:, :, ii);
 
         % Loop through this volume, calculating local character
         % may want to do this in larger steps
diff --git a/zUtil_Indexter/gtINDEXMatchGrains.m b/zUtil_Indexter/gtINDEXMatchGrains.m
index 3cde44dc8d4f8e1dca0e8e44a62ce0a70a039202..f580d9c3cc311eef3dbeef8feaa77a124f12a9bf 100644
--- a/zUtil_Indexter/gtINDEXMatchGrains.m
+++ b/zUtil_Indexter/gtINDEXMatchGrains.m
@@ -304,7 +304,7 @@ gr2.Rlines = cell(nof_grains2,1);
 gr2.Rids = cell(nof_grains2,1);
 
 % All orientation matrices in dataset2
-all_g2 = gtMathsRod2OriMat(gr2.R_vector);
+all_g2 = gtMathsRod2OriMat(gr2.R_vector.');
 
 if (any(gtIndexAllGrainValues(grain1,'R_onedge',[],1,[])) || ~isempty(drot))
    
@@ -368,7 +368,7 @@ gr1.Rlines = cell(nof_grains1,1);
 gr1.Rids   = cell(nof_grains1,1);
 
 % All orientation matrices in dataset1
-all_g1 = gtMathsRod2OriMat(gr1.R_vector);
+all_g1 = gtMathsRod2OriMat(gr1.R_vector.');
 
 % Plane normals in the Sample ref. of dataset2
 all_pl_sam1 = gtVectorCryst2Lab(pl_cry.', all_g1);
diff --git a/zUtil_Indexter/gtIndexRefineGrains.m b/zUtil_Indexter/gtIndexRefineGrains.m
index f24816a291de3519babbd38912882996ec9f8793..a751c07dbc0ee5be2f22a3d3d16457a2a0979071 100644
--- a/zUtil_Indexter/gtIndexRefineGrains.m
+++ b/zUtil_Indexter/gtIndexRefineGrains.m
@@ -43,8 +43,11 @@ gtodel = [];
 indshkldirs_cry = gtCrystHKL2Cartesian(vertcat(cryst.shklfam{:})', ...
                                        cryst.Bmat)';
 
+% Extract all the Rodrigues vectors into column vectors
+R_vectors = gtIndexAllGrainValues(grain, 'R_vector', [], 1, 1:3);
+
 % Compute all orientation matrices g
-all_g = gtMathsRod2OriMat({grain.R_vector});
+all_g = gtMathsRod2OriMat(R_vectors.');
 
 % Express normalised signed hkl vectors in cartesian SAMPLE CS
 all_shkldirs = gtVectorCryst2Lab(indshkldirs_cry, all_g); 
diff --git a/zUtil_Indexter/gtIndexStatisticalFit.m b/zUtil_Indexter/gtIndexStatisticalFit.m
index 433ad3954d562eada8db6124f41fb3537b9c3132..9f4f55e51432d3eead31598be57f206f5617e9cd 100644
--- a/zUtil_Indexter/gtIndexStatisticalFit.m
+++ b/zUtil_Indexter/gtIndexStatisticalFit.m
@@ -51,11 +51,14 @@ gnesslist  = cell(nof_srefs,1);
 % All possible normalised signed hkl vectors in cartesian CRYSTAL coordinates
 shkldirs_cry = gtCrystHKL2Cartesian(vertcat(cryst.shklfam{:})', cryst.Bmat)';
 
-% Compute orientation matrix g
-all_g = gtMathsRod2OriMat({grain.R_vector});
+% Extract all the Rodrigues vectors into column vectors
+R_vectors = gtIndexAllGrainValues(grain, 'R_vector', [], 1, 1:3);
+
+% Compute all orientation matrices g
+all_g = gtMathsRod2OriMat(R_vectors.');
 
 % Express normalised signed hkl vectors in cartesian SAMPLE CS
-all_shkldirs = gtVectorCryst2Lab(indshkldirs_cry, all_g); 
+all_shkldirs = gtVectorCryst2Lab(shkldirs_cry, all_g); 
 
 for ii = 1:nof_grains
     % Store normalised signed hkl vectors in cartesian SAMPLE CS
diff --git a/zUtil_Strain/gtStrainDrawPrincipalStrainVectors.m b/zUtil_Strain/gtStrainDrawPrincipalStrainVectors.m
index 23127ac8773759b9967f59a6e861d2c3fc1c0b56..281c5590fa81779c2033e26965ada9748c047e3b 100644
--- a/zUtil_Strain/gtStrainDrawPrincipalStrainVectors.m
+++ b/zUtil_Strain/gtStrainDrawPrincipalStrainVectors.m
@@ -44,9 +44,13 @@ set(h,'FaceColor','b','EdgeColor','w')
 h = findobj(gca,'Type','line');
 set(h,'Color','b','Linewidth',2)
 
+% Extract all the Rodrigues vectors into column vectors
+R_vectors = gtIndexAllGrainValues(grains, 'R_vector', [], 1, 1:3);
 
-figure('name','Hist grain orientation')
-all_g = gtMathsRod2OriMat({grains.R_vector});
+% Compute all orientation matrices g
+all_g = gtMathsRod2OriMat(R_vectors.');
+
+figure('name','Hist grain orientation');
 for ii = 1:nof_grains
     R = all_g(:, :, ii).';
     v1 = R(1,3);
diff --git a/zUtil_Strain/gtStrainFitElasticConstants.m b/zUtil_Strain/gtStrainFitElasticConstants.m
index 6fccd921da65903aac32a29f0fb73a24ee5db17e..33b58b1edc0af70d44094325e439ae574fc94aa5 100644
--- a/zUtil_Strain/gtStrainFitElasticConstants.m
+++ b/zUtil_Strain/gtStrainFitElasticConstants.m
@@ -58,8 +58,11 @@ AI=AI*ps;
 % end
 % toc
 
+% Extract all the Rodrigues vectors into column vectors
+R_vectors = gtIndexAllGrainValues(grain, 'R_vector', [], 1, 1:3);
+
 % Compute all orientation matrices g
-all_g = gtMathsRod2OriMat({grain.R_vector}); % transform from LAB to CRYSTAL
+all_g = gtMathsRod2OriMat(R_vectors.');
 
 % Express zlab in LAB CS
 all_zcryst = gtVectorLab2Cryst(zlab, all_g);
diff --git a/zUtil_Strain/gtStrainFitElasticConstants_test.m b/zUtil_Strain/gtStrainFitElasticConstants_test.m
index cfd92386bd3295b72bb2c924f43c3547ddfaa078..1a23b47ee8b01a3e466a51e8c7f3e2ad430650ce 100644
--- a/zUtil_Strain/gtStrainFitElasticConstants_test.m
+++ b/zUtil_Strain/gtStrainFitElasticConstants_test.m
@@ -52,8 +52,11 @@ AI=AI*ps;
 % end
 % toc
 
+% Extract all the Rodrigues vectors into column vectors
+R_vectors = gtIndexAllGrainValues(grain, 'R_vector', [], 1, 1:3);
+
 % Compute all orientation matrices g
-all_g = gtMathsRod2OriMat({grain.R_vector}); % transform from LAB to CRYSTAL
+all_g = gtMathsRod2OriMat(R_vectors.');
 
 % Express zlab in LAB CS
 all_zcryst = gtVectorLab2Cryst(zlab, all_g);
diff --git a/zUtil_Taper/gtTaperConvertFableOutputToDCT.m b/zUtil_Taper/gtTaperConvertFableOutputToDCT.m
index c94e94a4f3f45296c5db1592ada241db2d5368cc..51775e15c7f2d018ff18d718d86c37863ea20120 100644
--- a/zUtil_Taper/gtTaperConvertFableOutputToDCT.m
+++ b/zUtil_Taper/gtTaperConvertFableOutputToDCT.m
@@ -126,7 +126,7 @@ for ii = 1:length(cen)
     grain{ii}.stat.bbysmean = 1;
 end
 
-rvec = gtIndexAllGrainValues(grain,'R_vector',[],1:3,[]);
+rvec = gtIndexAllGrainValues(grain,'R_vector',[],1,1:3);
 rotmat = gtMathsRod2OriMat( rvec.' );
 
 for ii = 1:length(grain)