diff --git a/6_rendering/gtHexagonalUnitCell.m b/6_rendering/gtHexagonalUnitCell.m
index 8a4862d7cbc74477737a2742eb31d3eca8d0c271..89260cc469f4f6bb4fb2a1f06f883b80c3b17544 100644
--- a/6_rendering/gtHexagonalUnitCell.m
+++ b/6_rendering/gtHexagonalUnitCell.m
@@ -1,6 +1,18 @@
-function vertices = gtHexagonalUnitCell(a,c,draw)
+function data = gtHexagonalUnitCell(varargin)
 
-% create a hexagonal grid (radius <r>)
+app.centered = true;
+app.draw     = true;
+app.ratio    = 1.5857;
+app.ind      = 0;
+app.vertices = false;
+app.numbers  = false;
+app.patch    = false;
+app.axes     = false;
+
+app = parse_pv_pairs(app, varargin);
+
+c=1;
+a=1/app.ratio;
 
 % apothem = c2 = r*sqrt(3)/2
 cx=a*sqrt(3);
@@ -10,83 +22,111 @@ cv=a*cosd(v);
 sv=a*sind(v);
 x0 = 0*cx+cv;
 y0 = 0*cx+sv;
-z0 = repmat(-c/2,1,length(v));
-z1 = repmat(+c/2,1,length(v));
-
+if app.centered
+    z0 = repmat(-c/2,1,length(v));
+    z1 = repmat(+c/2,1,length(v));
+else
+    z0 = zeros(1,length(v));
+    z1 = repmat(c,1,length(v));
+end
 vertices = [[x0;y0;z0],[x0;y0;z1]]';
 vertices(7,:)=[];
 vertices(end,:)=[];
 
-if draw
+if app.draw
     figure;
     line(x0,y0,z0);
     line(x0,y0,z1);
-    
     hold on;
-    plot3(x0,y0,z0,'r*')
-    plot3(x0(1),y0(1),z0(1),'b*')
-    plot3(x0,y0,z1,'r*')
-    plot3(x0(1),y0(1),z1(1),'b*')
+    if app.vertices
+        plot3(x0,y0,z0,'.k')
+        plot3(x0(1),y0(1),z0(1),'.b')
+        plot3(x0,y0,z1,'.k')
+        plot3(x0(1),y0(1),z1(1),'.b')
+    end
     for ii=1:length(x0)
         plot3([x0(ii) x0(ii)],[y0(ii) y0(ii)],[z0(ii) z1(ii)], 'b-', 'LineWidth', 1)
     end
