%==========================================================================
%==========================================================================
%__                        
%                                                |  \                       
%                     ______    ______   ______   \$$ _______               
%                    /      \  /      \ |      \ |  \|       \              
%                   |  $$$$$$\|  $$$$$$\ \$$$$$$\| $$| $$$$$$$\             
%                   | $$  | $$| $$   \$$/      $$| $$| $$  | $$             
%                   | $$__| $$| $$     |  $$$$$$$| $$| $$  | $$             
%                    \$$    $$| $$      \$$    $$| $$| $$  | $$             
%                    _\$$$$$$$ \$$       \$$$$$$$ \$$ \$$   \$$             
%                   |  \__| $$                                              
%                    \$$    $$                                              
%                     \$$$$$$                                               
%                                                    __                     
%                                                   |  \                    
%  ______   __    __   ______    ______    _______  \$$  _______   ______  
% /      \ |  \  /  \ /      \  /      \  /       \|  \ /       \ /      \ 
%|  $$$$$$\ \$$\/  $$|  $$$$$$\|  $$$$$$\|  $$$$$$$| $$|  $$$$$$$|  $$$$$$\
%| $$    $$  >$$  $$ | $$    $$| $$   \$$| $$      | $$ \$$    \ | $$    $$
%| $$$$$$$$ /  $$$$\ | $$$$$$$$| $$      | $$_____ | $$ _\$$$$$$\| $$$$$$$$
% \$$     \|  $$ \$$\ \$$     \| $$       \$$     \| $$|       $$ \$$     \
%  \$$$$$$$ \$$   \$$  \$$$$$$$ \$$        \$$$$$$$ \$$ \$$$$$$$   \$$$$$$$
%                                                                           
%                                                                           
%                                                                           
%                     ______    ______    ______     __                     
%                    /      \  /      \  /      \  _/  \                    
%                   |  $$$$$$\|  $$$$$$\|  $$$$$$\|   $$                    
%                    \$$__| $$| $$$\| $$ \$$__| $$ \$$$$                    
%                    /      $$| $$$$\ $$ /      $$  | $$                    
%                   |  $$$$$$ | $$\$$\$$|  $$$$$$   | $$                    
%                   | $$_____ | $$_\$$$$| $$_____  _| $$_                   
%                   | $$     \ \$$  \$$$| $$     \|   $$ \                  
%                    \$$$$$$$$  \$$$$$$  \$$$$$$$$ \$$$$$$                  
%                                                                           
%
%==========================================================================
%==========================================================================
%% import a simple ebsd dataset using the import_wizard
% use sample40.ang
import_wizard('EBSD')

%% we can put the output of import_wizard here

% crystal symmetry
CS = {... 
  'notIndexed',...
  crystalSymmetry('-3m', [4.9 4.9 5.4], 'X||a', 'Y||b*', 'Z||c*', ...
                 'mineral', 'Quartz', 'color', [0.53 0.81 0.98])};

% plotting convention
setMTEXpref('xAxisDirection','east');
setMTEXpref('zAxisDirection','intoPlane');

% path to files
pname = '/home/rk/mtex/workshop/workshop21/exercises/';

% data from Giulia degli Alessandrini MSc thesis,
% "THE EITA SHEAR ZONE; Val Grosina (N Italy), 2013
% supervised by Neil Mancktelow, Eric Reusser
% data by Karsten Kunze, ETHZ

%fname = [pname '/sample13.ang'];
fname = [pname '/sample40.ang'];

% create an EBSD variable containing the data
ebsd = EBSD.load(fname,CS,'interface','ang',...
  'convertEuler2SpatialReferenceFrame','setting 2');

% note that when we use the import_wizard the "setting" option form OIM is 
% not set, meaning 'setting 1' by default (for historic reasons). However, 
% 'setting 2' seems to be the default in many cases as is also the case here.
%% 1) compute grains and do an ipf colorcoding
% useful command: calcGrains, ipfHSVKey
[grains, ebsd.grainId] = ebsd.calcGrains('angle',10*degree);

plot(grains)

cs = ebsd('q').CS;
ck = ipfHSVKey(cs);
ck.inversePoleFigureDirection = xvector;

plot(grains('q'),ck.orientation2color(grains('q').meanOrientation))

%% 2) colorcode the c-axis direction of grains and compute grains 
% based on c-axis direction
% useful command: calcGrains, HSVDirectionKey

ck2 = HSVDirectionKey(specimenSymmetry('-1'));

cdir = grains('q').meanOrientation.* Miller(0,0,0,1,cs);
color = ck2.direction2color(cdir);

nextAxis
plot(grains('q'),color)

% Do you observe a difference to the IPF mapping? Which one

% compute c-axes directions for each point
cdir_ebsd = ebsd('q').orientations.*Miller(0,0,0,1,cs);

% in case of notIndexed, assign nan vectors to notIdnexed pixels
v = vector3d.nan(length(ebsd),1);
v(ebsd.isIndexed)=cdir_ebsd;
v.antipodal=1;
% compute grains using the custom option

cgrains = ebsd.calcGrains(ebsd,'custom',v,'delta',10*degree)
nextAxis
plot(cgrains('q'),ck2.direction2color(cgrains('q').meanOrientation*Miller(0,0,0,1,cs)))


