Overlay HCP or cubic unit cells on EBSD Map #mtexScript

UnitCellOverlay

%% Draw HCP or cubic unit cells overlaid on EBSD Map rotated by Euler angles
%rev 1.
%file assumes the following; -you select the crystal unit cell to draw
%-that the default for HCP is z//c axis, x//to [1,0,-1,0],
%that the default for cubic is faces aligned with axes
%no promises at all, but a pole figure plotting utility included to check your annotations (based on mtex and not my work)

%IF YOUR DATA HAS MORE THAN ONE PHASE you'll need to edit the plot(ebsd...line)

%Rev1; fixed problem with cases where maps plotted with z into plane.
%start by plotting the map
figure; plot(ebsd('Indexed'),ebsd('Indexed').orientations)
%click where you want the cells, hit 'enter' when done.
[x,y]=ginput;

%set scale of drawn unit cells and line width
ScaleFactor=1;
LineWdth=2;
%fill drawn faces 1=fill
FillFace=1;

%measurement location markers.  1=on
LocMarker=0;
FontSzMkr=15;

%measurement location numbers 1=on
LocNum=0;
FontSzNum=15;

%plot HCP(1) or Cubic(2)
UnitCellType=1;

hold on
% Part 1 Draw the unit cell

if UnitCellType==1
%since MTEX default is X along [10-10], picture hexagon with pointy end up
%(in y) length of vertex from center = 6 units.  Label vertexes A-F.  Top
%is upper vertexes, B is bottom vertexes

%clockwise from top.  Note that we use vector3d instead of standard arrays
%to allow use of euler rotation functions
vertexAT=vector3d(0,6,12);
vertexBT=vector3d(5.196,3,12);
vertexCT=vector3d(5.196,-3,12);
vertexDT=vector3d(0,-6,12);
vertexET=vector3d(-5.196,-3,12);
vertexFT=vector3d(-5.196,3,12);
vertexAB=vector3d(0,6,0);
vertexBB=vector3d(5.196,3,0);
vertexCB=vector3d(5.196,-3,0);
vertexDB=vector3d(0,-6,0);
vertexEB=vector3d(-5.196,-3,0);
vertexFB=vector3d(-5.196,3,0);
%combine in vertex array 
v=[vertexAT;vertexBT;vertexCT;vertexDT;vertexET;vertexFT;vertexAB;vertexBB;vertexCB;vertexDB;vertexEB;vertexFB];
%provide information on the order in which to connect faces of hex prism.
%Note, sides only.
f=[1 2 8 7;2 3 9 8;3 4 10 9;4 5 11 10;5 6 12 11;6 1 7 12];
else
  
    %clockwise from top.  Note that we use vector3d instead of standard arrays
%to allow use of euler rotation functions
vertexAT=vector3d(12,0,12);
vertexBT=vector3d(0,0,12);
vertexCT=vector3d(0,12,12);
vertexDT=vector3d(12,12,12);
vertexAB=vector3d(12,0,0);
vertexBB=vector3d(0,0,0);
vertexCB=vector3d(0,12,0);
vertexDB=vector3d(12,12,0);

%combine in vertex array 
v=[vertexAT;vertexBT;vertexCT;vertexDT;vertexAB;vertexBB;vertexCB;vertexDB];
%provide information on the order in which to connect faces of hex prism.
%Note, sides only.
f=[1 2 6 5;2 3 7 6;3 4 8 7;4 1 5 8];
end

for i=1:size(x,1)
%enter rotation to transform the unit cell
oriP=ebsd(x(i),y(i)).orientations;
o = rotation('Euler',oriP.phi1,oriP.Phi,oriP.phi2);

%combine in vertex array and rotate
v_T3d=o*v;

if UnitCellType==1
%convert rotated vector3d back to matrix for plotting
v_T=[v_T3d(1).x v_T3d(1).y v_T3d(1).z;v_T3d(2).x v_T3d(2).y v_T3d(2).z;v_T3d(3).x v_T3d(3).y v_T3d(3).z;v_T3d(4).x v_T3d(4).y v_T3d(4).z;v_T3d(5).x v_T3d(5).y v_T3d(5).z;v_T3d(6).x v_T3d(6).y v_T3d(6).z;v_T3d(7).x v_T3d(7).y v_T3d(7).z;v_T3d(8).x v_T3d(8).y v_T3d(8).z;v_T3d(9).x v_T3d(9).y v_T3d(9).z;v_T3d(10).x v_T3d(10).y v_T3d(10).z;v_T3d(11).x v_T3d(11).y v_T3d(11).z;v_T3d(12).x v_T3d(12).y v_T3d(12).z];
else
 v_T=[v_T3d(1).x v_T3d(1).y v_T3d(1).z;v_T3d(2).x v_T3d(2).y v_T3d(2).z;v_T3d(3).x v_T3d(3).y v_T3d(3).z;v_T3d(4).x v_T3d(4).y v_T3d(4).z;v_T3d(5).x v_T3d(5).y v_T3d(5).z;v_T3d(6).x v_T3d(6).y v_T3d(6).z;v_T3d(7).x v_T3d(7).y v_T3d(7).z;v_T3d(8).x v_T3d(8).y v_T3d(8).z];
end


%scale v_T if required
v_T=v_T*ScaleFactor;
%check which way axes are positive

%move v_T to appropriate location
v_T(:,1)=v_T(:,1)+x(i);
v_T(:,2)=v_T(:,2)+y(i);

%if z plot direction is into plane, Ensure v_T(z) is all <0 to ensure all 
%lines above ebsd map surface and visible
if strcmpi(getMTEXpref('zAxisDirection'),'intoPlane')
test1=max(v_T(:,3));
if test1>0
v_T(:,3)=v_T(:,3)-test1;
end

else %if z out of plane Ensure v_T(z) is all >0 to ensure all lines above 
    %ebsd map surface and visible
test1=min(v_T(:,3));
if test1<0
v_T(:,3)=v_T(:,3)-test1;
end
end

if FillFace==1
h1=patch('Faces',f,'Vertices',v_T,'FaceColor','red','LineWidth',LineWdth); % plotting the prism with face and vertex data

else %plot wireframe
h1=patch('Faces',f,'Vertices',v_T,'FaceAlpha',0,'LineWidth',LineWdth); 
end

if LocMarker==1
     text(x(i),y(i),15,'x','FontSize',FontSzMkr,'FontWeight','bold');
end

%annotate measurement location numbers 1=on
if LocNum==1
    strText=int2str(i);
    text(x(i),y(i),15,strText,'FontSize',FontSzNum,'FontWeight','bold');
end


end
hold off
%view([0 0 1])

%% verifying with pf plots

%note, this must match the symmetry of the phase you selected
cs=ebsd('Indexed').CS;
h=[Miller(0,0,1,cs), Miller(1,0,0,cs)];

for i=1:size(x,1)
    FigTitle=(['Measurement point ',num2str(i)]);
    oriP=ebsd(x(i),y(i)).orientations;
    %if z plot direction is into plane, plot lower hemisphere
if strcmpi(getMTEXpref('zAxisDirection'),'intoPlane')
    hfig=figure; plotPDF(oriP,h,'lower','grid','projection','eangle','MarkerSize',10);
else
    hfig=figure; plotPDF(oriP,h,'upper','grid','projection','eangle','MarkerSize',10);
end
    %note, you can't do this in one step bc mtex retitles when it plots
    %pole figs
    set(hfig,'Name',FigTitle,'NumberTitle','off');
end