-    % lab system : x y z vectors
-    quiver3(-1,0,0,2,0,0,'r');
-    quiver3(0,-1,0,0,2,0,'k');
-    quiver3(0,0,-1,0,0,2,'b');
-    axis equal;
-    xlabel('x (100)')
-    ylabel('y (010)')
-    zlabel('z (001)')
+    quiver3(0,0,z0(1),1,0,0,0,'k')
+    quiver3(0,0,z0(1),-0.5,0.866,0,0,'k')
+    quiver3(0,0,z0(1),-0.5,-0.866,0,0,'k')
+    quiver3(0,0,z0(1),0,0,1+0.4,0,'k')
+    text(1+0.1,0,z0(1),'a1','Color',[0 0 0])
+    text(-0.5-0.1,0.866+0.1,z0(1),'a2','Color',[0 0 0])
+    text(-0.5-0.1,-0.866-0.1,z0(1),'a3','Color',[0 0 0])
+    text(0,0,z1(1)+0.5,'c','Color',[0 0 0])
+    if app.numbers
+        for ii=1:12
+            text(vertices(ii,1),vertices(ii,2),vertices(ii,3)+0.1,num2str(ii),'Color',[0 0 0])
+        end
+    end
+    view(103,20)
+    if app.axes
+        % lab system : x y z vectors
+        quiver3(-1,0,0,2,0,0,0,'r');
+        quiver3(0,-1,0,0,2,0,0,'k');
+        quiver3(0,0,-1,0,0,2,0,'b');
+        xlabel('x')
+        ylabel('y')
+        zlabel('z')
+    end
+end
+
+% hexagon face directions
+%                   +yz      -yz        -xz       +xz     +xy      -xy
+%                   +x       -x         -y        +y       +z       -z
+faces_sides = [1 2 8 7; 4 5 11 10; 5 6 12 11; 2 3 9 8; 3 4 10 9; 6 1 7 12];
+faces_end   = [1 2 3 4 5 6; 7 8 9 10 11 12]; % -
+hkil        = [1 0 -1 0; -1 0 1 0; 0 -1 1 0; 0 1 -1 0; -1 1 0 0; 1 -1 0 0; ...
+               0 0 0 -1; 0 0 0 1]'; % column vectors
+         
+if app.draw && app.patch
+    if app.ind >=1 && app.ind <=6
+      patch('Vertices', vertices, 'Faces', faces_sides(app.ind,:), ...
+        'FaceVertexCData', repmat([0 0.8 0],length(app.ind),1), 'FaceColor', 'flat');
+    elseif app.ind==7 || app.ind==8
+      patch('Vertices', vertices, 'Faces', faces_end(app.ind-6,:), ...
+        'FaceVertexCData', repmat([0 0.8 0],length(app.ind-6),1), 'FaceColor', 'flat');
+    elseif app.ind==0
+      patch('Vertices', vertices, 'Faces', faces_sides, ...
+        'FaceVertexCData', repmat([0 0.8 0],6,1), 'FaceColor', 'flat');
+      patch('Vertices', vertices, 'Faces', faces_end, ...
+        'FaceVertexCData', repmat([0 0.8 0],2,1), 'FaceColor', 'flat');
+    end
 end
-% 
-% 
-% % hexagon made of hexagons
-% figure;
-% axis equal;
-% for f=[1,2] % we test <f>!
-%     for a=0:60:360
-%         line(f*cx*cosd(a)+cv,f*cx*sind(a)+sv);
-%     end
-%     if f==2
-%         for a=30:60:360
-%             line(c2*cx*cosd(a)+cv,c2*cx*sind(a)+sv);
-%         end
-%     end
-% end
 
-% 
-% xy=unique([x,y],'rows')
-% 
-% x=xy(:,1);
-% y=xy(:,2);
-% xx=[];
-% yy=[];
-% for i=1:3
-%     p=convhull(x,y);
-%     xx=[xx;x(p)]; % quick and dirty!
-%     yy=[yy;y(p)];
-%     x(p)=[];
-%     y(p)=[];
-% end
-% at=atan2(xx,yy);
-% [ax,ax]=sort(at);
-% xx=xx(ax);
-% yy=yy(ax);
-% line(xx,yy,...
-%     'marker','s',...
-%     'color',[0,0,0]);
-% % test some data points
-% px=20*rand(500,1)-10;
-% py=20*rand(500,1)-10;
-% line(px,py,...
-%     'marker','.',...
-%     'markersize',4,...
-%     'linestyle','none');
-% ix=inpolygon(px,py,xx,yy);
-% line(px(ix),py(ix),...
-%     'marker','.',...
-%     'linestyle','none',...
-%     'color',[1,0,0]);
 
+
+vertices(end+1,:) = [0 0 z0(1)];
+vertices(end+1,:) = [0 0 z1(1)];
+
+vertices(end+1,:) = [0 0 (z1(1)-z0(1))/3];
+vertices(end+1,:) = [0 0 (z1(1)-z0(1))/3*2];
+
+vertices(end+1,:) = [0 0 (z1(1)-z0(1))/4];
+vertices(end+1,:) = [0 0 (z1(1)-z0(1))/4*2];
+vertices(end+1,:) = [0 0 (z1(1)-z0(1))/4*3];
+
+if app.draw
+  axis equal
+  set(gca,'Visible','off')
+    if app.numbers
+        for ii=13:14
+            text(vertices(ii,1),vertices(ii,2),vertices(ii,3)+0.1,num2str(ii),'Color',[0 0 0])
+        end
+    end
+    if app.vertices
+        for ii=13:length(vertices)
+            plot3(vertices(ii,1),vertices(ii,2),vertices(ii,3),'.k')
+        end
+    end
 end
