%==========================================================================
%==========================================================================
%
%                                    /$$                                    
%                                   | $$                                    
%                    /$$$$$$/$$$$  /$$$$$$    /$$$$$$  /$$   /$$            
%                   | $$_  $$_  $$|_  $$_/   /$$__  $$|  $$ /$$/            
%                   | $$ \ $$ \ $$  | $$    | $$$$$$$$ \  $$$$/             
%                   | $$ | $$ | $$  | $$ /$$| $$_____/  >$$  $$             
%                   | $$ | $$ | $$  |  $$$$/|  $$$$$$$ /$$/\  $$            
%                   |__/ |__/ |__/   \___/   \_______/|__/  \__/            
%                                                                           
%                                                                           
%                                                                           
%                              /$$                       /$$                
%                             | $$                      | $$                
%                     /$$$$$$ | $$$$$$$   /$$$$$$$  /$$$$$$$                
%                    /$$__  $$| $$__  $$ /$$_____/ /$$__  $$                
%                   | $$$$$$$$| $$  \ $$|  $$$$$$ | $$  | $$                
%                   | $$_____/| $$  | $$ \____  $$| $$  | $$                
%                   |  $$$$$$$| $$$$$$$/ /$$$$$$$/|  $$$$$$$                
%                    \_______/|_______/ |_______/  \_______/                
%                                                                           
%                                                                           
%                                                                           
%                                                    /$$                    
%                                                   |__/                    
%   /$$$$$$  /$$   /$$  /$$$$$$   /$$$$$$   /$$$$$$$ /$$  /$$$$$$$  /$$$$$$ 
%  /$$__  $$|  $$ /$$/ /$$__  $$ /$$__  $$ /$$_____/| $$ /$$_____/ /$$__  $$
% | $$$$$$$$ \  $$$$/ | $$$$$$$$| $$  \__/| $$      | $$|  $$$$$$ | $$$$$$$$
% | $$_____/  >$$  $$ | $$_____/| $$      | $$      | $$ \____  $$| $$_____/
% |  $$$$$$$ /$$/\  $$|  $$$$$$$| $$      |  $$$$$$$| $$ /$$$$$$$/|  $$$$$$$
%  \_______/|__/  \__/ \_______/|__/       \_______/|__/|_______/  \_______/
%                                                                           
%                                                                           
%==========================================================================
%==========================================================================
 
% The aim of this exercise is to get familiar with mtex look at some ebsd
% data

