%**************************************************************************
%**************************************************************************

%                          _                               
%                         (_)                              
%          ____  ____ ____ _ ____                          
%         / _  |/ ___) _  | |  _ \                         
%        ( ( | | |  ( ( | | | | | |                        
%         \_|| |_|   \_||_|_|_| |_|                        
%        (_____|                                           
%    _                           _             _            
%   | |                         | |           (_)           
%   | | _   ___  _   _ ____   _ | | ____  ____ _  ____  ___ 
%   | || \ / _ \| | | |  _ \ / || |/ _  |/ ___) |/ _  )/___)
%   | |_) ) |_| | |_| | | | ( (_| ( ( | | |   | ( (/ /|___ |
%   |____/ \___/ \____|_| |_|\____|\_||_|_|   |_|\____|___/ 
%                                                         
%                                    _                        
%                                   (_)                       
%         ____ _   _ ____  ____ ____ _  ___  ____             
%       / _  | \ / ) _  )/ ___) ___) |/___)/ _  )            
%      ( (/ / ) X ( (/ /| |  ( (___| |___ ( (/ /             
%       \____|_/ \_)____)_|   \____)_(___/ \____)            
%                                                         
%           ______   ______ ______   __                      
%          (_____ \ / __   (_____ \ /  |                     
%            ____) ) | //| | ____) )_/ |                     
%           /_____/| |// | |/_____/  | |                     
%           _______|  /__| |_______  | |                     
%          (_______)\_____/(_______) |_|                     
%                                                         
%                                                         
%                                                                                                           
%**************************************************************************
%**************************************************************************
    

% Exercise 5: grain boundaries and misorientations
% we will load a crop from a larger map, previously saved in mtex
dpath = '/home/rk/mtex/workshop/workshop21/exercises/'
load([dpath 'exp_qtz_crop.mat'])
% let's set the symmetry to the Laue group
ebsd('q').CS =  ebsd('q').CS.Laue;
ebsd('q').CSList{2}.color=[0.2 0.1 0.8] % this fix is just needed bc the 
% data was saved under an older version of mtex

%% 1) Compute grains 
% and determine the total grain boundary length. How much does it change
% when smoothign grains?

[grains ebsd.grainId ebsd.mis2mean] = ebsd.calcGrains('angle',[10 1]*degree);
l1 = sum(grains.boundary('q','q').segLength)

grains = smooth(grains,7)
l2 = sum(grains.boundary('q','q').segLength)

% change
(l1-l2)/l1

%% 2) Make a nice plot of the grain boundary directions.

% all qtz-qtz boundaries
gbq = grains.boundary('q','q')

gbdir = gbq.direction

plot(gbq,mod(gbdir.rho,pi)/degree,'linewidth',2)
mtexColorMap hsv % hsv isn't so nice but does the job,
% a nice colormap for directions is 'C2' from
% https://www.peterkovesi.com/matlabfns/Colourmaps/colorcet.m

% What is the preferred direction?

% one option
vm= mean(gbdir);
vm.rho/degree


% another option
omega= 0:0.1:180;
cumpl = surfor(gbq,omega*degree);
plot(omega,cumpl)

[~,idmax]=max(cumpl);
[~,idmin]=min(cumpl);
omega(idmax)
omega(idmin);



%% 3) Determine low angle (LA) boundaries
% since we computed grains with a second argument (1 degree), we habe LAB down to 1
% degree, however they go up to... high angles

plot(ebsd('q'),ebsd('q').orientations)
hold on 
plot(gbq)
hold on
plot(grains('q').innerBoundary,'linecolor','b') %innerBoundary 
hold off

% max angle of inner boundaries
max(grains('q').innerBoundary.misorientation.angle)/degree

% condition for misorientation angles of inenr boundaries smaller ...
cond = grains('q').innerBoundary.misorientation.angle < 9*degree;

% those are the low angle boundaries now
LAB = grains('q').innerBoundary(cond)

plot(ebsd('q'),ebsd('q').orientations,'FaceAlpha',0.3)
hold on 
plot(LAB,'linecolor','r')
hold off


% this is the misorientation of the low angle boundaries
mori = LAB.misorientation;

%--------------------------------------------------------------------------
% a) Determine the distribution of low angle boundaries (1-9 degree).
plotAngleDistribution(mori)



%--------------------------------------------------------------------------
% b) Colorcode misorientation axes in crystal coordiantes and plot them in
% a map.


% misorientation axes in crystal coordinates
axC =  mori.axis;
% colorcoding
cKC = HSVDirectionKey(ebsd('q').CS)


figure
plot(ebsd('q'),ebsd('q').orientations,'FaceAlpha',0.3)
hold on
plot(LAB,cKC.direction2color(axC),'linewidth',2)
hold off
figure
plot(cKC)