+if nargout == 1
+    data.vertices = vertices;
+    data.faces_sides = faces_sides;
+    data.faces_end = faces_end;
+    data.hkil = hkil;
+    data.hf  = gcf;
+end
+
+end % end of function
diff --git a/6_rendering/gtPlotHexagon.m b/6_rendering/gtPlotHexagon.m
index afea72947095cc7ea86ec5b766aecfdfc6f9a444..cc9636e110100a2ad7fc67bbbd0ad14a59e5bb86 100644
--- a/6_rendering/gtPlotHexagon.m
+++ b/6_rendering/gtPlotHexagon.m
@@ -1,45 +1,43 @@
 function hf = gtPlotHexagon(grains, r_vectors, angle, cmap)
 %     hf = gtPlotHexagon(grains, r_vectors, angle, cmap)
 %     --------------------------------------------------
-%     make a figure containing subplots of the unit cells of the selected grains
-%     angle is the angle of the section plane, defined as in
-%     gtShowSampleSurface
-%     can pass in a colormap as well, so colours match the section
+%     Makes a figure containing subplots of the unit cells of the selected grains
 %
-%     as in gtShowSampleSurface
-%     Deal with coordinate systems.  R-vector orientation data in the instrument
+%     INPUT:
+%       grains    = IDs of grains of interest <double 1xN>
+%       r_vectors = orientations for all the grains <double Nx4>
+%       angle     = rotation around z-axis to visualize orientations from the
+%                   section plane as defined in gtShowSampleSurface
+%       cmap      = colour map for all the grains (rgb) <double N+1x3>
+%
+%     OUTPUT:
+%       hf = figure handle
+%
+%     Deals with coordinate systems.  R-vector orientation data in the instrument
 %     system: x//beam, z upwards, y towards door.
 %     In the reconstructed tomography/grain map data, z is downwards, y//x
 %     instruments, x//y instrument.
-%     Finally present all information in the frame of the tomo data.
+%
+%     Reference hexagon coordinates see gtHexagonalUnitCell
 
 
-if ~exist('cmap', 'var')
+if ~exist('angle','var') || isempty(angle)
+    angle = 0;
+end
+if ~exist('cmap', 'var') || isempty(cmap)
     %if no colourmap passed in, use red
-    cmap=zeros(size(r_vectors, 1), 3);
-    cmap(:,1)=1;
+    cmap = zeros(size(r_vectors, 1), 3);
+    cmap(:,1) = 1;
 end
 
 
 % patch object hexagonal prism
-% 12 corners
-% [u v t w]
-corners4=[ 2 -1 -1 3;...
-          -1  2 -1 3;...
-          -1 -1  2 3;...
-           1  1 -2 3;...
-           1 -2  1 3;...
-          -2  1  1 3];
-corners4=[corners4; corners4];
-corners4(7:12, 4)=-3;%lower corners
+data = gtHexagonalUnitCell('ratio',2,'draw',false);
+corners3 = data.vertices(1:12,:);
+rec_faces = data.faces_sides;
+hex_faces = data.faces_end;
 
-% turn these into 3 index coords
-corners3(:,3) = corners4(:,4);
-corners3(:,1) = (corners4(:,1)-corners4(:,3))*cosd(30);
-corners3(:,2) = corners4(:,2)-((corners4(:,1)+corners4(:,3))*sind(30));
-corners3 = corners3./repmat(sqrt(sum((corners3.*corners3),2)), 1,3);
-
-hf=figure;
+hf=figure();
 
 % apply the R-vector
 for i=1:length(grains)
@@ -51,28 +49,17 @@ for i=1:length(grains)
 
     % allow for the two coordinate systems.  Transform from instrument to
     % reconstructed tomo coordinates: x-->y y-->x z-->-z