%% 3) overcome a twin related problem
% if we inspect the orientation maps (c-axis and ipf colorcoding) 
% we observe many mottled grains which relate to a twinning related problem
% -> quartz has a popular twinning law which is a 60 degree rotation
% around the c-axis (Dauphiné law).

% How can we get around this and compute grains with a meaningful mean
% orientation which ignore those twins WITHOUT using merge (because the 
% solution is more elegant and we do not know about boundaries yet)

% What symmetry operation would correspond to the Dauphiné twinning?

% useful hint: we can change the crystal symmetry of ebsd and grain
% objects but we must take care to use unique names


cshex = crystalSymmetry('622','mineral','hexq');
ebsd('q').CS = cshex;

[grains, ebsd.grainId ebsd.mis2mean] = ebsd.calcGrains('angle',8*degree);

ebsd('hexq').CS = cs;
grains('hexq').CS = cs;

nextAxis
plot(grains('q'),ck.orientation2color(grains('q').meanOrientation))

% 3b) delete all grains <=5 pixels (let's asume they are bad measurements)

ebsd(grains(grains.grainSize<=5))=[];
% repeat step above
ebsd('q').CS = cshex;
[grains, ebsd.grainId ebsd.mis2mean] = ebsd.calcGrains('angle',8*degree);
ebsd('hexq').CS = cs;
grains('hexq').CS = cs;

%% 4) compute grain properties
% a) Plot a map of the GOS and map of grain size. Is there any correlation
% between those?

qgrains = grains('q')

plot(qgrains,qgrains.equivalentRadius.*2)
nextAxis
plot(qgrains,qgrains.GOS./degree)
mtexColorbar('title',{'eq_{dia}' 'GOS'})

% we could also compute the correlation coefficient
[R,P,RLO,RUP]=corrcoef(qgrains.area,qgrains.GOS)

%% 5) Is there a relation of texture and grain size?
% a) Create two groups, one with the 100 largest grains

gs = qgrains.equivalentRadius.*2;
[~,id] = sort(gs,'descend');
g1 = qgrains(id(1:100))
g2 = qgrains(id(101:end))

plot(g1,'FaceColor','Salmon','lineColor','k')
hold on
plot(g2,'FaceColor','CadetBlue','lineColor','k')
hold off

% b) plot polefigures for each class
h = Miller({0,0,0,1},{1,1,-2,0},{1,0,-1,0},{1,0,-1,1},{0,1,-1,1},ebsd('q').CS)
plotPDF(g1.meanOrientation,h([1,2,5]),'MarkerSize',3,'points',100)
nextAxis
plotPDF(g2.meanOrientation,h([1,2,5]),'MarkerSize',3,'points',100)

% c) estimate two odfs and plot PDFs from those
% note that the polulations have hugely different sizes
% useful commands: calcKernel, calcDensity
% Do the different options of kernel width estimators have an influence on
% your conclusion?

psi1 = calcKernel(g1.meanOrientation)
psi2 = calcKernel(g2.meanOrientation)

% psi1 = calcKernel(g1.meanOrientation,'method','magicrule')
% psi2 = calcKernel(g2.meanOrientation,'method','magicrule')

odf1= calcDensity(g1.meanOrientation,'kernel',psi1)
odf2= calcDensity(g2.meanOrientation,'kernel',psi2)

plotPDF(odf1,h([1,2,5]),'contourf','contours',1:10,'minmax')
nextAxis
plotPDF(odf2,h([1,2,5]),'contourf','contours',1:10,'minmax')
CLim(gcm,'equal')


%% 7) Do equiaxed, small grains have a different texture than 
% elongated, small grains
% Fit an ellipse to grains and inspect the results
% We should probably cut away grains at the boundaries.
% useful commands: isBoundary, principalComponents,smooth

gs = grains('q').equivalentRadius.*2;
[~,id] = sort(gs,'descend');
g2 = qgrains(id(101:end))
% get rif of the boudnary grains
g2(g2.isBoundary)=[];
plot(g2)

g2 = smooth(g2,5)
[omega, a,b] = principalComponents(g2);

plot(g2,a./b)
hold on
plotEllipse(g2.centroid,a,b,omega,'lineColor','r')
hold off

% some condition for equiaxed grains
cond = a./b < 1.3;

psi1 = calcKernel(g2(cond).meanOrientation,'method','magicrule')
psi2 = calcKernel(g2(~cond).meanOrientation,'method','magicrule')

odf1= calcDensity(g2(cond).meanOrientation,'kernel',psi1)
odf2= calcDensity(g2(~cond).meanOrientation,'kernel',psi2)

plotPDF(odf1,h([1,2,5]),'contourf','contours',1:10,'minmax')
nextAxis
plotPDF(odf2,h([1,2,5]),'contourf','contours',1:10,'minmax')
CLim(gcm,'equal')


%% 8) What can we say about the preferred shape fabric of small grains
% plot a length and area weighted histograms of grain long axes
% useful commands: longAxis, paror

lax = g2.longAxis;
subplot(1,2,1)
histogram(lax,'weights',lax.norm,180)
subplot(1,2,2)
histogram(lax,'weights',g2.area,180)

omega = [0:360];
cumpl = g2.paror(omega*degree)
plot(omega,cumpl)
xlim([0,180])



%% 9) Question session: What else would you want to do with grains?