%--------------------------------------------------------------------------
% c) Colorcode misorientation axes in specimen coordiantes and plot them in
%    a map.

% extract orientations from both sides of the LAB
o = ebsd('id',LAB.ebsdId).orientations

% compute the misorientation axis in specimen coordinates
axS = axis(o(:,1),o(:,2),'antipodal')
% colorcoding
cKS= HSVDirectionKey(specimenSymmetry('-1'))


figure
plot(ebsd('q'),ebsd('q').orientations,'FaceAlpha',0.3)
hold on
plot(LAB,cKS.direction2color(axS),'linewidth',2)
hold off
figure
plot(cKS)



%--------------------------------------------------------------------------
% d) Plot the misorientation axes in crystal coordiantes in a point plot
%    and colrocode them according to the corresponding direction in
%    specimen coordinates.
%    AND
%    plot the misorientation axes in specimen coordinates and color code
%    with respect to the misorientatin axis in crytsal coordinates.
%
plot(axC,cKS.direction2color(axS),'MarkerSize',4,'fundamentalRegion')
nextAxis
plot(cKS)
nextAxis
plot(axS,cKC.direction2color(axC))
nextAxis
plot(cKC)


%% Special
% since the above may be confusing, let us select just a small region
xy = [14.6 77.5 1.5 1.9];% 
ebsd_sel =  ebsd(inpolygon(ebsd,xy));
id_sel = inpolygon(LAB,xy);

% make a sharp colorkey
aack = axisAngleColorKey
aack.oriRef=mean(ebsd_sel('q').orientations);
colorm=aack.orientation2color(ebsd_sel('q').orientations);

% map w boundary in crystal coords.
plot(ebsd_sel('q'),colorm,'FaceAlpha',0.9)
hold on
plot(LAB(id_sel),'linecolor','k','linewidth',8)
hold on
plot(LAB(id_sel),cKC.direction2color(axC(id_sel)),'linewidth',3)
hold off
mtexTitle('LAB crystal coords')
nextAxis

% map w boundary in specimen coords.
plot(ebsd_sel('q'),colorm,'FaceAlpha',0.9)
hold on
plot(LAB(id_sel),'linecolor','k','linewidth',8)
hold on
plot(LAB(id_sel),cKS.direction2color(axS(id_sel)),'linewidth',3)
hold off
mtexTitle('LAB specim coords')
nextAxis

% axes in crystal coords.
plot(cKC)
hold on
plot(axC(id_sel).project2FundamentalRegion,'fundamentalRegion', 'MarkerSize',4)
hold off
mtexTitle('LAB crystal coords')
nextAxis

% axes in specimen coords.
plot(cKS)
hold on
plot(axS(id_sel), 'MarkerSize',4)
mtexTitle('LAB specim coords')

% add boundary traces
% cluster the boundary direction - we just do this sicne there are two
% boundaries here
gbdir=LAB(id_sel).direction;
cid=calcCluster(gbdir);
mb1 = mean(gbdir(cid ==1));
mb2 = mean(gbdir(cid ==2));
hold on
plot(cross([mb1 mb2],zvector),'plane','linestyle','--','linewidth',2)
hold off
nextAxis

% orientations as polefigure
plotPDF(ebsd_sel('q').orientations,colorm,Miller({0,0,0,1},{1,1,-2,0},ebsd('q').CS), ...
        'MarkerSize',4)


%% 4) Select low angle (LA) boundaries
% a) Select LA boundaries with a misorientation angle of 3 - 5 degree.

LAB(mori.angle>3*degree & mori.angle<5*degree) 

%--------------------------------------------------------------------------
% b) Select all LA boundaries with misorientation axes which are parallel 
%    to [1,0,-1,0] (within some threshold).

threshold = 5*degree
refDirC = Miller(1,0,-1,0,ebsd('q').CS);
cond = angle(axC,refDirC) < threshold;
LAB(cond);

plot(LAB,'linewidth',2)
hold on
plot(LAB(cond),'LineColor','r','linewidth',2)
hold off
nextAxis

plot(axC(cond).project2FundamentalRegion, ...
     'MarkerFaceColor','r','fundamentalRegion')
hold on
plot(axC(~cond).project2FundamentalRegion, ...
     'MarkerFaceColor','k','fundamentalRegion')
hold on
plot(refDirC.symmetrise('unique'),'MarkerFaceColor','b')
circle(refDirC.symmetrise('unique'),threshold,'linewidth',2,'LineColor','b')
hold off

%--------------------------------------------------------------------------
% c) Select all LA boundaries with misorientation axes which are parallel to 
%    the z direction. Display with a proper color coding.

threshold = 15*degree;
refDirS =zvector;
cond = angle(axS,refDirS) < threshold;
LAB(cond);


