diff --git a/6_rendering/gtAaxisCmap.m b/6_rendering/gtAaxisCmap.m
index 4621c8e8c16689872a8efe7db792d3cd84e32f3e..0a820f2d88480cee2065a52f40ddf856ab0ed192 100644
--- a/6_rendering/gtAaxisCmap.m
+++ b/6_rendering/gtAaxisCmap.m
@@ -16,6 +16,8 @@ function [a_axis, a_axismap] = gtAaxisCmap(phaseid, reflection)
 %     r_vectors -> a-axis -> phi/psi -> hsl -> rgb colormap
 %     use psi to give brightness; phi to give colour
 %
+%     Version 003 19-10-2012 by YGuilhem
+%       Factorized using gtVectorOrientationColor
 %
 %     Version 002 28-06-2012 by LNervo
 %       Update to version 2 of parameters; cleaning formatting
@@ -55,10 +57,9 @@ normalised_hkls = all_hkls./(repmat(tmp,1,3));
 all_hkls = [];
 all_hkls = all_hkils;
 
-
 for i=1:size(r_vectors,1)
-    g = Rod2g(r_vectors(i,2:4));
-    all_normals = (g * normalised_hkls')';
+    crys2sam = Rod2g(r_vectors(i,2:4));
+    all_normals = (crys2sam * normalised_hkls')';
     a_axis(i,1)=r_vectors(i,1);
     if all_normals(6,1)>=0
         a_axis(i,2)=all_normals(6,1);
@@ -71,46 +72,18 @@ for i=1:size(r_vectors,1)
     end
 end
 
+save('a_axis.mat', 'a_axis'); % list of the a axis
 
-save a_axis a_axis % list of the a axis
-
-
-psi=acosd(a_axis(:,4));
-phi=atan2(a_axis(:,3), a_axis(:,2));
-
-
-
-%get everything in degrees, and on the right interval
-phi=phi*180/pi;
-% %want all psi between 0 and 90
-% psi(find(psi>90))=180-psi(find(psi>90));
-% %if we change psi, should also change phi
-% phi(find(psi>90))=phi(find(psi>90))+180;
+% Get the orientation color defined by HSL color code, translated to RGB
+rgb = gtVectorOrientationColor(a_axis(:, 2:4));
 
-dummy=find(psi>90);
-psi(dummy)=180-psi(dummy);
-%if we change psi, should also change phi
-phi(dummy)=phi(dummy)+180;
-
-phi(find(phi<0))=phi(find(phi<0))+360;
-phi(find(phi>360))=phi(find(phi>360))-360;
-
-%use hsl2rgb.m
-%hue
-h=phi/360;
-%saturation
-s=ones(size(h));
-%lightness (white //z, black in x-y plane)
-l=0.9-(0.8*(psi/90));
-
-rgb=hsl2rgb([h s l]);
-
-%add black background
-a_axismap=zeros(size(rgb, 1)+1, 3);
-a_axismap(2:end, :)=rgb;
+% Add black background and expand list length of max(grainid)+1
+a_axismap = zeros(max(a_axis(:, 1))+1, 3);
+for i=1:length(a_axis)
+    a_axismap(a_axis(i,1)+1, :) = rgb(i, :);
+end
 
-%save the a-axis colourmap
-save a_axismap a_axismap
+% Save the a-axis colourmap
+save('a_axismap.mat', 'a_axismap');
 
 end % end of function
-
diff --git a/6_rendering/gtCaxisCmap.m b/6_rendering/gtCaxisCmap.m
index 69d2b08d67bf652c1faa6b6d8850111f4bc135c1..c5a46977fe103068bffaf99e5d2f53ffe5d26d05 100644
--- a/6_rendering/gtCaxisCmap.m
+++ b/6_rendering/gtCaxisCmap.m
@@ -16,6 +16,8 @@ function [c_axis, c_axismap] = gtCaxisCmap(phaseid, reflection)
 %     r_vectors -> c-axis -> phi/psi -> hsl -> rgb colormap
 %     use psi to give brightness; phi to give colour
 %
+%     Version 003 19-10-2012 by YGuilhem
+%       Factorized using gtVectorOrientationColor
 %
 %     Version 002 28-06-2012 by LNervo
 %       Update to version 2 of parameters; cleaning formatting
@@ -55,11 +57,9 @@ normalised_hkls = all_hkls./(repmat(tmp,1,3));
 all_hkls = [];
 all_hkls = all_hkils;
 
-
 for i=1:size(r_vectors,1)
-    g = Rod2g(r_vectors(i,2:4));
-   
-    all_normals = (g * normalised_hkls')';
+    crys2sam = Rod2g(r_vectors(i,2:4));
+    all_normals = (crys2sam * normalised_hkls')';
     c_axis(i,1)=r_vectors(i,1);
     if all_normals(2,3)>=0
         c_axis(i,2)=all_normals(2,1);
@@ -72,45 +72,18 @@ for i=1:size(r_vectors,1)
     end
 end
 
-
 save('c_axis.mat','c_axis'); % list of the c axis
 
+% Get the orientation color defined by HSL color code, translated to RGB
+rgb = gtVectorOrientationColor(c_axis(:, 2:4));
 
-psi=acosd(c_axis(:,4));
-phi=atan2(c_axis(:,3), c_axis(:,2));
-% phi=atan2(c_axis(:,2), c_axis(:,3));
-%get everything in degrees, and on the right interval
-phi=phi*180/pi;
-% %want all psi between 0 and 90
-% psi(find(psi>90))=180-psi(find(psi>90));
-% %if we change psi, should also change phi
-% phi(find(psi>90))=phi(find(psi>90))+180;
-
-dummy=find(psi>90);
-psi(dummy)=180-psi(dummy);
-%if we change psi, should also change phi
-phi(dummy)=phi(dummy)+180;
-
-phi(find(phi<0))=phi(find(phi<0))+360;
-phi(find(phi>360))=phi(find(phi>360))-360;
-
-%use hsl2rgb.m
-%hue
-h=phi/360;
-%saturation
-s=ones(size(h));
-%lightness (white //z, black in x-y plane)
-l=0.9-(0.8*(psi/90));
-
-rgb=hsl2rgb([h s l]);
-
-%add black background and expand list length of max(grainid)+1
-c_axismap=zeros(max(c_axis(:,1))+1, 3);
+% Add black background and expand list length of max(grainid)+1
+c_axismap = zeros(max(c_axis(:, 1))+1, 3);
 for i=1:length(c_axis)
-  c_axismap(c_axis(i,1)+1, :)=rgb(i,:);
+    c_axismap(c_axis(i,1)+1, :) = rgb(i, :);
 end
 
-%save the c-axis colourmap
+% Save the c-axis colourmap
 save('c_axismap.mat','c_axismap');
 
 end % end of function
diff --git a/zUtil_Analysis/gtAnalyseCrackPlane.m b/zUtil_Analysis/gtAnalyseCrackPlane.m
index cfb91ca6af8499a4ce518e4d66ea9ffbc9621bbf..5866c39b0e623755427a0ffaf61175dce233d1d5 100644
--- a/zUtil_Analysis/gtAnalyseCrackPlane.m
+++ b/zUtil_Analysis/gtAnalyseCrackPlane.m
@@ -1,11 +1,8 @@
-function [data, vol, map, surface_view]=gtAnalyseCrackPlane(crack);
-
+function [data, vol, map, surface_view] = gtAnalyseCrackPlane(crack)
 %returns the crack plane data as a list [x y z nx ny nz]
 %also returns a labelled volume and tle corresponding colourmap
 %surface_view is a top down image
 
-
-
 %read the x,y,z points of the crack
 [y,x,z]=ind2sub(size(crack),find(crack)); %voxels of crack
 
@@ -13,80 +10,54 @@ boxsize=5;
 %data out =[x y z nx ny nz]
 data=zeros(length(find(crack)), 6);
 
-
 for i=1:length(x)
-  
-  xx=x(i); yy=y(i); zz=z(i);
-  
-
-dummy=find(x>xx-boxsize & x<xx+boxsize & y>yy-boxsize & y<yy+boxsize & z>zz-boxsize & z<zz+boxsize);
-
-if length(dummy)>3
-  
-%this is right in tomography coordinates (ie as in volview, x l->r, y
-%top->bot, z up->down)
-
-  [origin, a]=lsplane([x(dummy) y(dummy) z(dummy)]);
-  
-%convert to instrument coordinates to agree with gtCaxisCmap conventions
-a=a([2 1 3]);
-a(3)=-a(3);
+    xx=x(i); yy=y(i); zz=z(i);
+    dummy=find(x>xx-boxsize & x<xx+boxsize & y>yy-boxsize & y<yy+boxsize & z>zz-boxsize & z<zz+boxsize);
 
-data(i, :)=[xx yy zz a(1) a(2) a(3)];
+    if length(dummy)>3
+        %this is right in tomography coordinates (ie as in volview, x l->r, y
+        %top->bot, z up->down)
+        [origin, a]=lsplane([x(dummy) y(dummy) z(dummy)]);
 
-if mod(i, 1000)==0
-  100*i/length(x)
-end
+        %convert to instrument coordinates to agree with gtCaxisCmap conventions
+        a=a([2 1 3]);
+        a(3)=-a(3);
 
-end
+        data(i, :)=[xx yy zz a(1) a(2) a(3)];
 
+        if mod(i, 1000)==0
+            100*i/length(x)
+        end
+    end
 end
 
 %cut the empty bits of data
-data(find(data(:,1)==0), :)=[]; 
+data(find(data(:,1)==0), :)=[];
 
 %now produce a labelled volume, and the corresponding colormap
 %volume
 vol=zeros(size(crack));
 for i=1:length(data)
-  vol(data(i,2), data(i,1), data(i,3))=i;
+    vol(data(i,2), data(i,1), data(i,3))=i;
 end
-%colormap
-%as in gtCaxisCmap
-psi=acosd(data(:,6));
-phi=atan2(data(:,5), data(:,4));
 
-%get everything in degrees, and on the right interval
-phi=phi*180/pi;
-dummy=find(psi>90);
-psi(dummy)=180-psi(dummy);
-phi(dummy)=phi(dummy)+180;
-phi(find(phi<0))=phi(find(phi<0))+360;
-phi(find(phi>360))=phi(find(phi>360))-360;
+% Colormap
+rgb = gtVectorOrientationColor(data(:, 4:6));
 
-%hue
-h=phi/360;
-%saturation
-s=ones(size(h));
-%lightness (white //z, black in x-y plane)
-l=0.9-(0.8*(psi/90));
-rgb=hsl2rgb([h s l]);
 %add black background
 map=zeros(size(data,1)+1, 3);
 map(2:end, :)=rgb; %ie rgb plus zero background
 
-
 %read out "bird's eye" view of coloured crack;
 surface_view=zeros(size(crack, 1), size(crack, 2));
 for y=1:size(crack, 1)
     for x=1:size(crack, 2)
-a=squeeze(vol(y,x,:));
-a=a(find(a, 1, 'first'));
-if ~isempty(a)
-    surface_view(y,x)=a;
-end
+        a=squeeze(vol(y,x,:));
+        a=a(find(a, 1, 'first'));
+        if ~isempty(a)
+            surface_view(y,x)=a;
+        end
     end
 end
 
-
-
+end % end of function
diff --git a/zUtil_Imaging/gtVectorOrientationColor.m b/zUtil_Imaging/gtVectorOrientationColor.m
new file mode 100644
index 0000000000000000000000000000000000000000..cf6909401683ee80fcbd3b325e3d747e759bab08
--- /dev/null
+++ b/zUtil_Imaging/gtVectorOrientationColor.m
@@ -0,0 +1,35 @@
+function rgb = gtVectorOrientationColor(vectors)
+% GTVECTORORIENTATIONCOLOR Convert a list of unit vectors into RGB color values
+%
+%     rgb = gtVectorOrientationColor(vectors)
+%     -------------------------------------------------------------------------
+%     The output RGB color values describe roientation with HSL code which use:
+%     - psi to give brightness
+%     - phi to give color
+%
+%     INPUT:
+%       vectors = <double>  M row vectors (Mx3 matrix)
+%
+%     OUTPUT:
+%       rgb     = <double>  Red-Green-Blue color values 
+%
+%     Version 001 16-10-2012 by YGuilhem
+
+    psi = acos(vectors(:, 3));
+    phi = atan2(vectors(:, 2), vectors(:, 1));
+
+    % Get everything in the right interval
+    dummy = find(psi > pi/2);
+    psi(dummy) = pi - psi(dummy);
+
+    % If we change psi, should also change phi
+    phi(dummy) = phi(dummy) + pi;
+    phi = mod(phi, pi);
+
+    h = phi/pi;                   % Hue
+    s = ones(size(h));            % Saturation
+    l = 0.9 - (0.8*(psi/(pi/2))); % Lightness (white //z, black in x-y plane)
+
+    rgb = hsl2rgb([h s l]);   % Turn HSL into RGB
+
+end % end of function