-    %corners3a=[corners3a(:,2) corners3a(:,1) -corners3a(:,3)]; 
+%     corners3a=[corners3a(:,2) corners3a(:,1) -corners3a(:,3)]; 
     
     % do the rotate to follow the section plane
     corners3b = ([cosd(angle) sind(angle) 0; -sind(angle) cosd(angle) 0;0 0 1]*corners3a')';
 
-    % write these corners into the correct format for patch
-
-    % corners3a is "vertices"
-    hex_faces=[1 4 2 6 3 5; 7 11 9 12 8 10];
-    rec_faces=[1 7 10 4;...
-        4 2 8 10;...
-        2 6 12 8;...
-        6 3 9 12;...
-        3 5 11 9;...
-        5 1 7 11];
-
     subplot(1, length(grains), i);
     % reference hex
-    % patch('Vertices', corners3, 'Faces', hex_faces, 'FaceColor', 'b', 'FaceAlpha', 0.4)
-    % patch('Vertices', corners3, 'Faces', rec_faces, 'FaceColor', 'b', 'FaceAlpha', 0.4)
+    patch('Vertices', corners3, 'Faces', hex_faces, 'FaceColor', 'b', 'FaceAlpha', 0.4)
+    patch('Vertices', corners3, 'Faces', rec_faces, 'FaceColor', 'b', 'FaceAlpha', 0.4)
     % reference background like the section plane
-    xpos=min(mean(corners3b([1 4 2 6 3 5], 1)), mean(corners3b([7 11 9 12 8 10], 1)));
+    xpos=min(mean(corners3b(hex_faces(1,:), 1)), mean(corners3b(hex_faces(2,:), 1)));
     section_plane_vertices=[xpos 1 1; xpos 1 -1; xpos -1 1; xpos -1 -1];
     section_plane_face=[1 2 4 3];
     patch('Vertices', section_plane_vertices, 'Faces', section_plane_face, 'FaceColor', 'k', 'FaceAlpha', 0.4)
@@ -84,19 +71,13 @@ for i=1:length(grains)
     set(get(get(hf, 'CurrentAxes'), 'xLabel'),'String', 'x axis')
     set(get(get(hf, 'CurrentAxes'), 'yLabel'),'String', 'y axis')
     set(get(get(hf, 'CurrentAxes'), 'zLabel'),'String', 'z axis')
-    %set(get(hf, 'CurrentAxes'), 'zdir','reverse')
-    %set(get(hf, 'CurrentAxes'), 'xdir','reverse')
     % need to mirror twice to stay consistant with tomo
     set(get(hf, 'CurrentAxes'), 'view', [90 0])
 
     axis tight
     axis equal
+    axis square
     axis vis3d
 end
 
-
-
-
-
-
-
+end % end of function
diff --git a/zUtil_Drawing/gtDrawGrainUnitCells.m b/zUtil_Drawing/gtDrawGrainUnitCells.m
index 248deafb7884230664c58f2ed4c3b926310e1f4b..1d0c88eac76316ba82e16768e20dd834753336ec 100644
--- a/zUtil_Drawing/gtDrawGrainUnitCells.m
+++ b/zUtil_Drawing/gtDrawGrainUnitCells.m
@@ -124,14 +124,11 @@ elseif strcmpi(app.type, 'hexagonal')
     c = app.latticepar(3);
     maxx = max(a,c);
     % normalize the biggest axis to 1
-    vertices = gtHexagonalUnitCell(a/maxx,c/maxx,false);
-    % hexagon face directions
-    %                   +yz      -yz        -xz       +xz                          +xy                -xy
-    %                   +x       -x         -y        +y                           +z                 -z
-    faces_sides = [1 2 8 7; 4 5 11 10; 5 6 12 11; 2 3 9 8; 3 4 10 9; 6 1 7 12];
-    faces_end   = [1 2 3 4 5 6; 7 8 9 10 11 12];
-    hkl         = [1 0 -1 0; -1 0 -1 0; 0 -1 1 0; 0 1 -1 0; -1 1 0 0; 1 -1 0 0; ...
-                   0 0 0 1; 0 0 0 -1]'; % column vectors
+    data = gtHexagonalUnitCell('ratio',c/a,'draw',false);
+    vertices = data.vertices(1:12,:);
+    faces_sides = data.faces_sides;
+    faces_end = data.faces_end;
+    hkl = data.hkil;
 end
 
 [Bmat, Amat] = gtCrystHKL2CartesianMatrix(app.latticepar);