plot(LAB,'linewidth',2)
hold on
plot(LAB(cond),'LineColor','r','linewidth',2)
hold off
nextAxis


plot(axS(cond),'MarkerFaceColor','r')
hold on
plot(axS(~cond),'MarkerFaceColor','k')
hold on
plot(refDirS,'antipodal','MarkerFaceColor','b')
circle(refDirS,threshold,'linewidth',2,'antipodal','LineColor','b')
hold off


%--------------------------------------------------------------------------
% d) Plot all boundaries which have a c-axis misrotation angle > 3*degree and 
%    <9degree. Plot the corresponding rotation axes.

% derive caxis direction
hc = Miller(0,0,0,1,ebsd('q').CS);
cdir = o.*hc

% rotation axis of c-axis
ctilt = cross(cdir(:,1),cdir(:,2));
% rotation angle of c-axis
ctilta = angle(cdir(:,1),cdir(:,2));

cond = ctilta > 3*degree & ctilta < 9*degree;

plot(ctilt(cond),ctilta(cond)/degree)

% plot those rotation axes in a map

% color key for caxis rotation axis
sk  = HSVDirectionKey(specimenSymmetry('-1'));

plot(ebsd('q'),sk.direction2color(ebsd('q').orientations.*hc))
hold on
plot(LAB(cond),'lineColor','k','linewidth',4)
hold on
plot(LAB(cond),sk.direction2color(ctilt(cond)),'linewidth',2)
hold off
nextAxis

plot(ctilt(cond),sk.direction2color(ctilt(cond)),'antipodal','MarkerSize',4)

% this looks actuially quite funny. Can someone explain this?

%% 5) Can we check if some of these LABs are tilt boundaries
%  A tilt boundary should have the rotation axis within the boundary plane.
%  What is possible with 2D EBSD?

% compute the angle between the grain boundary trace and the misorientation
% axis in specimen coordiantes. Only if those are "parallel" we can be sure
% for this given boudnary that it might be a tilt boundary. Since we lack
% any other 3D information on the grain boudnary, any assumption for tilt
% boundaries would involve some additional infor on it's inclination.

gbta = angle(axS,LAB.direction);

plot(LAB,gbta/degree,'linewidth',2)
hold on
quiver(LAB(1:5:end),axS(1:5:end),'color','k')
hold off



%% Bonus exercise gb / misorientation
%**************************************************************************
% Question (not really boundary related):
% Some people like to click somewhere in a grain and then get a
% map with the points color-coded according to the misorientation angle 
% with respect to this reference point.
%
% Write a little script which does exactly provide this functionality.
% Bonus: Do not only display a map with the misorientation angle but also
% with the misorientation axes for each point in the grain.
%**************************************************************************
dpath = '/home/rk/mtex/workshop/workshop21/exercises/'
load([dpath 'exp_qtz_crop.mat'])
ebsd('q').CS =  ebsd('q').CS.Laue;

% calc grains
[grains, ebsd.grainId]=calcGrains(ebsd,'angle',10*degree);
% plot it with default ipdf parallel Z (no OM specified)
plot(ebsd('q'),ebsd('q').orientations)

%select a point
xy=ginput(1);
% get just the grains of the point
g_sel = grains(findByLocation(grains,xy));
e_sel = ebsd(g_sel)

% get the orientation at point xy
oatpoint=e_sel(findByLocation(e_sel,xy)).orientations;
% calc the misorientation to point 
mTP= inv(e_sel.orientations).*oatpoint;

% plot angles
plot(e_sel,mTP.angle/degree)
hold on
plot(g_sel.boundary)
hold on
plot(xy(1),xy(2),'Marker','p','MarkerSize',10,'MarkerFaceColor','r')
hold off
mtexColorbar('Title','mis. angle2point')

% plot axis
% axis in specimen coordinates
axS = axis(e_sel.orientations, oatpoint,'antipodal');

dmS=HSVDirectionKey(specimenSymmetry('-1'));
% axis in crystal coordinates
axC = mTP.axis; 
dmC=HSVDirectionKey(axC.CS);

figure
plot(e_sel,dmC.direction2color(axC))
hold on
plot(g_sel.boundary)
hold on
plot(xy(1),xy(2),'Marker','p','MarkerSize',10,'MarkerFaceColor','r')
hold off
nextAxis

plot(e_sel,dmS.direction2color(axS))
hold on
plot(g_sel.boundary)
hold on
plot(xy(1),xy(2),'Marker','p','MarkerSize',10,'MarkerFaceColor','r')
hold off

% plot colorcodings
figure
plot(dmC)
hold on
plot(axC,'points',500,'MarkerSize',4)
hold off

nextAxis
plot(dmS)
hold on
plot(axS,'points',500,'MarkerSize',4)
hold off