% Notes:
% - usually you can find help by typing `help` command
%   if a command with an identical name might be available for different
%   data types try e.g. `help orientation/calcMDF` or `help odf/calcMDF` or
%   help orientation/mean` (`help mean` will be Matlabs mean for numbers)
% - using TAB-completion makes many tasks faster, prevents typos and helps
%   to explore functions and properties available to an object
%   e.g. ebsd.[hit the tabulator key] will present a list of features
%  available with the ebsd object
%% install/startup mtex (in case it isn't already)
addpath /Path/to/mtex 
startup_mtex

%% 0) load some data -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% we can use the import_wizard and generate a script which imports the ebsd
% data


import_wizard
% selcte EBSD and click Next ...
% when asked for the reference frame, x should be east and y south when
% using convertEuler2SpatialReferenceFrame


%% 1) explore the ebsd data *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% hint: using ebsd.[TAB] will show you some properties and functions

% %########################################################################
% Question: What can you tell about the ebsd data?
% %########################################################################

ebsd
% how many indexed pixels ,how many non-indexed?


ebsd.prop

%       mad = mean angular deviation:misfit between detected bands and
%       simulation (in degree, high mad = 'bad' solution)
%       bc = band contrast: ratio of actual intensity in hough space vs
%       average intensity of kikuchi bands (high value = high contrast)
%       bs =band slope : higest intensity gradient at the morgin of kikuchi
%       bands (high value = sharp band)

% access the first 1000 pixels

% access the  first 1000 indexed pixels

%% % 2) look at the map -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

% use the plot command to plot almost anything in mtex
plot(ebsd)
% maybe not so spectacular at first sight, it will just plot the phase

% mtex will always honor the convention
% plot(where,what)
% where in the case of EBSD data, the "where" are x,y coordiantes of each
% measurement point and the "what" can be anything as logn as it has the 
% same size/length as the first argument.

% %########################################################################
% Question:
% a) Colorcode the band contrast in a grey-scale color-coding
% and plot non-indexed pixels in red.
% hint1: when plotting (where,what) values, both need to have the same
% length
% hint2: to overlay two plots, use 'hold on' and 'hold off'
% %########################################################################

plot(ebsd,ebsd.prop.bc)
mtexColorMap white2black
mtexColorbar('Title','band contrast')

% xlabel(c,'bands') % on older matlab versions
hold on  % this keeps the present figure axis
plot(ebsd('notIndexed'),'FaceColor',[1 0.2 0],'FaceAlpha',0.6)
hold off 


%% 2) inspect orientations -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

% a) access orientations and the first 100 Euler angles
ebsd('q').orientations
ebsd('q').orientations(1:100).Euler


% b) plot a polefigure (point plot)
%    - define crystal directions (@Miller)
%    - either use plotPDF or remember that orientiation*Miller -> specimen direction

cs = ebsd('q').CS
hs= [... % 
 Miller(0,0,0,1,cs),... %c
 Miller(1,1,-2,0,cs),... %a
 Miller(1,0,-1,0,cs),... %m
 Miller(0,1,-1,1,cs),... %z
 Miller(1,0,-1,1,cs)... %r
 ];
plotPDF(ebsd('q').orientations,hs)


% c) plot an inverse pole figure (point plot)
%    - define a specimen direction (@vector3d)
%    - either use plotIPDF or remember that orientation*vector3d -> crystal direction    
refv = [zvector vector3d(-5,1,0) vector3d(2,-3,0)]
plotIPDF(ebsd('q').orientations,refv,'MarkerSize',2)


% d) plot an orientation mapping  
%    - Is one mapping sufficient to understand the texture? 

oM = ipfHSVKey(ebsd('q').CS);
oM.inversePoleFigureDirection = refv(1);
col=oM.orientation2color(ebsd('q').orientations);
plot(ebsd('q'),col)

% bonus:
% e) colorcode points in a polefigure and in an inverse polefigure indentical
%    to the colorcoding in the map
for i=1:length(refv)
    oM.inversePoleFigureDirection = refv(i);
    col=oM.orientation2color(ebsd('q').orientations);
    figure
    plotPDF(ebsd('q').orientations,col, hs([1 2 4]))
    nextAxis
    plotIPDF(ebsd('q').orientations, col, refv,'MarkerSize',2)
    nextAxis 
    plot(ebsd('q'),col)
end


% f) plot a c-axis mapping
%    - the c-axis is just a direction, use .Miller2color

% h) plot orientations in orientations space


%% 3) Further questions ---------------------------------------------------
% %########################################################################
% a) Question:
% Do pixels indexed with just few bands (e.g. <=6) show a different 
% polefigure compared to those with more bands (e.g. >9) ?
% %########################################################################

ebsdQ=ebsd('q')
ind6=ebsdQ.prop.bands <= 6;
ind9=ebsdQ.prop.bands > 9;

figure
plotPDF(ebsdQ(ind6).orientations,ebsdQ(ind6).prop.mad,hs,'MarkerSize',3)
mtexColorbar('Title','mad')
CLim(gcm,[min(ebsd('q').prop.mad) max(ebsd('q').prop.mad)])
figure
plotPDF(ebsdQ(ind9).orientations,ebsdQ(ind9).prop.mad,hs,'MarkerSize',3)
mtexColorbar('Title','mad')
CLim(gcm,[min(ebsd('q').prop.mad) max(ebsd('q').prop.mad)])

figure
plotPDF(ebsdQ(ind6).orientations,hs,'MarkerSize',3)
f=gcm;
f.currentId = 1
hold on
plotPDF(ebsdQ(ind9).orientations,hs,'MarkerSize',3)
hold off

%%
% %########################################################################
% b) Question:
% Compare the <1,1,-2,0> distributions of those orientation which contain 
% c-axes at the periphery of the pole figure to those with c-axes in the center 
% %########################################################################

cs =   ebsdQ.CS;
cdir = ebsdQ.orientations.*Miller(0,0,0,1,cs);
cdir.antipodal =1;
cond1 = angle(cdir,zvector,'antipodal')<45*degree;
a=symmetrise(Miller(1,1,-2,0,cs));
aF1  =calcDensity(ebsdQ(cond1).orientations*a,'halfwidth',10*degree)
aF2  =calcDensity(ebsdQ(~cond1).orientations*a,'halfwidth',10*degree)

plot(aF1,'minmax','TR',{'pfJ:', num2str(round(10*mean(aF1.^2))/10)})
nextAxis
plot(aF2,'minmax','TR',{'pfJ:', num2str(round(10*mean(aF2.^2))/10)})

%%
% %########################################################################
% c) Question:
%  Plot a fibre mapping , color all orientations with <a> up to 25 degree
%  away from the a-maximum     
% %########################################################################
 
aF  =calcDensity(ebsdQ.orientations*a,'halfwidth',10*degree)
[~,maxpos] = max(aF);

om =ipfSpotKey(cs),
om.center=a.unique;
om.inversePoleFigureDirection=maxpos;
om.psi=deLaValeePoussinKernel('halfwidth',15*degree);

col = om.orientation2color(ebsdQ.orientations);
plot(ebsdQ,col)
figure
plotPDF(ebsdQ.orientations,col,qtzMiller(cs))
annotate(maxpos,'add2all')


%% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% 4) orientation/direction/angle mapping 
% -------------------------------------------------------------------------
% a) plot ebsd data in a map with ipdf colorcoding wrt to Z
% and a direction paralle to E-W (X)  (z in our example is the kinematic/strain y-direction,
% and the flow plane is oriented E-W
% plot also a pole figure of c,a,-a,m,r,z where individual points are colorcoded
% with the ipdf color-coding


%% -------------------------------------------------------------------------
% b) plot ebsd data color-coded with respect to a fibre e.g. with 0,0,0,1 
% parallel to the mean c-axis direction in a range of 20 degree
% hint: get the mean direction of the c-axis and use it as inversePoleFigureDirection
% in the orientation mapping
cs=ebsd('q').CS;
h=Miller(0,0,0,1,cs);

c_dir = ebsd('q').orientations.*h;
oM = ipfSpotKey(cs);
oM.inversePoleFigureDirection = mean(c_dir,'antipodal'); %specimendirection
oM.center = h; % fibre
oM.color = [1 0 1];
conewidth = 20*degree;
oM.psi = deLaValeePoussinKernel('halfwidth',conewidth);

figure;
plot(oM)
hold on
circle(oM.center,2*conewidth,'linewidth',2)
hold off
figure;
plot(ebsd('q'),oM.orientation2color(ebsd('q').orientations))

% how many of all orientations are in there?
fibreVolume(ebsd('q').orientations,h,mean(c_dir,'antipodal'),10*degree)
%% -------------------------------------------------------------------------
% c) plot the direction of the c-axis in a map in a "nice" color coding
% hint: the c-axis is just a direction, use .oM.Miller2color
% oMs = ipfHSVKey(crystalSymmetry('-1'))
% you will be informed that this color-coding is not topologically continuous

oMs = HSVDirectionKey;
plot(oMs,'zAxisDirection','intoPlane','xAxisDirection','east')
c_dir = ebsd('q').orientations.*Miller(0,0,0,1,ebsd('q').CS)
c_dir.antipodal = 1;

figure
ind= discretesample(1:length(ebsd('q')),500);
plot(c_dir(ind),oMs.direction2color(c_dir(ind)),'antipodal','MarkerSize',5)
% plot(c_dir(ind),oMs.Miller2Color(c_dir(ind)),'antipodal','MarkerSize',5)

figure
plot(ebsd('q'),oMs.direction2color(c_dir))
% plot(ebsd('q'),oMs.Miller2Color(c_dir))

%% -------------------------------------------------------------------------
% %########################################################################
% Question:
% d) plot the angle between flow pane (NW-SW) and an a-axis and plot only those
% within a range < 30 degree
% hint: first we need to get all eq a-axes and we can loop through all
% angles between the a-axes and the reference direction and find the minaml
% value
% %########################################################################
% hint: rotate the vector into crystal coordinates and calc the angle with the
% ref. crystal direction, inv(orientation*vector=Miller (inverse polefigure equation)

%%
 close all
 %% 5) manually select data - ebsd area
%select an ebsd area by clicking around it, then hit enter
plot(ebsd('q'),ebsd('q').orientations)

% old method
% ebsd_xy = ginput; % click around a polygon
% ebsd_sel = ebsd(inpolygon(ebsd,ebsd_xy))

% new methode
ebsd_sel = selectInteractive(ebsd)

plot(ebsd_sel('q'),ebsd_sel('q').orientations)
%%
 close all
%% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% 
%              /$$$$$$$$ /$$$$$$$   /$$$$$$  /$$$$$$$                       
%             | $$_____/| $$__  $$ /$$__  $$| $$__  $$                      
%             | $$      | $$  \ $$| $$  \__/| $$  \ $$                      
%             | $$$$$   | $$$$$$$ |  $$$$$$ | $$  | $$                      
%             | $$__/   | $$__  $$ \____  $$| $$  | $$                      
%             | $$      | $$  \ $$ /$$  \ $$| $$  | $$                      
%             | $$$$$$$$| $$$$$$$/|  $$$$$$/| $$$$$$$/                      
%             |________/|_______/  \______/ |_______/         

%                             /$$                                           
%                            | $$                                           
%                          /$$$$$$$$                                        
%                         |__  $$__/                                        
%                            | $$                                           
%                            |__/          

%         /$$$$$$  /$$$$$$$   /$$$$$$  /$$$$$$ /$$   /$$  /$$$$$$           
%        /$$__  $$| $$__  $$ /$$__  $$|_  $$_/| $$$ | $$ /$$__  $$          
%       | $$  \__/| $$  \ $$| $$  \ $$  | $$  | $$$$| $$| $$  \__/          
%       | $$ /$$$$| $$$$$$$/| $$$$$$$$  | $$  | $$ $$ $$|  $$$$$$           
%       | $$|_  $$| $$__  $$| $$__  $$  | $$  | $$  $$$$ \____  $$          
%       | $$  \ $$| $$  \ $$| $$  | $$  | $$  | $$\  $$$ /$$  \ $$          
%       |  $$$$$$/| $$  | $$| $$  | $$ /$$$$$$| $$ \  $$|  $$$$$$/          
%        \______/ |__/  |__/|__/  |__/|______/|__/  \__/ \______/       

%                                                    /$$                    
%                                                   |__/                    
%   /$$$$$$  /$$   /$$  /$$$$$$   /$$$$$$   /$$$$$$$ /$$  /$$$$$$$  /$$$$$$ 
%  /$$__  $$|  $$ /$$/ /$$__  $$ /$$__  $$ /$$_____/| $$ /$$_____/ /$$__  $$
% | $$$$$$$$ \  $$$$/ | $$$$$$$$| $$  \__/| $$      | $$|  $$$$$$ | $$$$$$$$
% | $$_____/  >$$  $$ | $$_____/| $$      | $$      | $$ \____  $$| $$_____/
% |  $$$$$$$ /$$/\  $$|  $$$$$$$| $$      |  $$$$$$$| $$ /$$$$$$$/|  $$$$$$$
%  \_______/|__/  \__/ \_______/|__/       \_______/|__/|_______/  \_______/
%                                                                           
%                                                                           
%                                                                           
%                                                                      
% The aim of this exercise is to become more familiar with ebsd data and
% grains, access and plot properties of grains

                                                       
%% 5) grains -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

% make grains inspect grains and access grains properties
[grains,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd,'angle',10*degree);
% grains are the grains
% grainId tells ev1ery ebsd pixel to which grain it belongs
% mis2mean tells each ebsd pixel it's misorientation to the mean grain orientation

% it is  also possible to use any other criterion to segment grains
% e.g. only the c-axis directions - not that it'll be particularly
% usefull but maybe intereting to compare with optical data
caxgrains = calcGrains(ebsd('q'),'custom',ebsd('q').orientations*Miller(0,0,0,1,cs),'delta',6*degree)

%% -------------------------------------------------------------------------
% a) explore the properties of grains:
% plot the frequency distribution of the grain orientation spread (.GOS) and
% colorcode each grains acording to the GOS
hist(grains('q').GOS./degree,20)
figure
plot(grains('q'),grains('q').GOS./degree)
mtexColorbar('Title','GOS')

%% -------------------------------------------------------------------------
% b) plot a histogram of the area equivalent diameter of grains, what is
% the frequency weighted mean/mode grain size?
figure
hist(grains('q').equivalentRadius.*2,20);
mean(grains('q').equivalentRadius.*2);

% to get an vague approximation of area (volume) fractions one can use
% frequncy.*bincenter.*r^2 (or ^3)
[binpop,bincenter] = hist(grains('q').equivalentRadius.*2,20);
vol = binpop.*bincenter.^2;
bar(bincenter,[binpop/sum(binpop); vol./sum(vol)]');

% So far we did not consider twinning, so keep those numbers in mind
% and go on to see if we were wrong with our gs estimate

%% ------------------------------------------------------------------------
% %########################################################################
% c) grains have boundaries: and we can access those with grains.boundary
% Question: How long is the total qtz-qtz boundary in mm
% hint: inspect grains.boundary('q','q').[TAB]
% %########################################################################

% note these are stepped boundaries so maybe not so close to true
% comarpe with smooth boundaries grains_smoothed=smooth(grains('q'),5,'exp');

%% -------------------------------------------------------------------------
% Before we proceed, we should make sure that we are looking at grains and not twins:
% A very common twin in quartz is the Dauphin?? twin, which is defined by
% a 60 degree totation around 0001

% it becomes pretty obvious when plotting e.g. the misorientation angles.
% axis-angle pairs etc. more on that later this week. For now we will just
% look at the angle histogram of grain boundary misorientations
% option 1)
hist(grains('q').boundary('q','q').misorientation.angle./degree,1:110)
% option 2)
plotAngleDistribution(grains('q').boundary('q','q').misorientation, 'resolution',1*degree)

%% -------------------------------------------------------------------------
% e) get rid of twins and merge grains which are connected about a boundary
% which obeys a twinnign relation

% option 1) merge the grains
% get the grain boundaries
gb_qtz = grains.boundary('q','q');
 %define a misorientation assosciated with dauphine twinning
rot = orientation('axis',Miller(0,0,0,1,ebsd('q').CS),'angle',60*degree,cs,cs);
% calculates the angle between two trotations - the misorientation and the twinning rotation 
ang_r = angle(gb_qtz('q','q').misorientation,rot);
% condition/index all angles below a certain (magic) threshold
ind = ang_r < 6.*degree; % what angle should we choose here?
%select the boundaries based on this condition                                                                        
twinBoundary = gb_qtz(ind);
% merge grains connected with that twin boundary
[grains_merged,parentID_grains]= merge(grains,twinBoundary);
% make a copy of the ebsd dataset
ebsd_merged = ebsd;
% update the grainIds to the parentIds, so each ebsd_merged pixel knows
% again to which grain it belongs to
ebsd_merged.grainId  = parentID_grains(ebsd.grainId);
% plot it
figure
plot(ebsd('q'),ebsd('q').orientations,'faceAlpha',0.15)
hold on
plot(grains_merged.boundary('q','q'),'linewidth',2)
hold on
plot(twinBoundary,'linecolor','r','linewidth',2)

% note that now, grains do not carry any meanOrientation - since it is
% probably nonsense!

% just for this particular twinning relation on can also temporarily
% increase the crystal symmetry to hexagonal - the twin dioes not exist - 
% segment grains and set back the symmetry. Now the meanOrientation esists
% and represents something like the modal orientation of the grains

% option 2) change the symmetry
csq= ebsd('q').CS
cshq = crystalSymmetry('622','mineral', 'Quartz')
ebsd('q').CS=cshq
[grainsD,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd,'angle',10*degree);
ebsd('q').CS=csq
grainsD('q').CS=csq
hold on
plot(grainsD.boundary('q','q'),'linecolor','b','linewidth',2)
hold off

% from here, go on with grains_merged / ebsd_merged or grainsD/ebsd

plotAngleDistribution(grainsD('q').boundary('q','q').misorientation, 'resolution',1*degree)


% compare grainsD and grains_merged
grainsD
grains_merged
% Question: Do they have the same number of grains and if not why and what
% could we do about it?

%%
close all
%% ------------------------------------------------------------------------
% e) plot "real" equivalent diameter in a map and in a histogram
plot(grainsD('q'),2.*grainsD('q').equivalentRadius)

bins=0:3:30;
[gs,~] = hist(grains('q').equivalentRadius.*2,bins)
[gsD,~] = hist(grainsD('q').equivalentRadius.*2,bins)

figure
bar(bins', [gs' gsD'])
legend(['no Dauphine twins considered: dmean =' num2str(mean(grains('q').equivalentRadius.*2))], ...
       ['Dauphine twins considered: dmean =' num2str(mean(grainsD('q').equivalentRadius.*2))])

%% -------------------------------------------------------------------------
% f) Explore grain parameters:
% Question: plot a map of grains which have an aspect ratio larger 2.5

% for simplicity create a variable just containing quartz grains
% and use interpolated shapes and discard grain which are too small
% note: with serious, work-on data, you probably wouldn't want to trust
% shapes of grains smaller 10,20,50 ... pixels, depending on the property

gmqtz=grainsD(grainsD.grainSize> 5, 'q')
gmqtz = smooth(gmqtz,5,'exp')


%look at the aspect ratios
plot(gmqtz,gmqtz.aspectRatio)
mtexColorbar('Title','aspect ratio')
CLim(gcm,[1 10])

% define a condition
condA = gmqtz.aspectRatio > 2.5;
plot(gmqtz(condA),gmqtz(condA).aspectRatio)
mtexColorbar('Title','aspect ratio')


%% -------------------------------------------------------------------------
% ########################################################################
% Question:
% plot only grains with an aspect ratio larger 2.5 AND with a long axis
% direction +/- 10 degree from ESE-WNW, plot ipdf standard colorcoding
% hint: use principalComponents to get the parameters for best fit ellipses
[theta,a,b]= gmqtz.principalComponents('hull'); %
theta =  mod(theta./degree,180);
% %########################################################################

% plot all laxdir with a reasonable colorcoding
plot(gmqtz,theta)
mtexColorMap(hsv)
mtexColorbar
condB = theta > 20 & theta < 40;

plot(gmqtz(condB & condA))


figure
plot(gmqtz,'FaceColor',[0.7 0.7 0.7])
hold on
plot(gmqtz(condA),'FaceColor',[1 0 0],'FaceAlpha',0.2,'noBoundary')
hold on
plot(gmqtz(condB),'FaceColor',[0 0 1],'FaceAlpha',0.2,'noBoundary')
%hold on
%plot(gmqtz(condA & condB).boundary,'linewidth',2)
hold off


%% -------------------------------------------------------------------------
% %########################################################################
% Question:
% compare the texture of grains with high aspect ratios and long axes +/-
% parallel to ESE-WNW to grains with lower aspect ratios. Also compare 
% number weighted (1 point per grains) and area weighted (1 point per pixel)
% measures of texture strength.
% hint: to access the ebsd of grains, you can index ebsd by grains, for
% example ebsd(grains([1 7 8])) returns all the ebsd points contrained in
% grains 1,7 and 8.

o1 = ebsd(gmqtz(condB & condA)).orientations
o2 = ebsd(gmqtz(gmqtz.aspectRatio<2)).orientations

og1 = gmqtz(condB & condA).meanOrientation
og2 = gmqtz(gmqtz.aspectRatio<1.5).meanOrientation

% we could compare odf based measures e.g. texture index
psi1=calcKernel(og1,'method','magicRule')
psi2=calcKernel(og2,'method','magicRule')

odf1=calcODF(o1,'kernel',psi1)
odf2=calcODF(o2,'kernel',psi2)

textureindex(odf1)
textureindex(odf2)

plotPDF(og1,hs,'MarkerFaceColor',[1 0 0],'MarkerSize',3,'')
f=gcm;
f.currentId = 1
hold on
plotPDF(og2,hs,'MarkerFaceColor',[0 0 1],'MarkerSize',3)
hold off
f.drawNow

plotPDF(og1,hs,'smooth')
f=gcm;
f.currentId = 1
hold on
plotPDF(og1,hs,'contour','linewidth',2,'linecolor',[0 0 0])
hold off
f.drawNow


%...
% %########################################################################
%% 
% some other things that  may come handy
%% -------------------------------------------------------------------------
% i)  manually select data - grains
%select a grain by clicking at it, when done, hit enter
plot(grainsD('q'),grainsD('q').meanOrientation)
grain_xy= ginput;% close on those grains you want to select, hit ENTER when done
selected_grains = grainsD(findByLocation(grainsD,grain_xy))
plot(ebsd(selected_grains),ebsd(selected_grains).orientations)

%%
 close all
%% 6) misorientations-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

% have a look at grans.boundary.misorientations.[TAB]

% a) colorcode grain boundaries according to the misorientation angle on
% top of a semi-tranparent map showing grain meanorientation
% 'FaceAlpha' 0-1 should take case of transparency
plot(grainsD('q'),grainsD('q').meanOrientation,'FaceAlpha',0.2)
hold on
plot(grainsD('q').boundary('q','q'),grainsD('q').boundary('q','q').misorientation.angle./degree,'linewidth',2)
hold off
c = mtexColorbar;ylabel(c,'misorientation angle')
%% -------------------------------------------------------------------------
% b) plot and colorcode the misorientation axes along boundaries with respect 
% to crystal ref. frame. Those will be of class Miller, so you can use an
% orinetation mapping with Miller2Color
oMC = TSLOrientationMapping(grainsD('q').CS.Laue);
plot(grainsD('q').boundary('q','q'),oMC.Miller2color(grainsD('q').boundary('q','q').misorientation.axis), 'linewidth',2)
figure
plot(oMC)
%% -------------------------------------------------------------------------
% %########################################################################
% Question:
% c) plot a map with only those boundaries which have misorientation axes 
% close to 10-10
% hint: calc misor angles and axes and use anf angle condition for logical indexing
% %########################################################################


%% -------------------------------------------------------------------------
% d) plot and colorcode the misorientation axes along boundaries with
% respect to sample ref. frame. Those will be of class vector3d
% get orientations along boundary
% hint: you can get the orientation along the boundary from
% ebsd(grains.boundary.ebsdID).orientations
% see also help orientation/axis
% note: misorientation is based in crystal coordinates, and
% misorientation.axis will be in crystalcords.
% orientations are based in specimen coordinates, hence axis(o1,o2) will be
% based in specimen coordinates

% get all orientations from ebsd along boundary
gbo = ebsd(grainsD.boundary('q','q').ebsdId,'q').orientations;    
% get the misorientation axis
misaxS = axis(gbo(:,1),gbo(:,2));
% what is the symmetry of misaxS
oMS = TSLOrientationMapping(crystalSymmetry('-1')); 
    plot(oMS,'zAxisDirection','intoPlane','xAxisDirection','east')
    hold on
    plot(misaxS,'antipodal','contour','contours',0:1:3)
    hold off
figure    
plot(grainsD('q').boundary('q','q'),oMS.Miller2color(misaxS), 'linewidth',2)
%% -------------------------------------------------------------------------
% %########################################################################
% Question:
% e) plot those misorientation axes in a pf/ipf with misorientation angles between 10 and
% 20degree
% hint: calc misor angles and axes and use anf angle condition for logical indexing
% %########################################################################


%%
close all
%% -------------------------------------------------------------------------
% f) comapre c-axis "misorientation" and real misorientation angles

% plot the angles
h=Miller(0,0,0,1,ebsd('q').CS)
gb= grainsD('q').boundary('q','q');
gbo = ebsd(gb.ebsdId,'q').orientations;  

c_axis_angle =   angle(gbo(:,1).*h,gbo(:,2).*h,'antipodal');

plot (c_axis_angle./degree,...
     gb.misorientation.angle./degree, 'o','MarkerSize',1)
xlabel('c-axis angle'); ylabel('misor angle');

% plot the c-axis angle
figure
plot(ebsd('q'),ebsd('q').orientations,'FaceAlpha',0.1)
hold on
plot(gb,c_axis_angle./degree,'linewidth',2)
hold off
c = mtexColorbar;
ylabel(c,'c axis "misor" angle')


% plot the difference betwen c-axis misorientation angles and real
% misorientation angles 
figure
plot(ebsd('q'),ebsd('q').orientations,'FaceAlpha',0.1)
hold on
plot(gb,(gb.misorientation.angle - c_axis_angle)./degree ,'linewidth',2)
hold off
c = mtexColorbar;
ylabel(c,'diff angle')

% but maybe more interesting;
% plot the misorientation angles for those boundaries 
% with c-axis "misorientations" smaller e.g. 3degree
% this can be used when comaring ebsd based segmentations with optically
% derived segmentations

lim=3
ind=c_axis_angle<lim*degree;

figure
plot(ebsd('q'),ebsd('q').orientations,'FaceAlpha',0.1)
hold on
plot(gb(ind),gb(ind).misorientation.angle./degree,'linewidth',2)
hold off
c = mtexColorbar;
ylabel(c,['true misor angle (for c-ax misor angle < ' num2str(lim) ')'])

%%
close all

%% %% -------------------------------------------------------------------------
% g) manually select data along lines -> misor profiles
%select ebsd along a line
plot(ebsd('q'),ebsd('q').orientations)
lineX = ginput(2) % select a line change (2) if your profile should have more than two points
[ebsd_lineX, dist ]= spatialProfile(ebsd('q'),lineX);
%plot pdf with points colorcoded after distance to origin
plotPDF(ebsd_lineX('q').orientations,dist, hs,'points','all','MarkerSize',5)

% ebsd along the line can be used for (mis)orientation profiles, eg
% misorinetation angles
o=ebsd_lineX('q').orientations
m = o(1:end-1).\o(2:end) %misor to neigbour
mo = o(1).\o(2:end) %misor to origin

%
figure
plot(dist(1:end-1),angle(mo)/degree,'linewidth',2)
hold on
plot(dist(1:end-1), angle(m)/degree,'color','r','linewidth',2)
xlabel('distance','FontSize',16); ylabel('misorientation angle in degree','FontSize',16)
legend('to reference orientation','to neighbour')
title(['misor angles (all angles)'],'FontSize',16)
%%
close all
%%
% -------------------------------------------------------------------------
% %########################################################################

% h) Big-Exercise: Compare misorientations of low angle boundaries (2-10 degree)
% of two types of grains, those with a flat (peripheral) and those with a steep
% (central) c-axis orientation, plot misorientation axes in sample and in 
% crystal reference frame, plot maps of grain boundaries, with a colorcoding
% for misorientation axes in crystal and sample reference frame

% Note: if you want to access boundaries with smaller misorientation angles, 
% you have to segment grains with a smaller threshold and select relevant angles above noise
% ebsd_s= ebsd
% [subgrains,ebsd_S.grainId,ebsd_s.mis2mean] = calcGrains(ebsd,'angle',0.01*degree);
%
% Suggestion: 
% 1) get caxis directio of all grains
% 2) discrimniate two groups of grains based on angle between c-axis direction and a vector
% 3) create 2 ebsd subsets from those grains
% 4) calculate two sets of new_grains with a very small segmentation angle
% 5) caluclate misorientation angle, axes etc from both new subsets of new_grains
% 6) use the misorientation angle to index reasonable values e.g. 2-10
% degree
% 7) plot grains of both c-axes subsets, plot boundary misorientations of new_grains
% ...

% %########################################################################




% 