%
%
%                  /$$$$$$$$ /$$$$$$$   /$$$$$$  /$$$$$$$                       
%                 | $$_____/| $$__  $$ /$$__  $$| $$__  $$                      
%                 | $$      | $$  \ $$| $$  \__/| $$  \ $$                      
%                 | $$$$$   | $$$$$$$ |  $$$$$$ | $$  | $$                      
%                 | $$__/   | $$__  $$ \____  $$| $$  | $$                      
%                 | $$      | $$  \ $$ /$$  \ $$| $$  | $$                      
%                 | $$$$$$$$| $$$$$$$/|  $$$$$$/| $$$$$$$/                      
%                 |________/|_______/  \______/ |_______/                       
% 
% 
% 
%                                                        /$$                    
%                                                       |__/                    
%       /$$$$$$  /$$   /$$  /$$$$$$   /$$$$$$   /$$$$$$$ /$$  /$$$$$$$  /$$$$$$ 
%      /$$__  $$|  $$ /$$/ /$$__  $$ /$$__  $$ /$$_____/| $$ /$$_____/ /$$__  $$
%     | $$$$$$$$ \  $$$$/ | $$$$$$$$| $$  \__/| $$      | $$|  $$$$$$ | $$$$$$$$
%     | $$_____/  >$$  $$ | $$_____/| $$      | $$      | $$ \____  $$| $$_____/
%     |  $$$$$$$ /$$/\  $$|  $$$$$$$| $$      |  $$$$$$$| $$ /$$$$$$$/|  $$$$$$$
%      \_______/|__/  \__/ \_______/|__/       \_______/|__/|_______/  \_______/
% 
% 
% The aim of this exercise is to become more familiar with ebsd data using
% mtex and addressing and plotting properties.

% Quite often, 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`


%% 0) load some data -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% change directory to the place where the sample data is located
% cd workshopfolder

% we will use a cropped part of a larger ebsd map, hence at the moment we
% will not use the import_wizard for now
load ebsd_workshop.mat

% the crystal symmetry as it would be produced by the import_wizard
% template
% CS = {... 
%   'notIndexed',...
%   crystalSymmetry('-3m1', [4.913 4.913 5.501], 'X||a*', 'Y||b', 'Z||c', 'mineral', 'Quartz', 'color', 'purple')};
% this dataset originates from an hkl/oxford system, hence X||a*,
% originally they set the crystalSymmetry to be -3m and the c-axis as 5.5,
% however it should be 321 and 5.401, we will change this later
% % we can retrieve the crystal symmetry as stored in the dataset from the ebsd 
% CS = ebsd.CSList

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

%% 1) explore the ebsd data *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% hint: using ebsd.[TAB] will show you some properties and functions
% available for the ebsd object
% a) Which properties does it have, how big is it (in pixel, in micron), ...
        ext = ebsd.extend
        ext(2)-ext(1)
        ext(4)-ext(3)
        % which is identical to...
        max(ebsd.x)-min(ebsd.x)
        max(ebsd.y)-min(ebsd.y)
        % number of x,y pixels
        (max(ebsd.y)-min(ebsd.y))/(2*ebsd.unitCell(1,1))
        % we seem to a little rounding issue with this dataset
        
% -------------------------------------------------------------------------
% b) What is the cystal symmetry, which parameters are used?
        ebsd('indexed').CS
% Quartz is trigonal with a symmetry of '321', and cell parameters of a,b= 4.913 c= 5.405
% Oxford/HKL decides decided to use -3m, is that a proper group?   
         ebsd('q').CS.isProper   
% Can we change the symmetry to a proper group?       
        ebsd('q').CS=ebsd('q').CS.properGroup
% The c-axis should be 5.401, define the correct cs with a correct correct
% c-length
        cs=crystalSymmetry('321', [4.913 4.913 5.401],...
            'X||a*', 'Y||b', 'Z||c', 'mineral', 'Quartz', 'color', 'light blue')
        ebsd('q').CS=cs
% -------------------------------------------------------------------------        
% c) What can you say about the quality of the data?
%       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)
%   -  What is the mean/min/max quartz mad?
       mean(ebsd('q').prop.mad)
%   -  What is the percentage of pixels with a less reliable indexing?
        % query points with a specific condition e.g. % of pixels with mad > 1 
        sum(ebsd('indexed').prop.mad >1)/length(ebsd('indexed'))*100
%   -  Is there any relation between mad and band contrast?
        plot(ebsd('indexed').prop.mad,ebsd('indexed').prop.bc,'*','MarkerSize',2)
        ylabel('band contrast'); xlabel('mad')
%% % 2) plot properties in a map -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% a) Colorcode the number of bands used for indexing in a grey-scale color-coding
% and non-indexed pixels in red
% hint: plotting ebsd data in mtex works in almost all cases like plot(ebsd,property)
% while ebsd and property obviously need to have the same length
plot(ebsd('q'),ebsd('q').prop.bands)
mtexColorMap white2black
c = mtexColorbar
xlabel(c,'bands')
hold on
plot(ebsd('notIndexed'),'FaceColor','r')
hold off

% -------------------------------------------------------------------------
% b) plot the band contrast and plot band contrast lower 50 as black
% hint: you could either plot two maps on top of each other using hold
% on/off, you can set the colorrange using CLim and you can get the
% colormap cm=colormap of the current plot and maniplulate the values
% another hint: use logical indexing: e.g. a=1:9; a(a > 4) returns 5 6 7 8 9
% ebsd data can be querried logical indexing


plot(ebsd,ebsd.bc)
%CLim(gcm,[50  max(reshape(ebsd.bc,length(ebsd),1))]) % in case it is EBSDsquare
CLim(gcm,[50  max(ebsd.bc)]) % in the usual case
hold on
ind = ebsd.bc < 50;
plot(ebsd(ind),'FaceColor','k','DisplayName','bc < 50')
hold off
mtexColorbar

% another option in case min(ebsd.bc)==0
plot(ebsd,ebsd.bc)
cm = colormap;
id50=round(length(cm)/max(ebsd.bc)*50)
cm((1:id50),:)= repmat([0 0 0],length((1:id50)),1);
colormap(cm)
c = mtexColorbar; ylabel(c,'band contrast')
% -------------------------------------------------------------------------
% c) many types of data be colorcoded e.g.
% colorcode all quartz points in the dataset as a function of distance from
% an arbitrary line
% hint: distance from a line defined by two points, 
% e.g. http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html
% using [pq]= ginput(2) will let you select two points on a map
% this data can also be added to the ebsd dataset by 
% ebsd.prop.distancefromline =distancefromline

% 
plot(ebsd)
[pq]= ginput(2)
xy=[ebsd('q').prop.x ebsd('q').prop.y];
p=pq(1,:); q=pq(2,:);
d =(abs((q(2)-p(2)).*xy(:,1) - (q(1)-p(1)).*xy(:,2) + q(1).*p(2) -q(2).*p(1)))./sqrt((q(2)-p(2)).^2 + (q(1)-p(1)).^2);
plot(ebsd('q'),d)
hold on
line(pq(:,1),pq(:,2),'linewidth',2)
hold off
%

%% 3) orientations -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

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

% ------------------------------------------------------------------------- 
% b) display orientations as polefigure plot for c (0,0,0,1) and a
%  (1,1,-2,0), learn that orientation*Miller gives a vector in sample
%  reference which can be plotted using plot() (see also help vector3d/plot)
%  and also use plotPDF() (see help orientation/plotPDF)
% 

%define Millers
h = [Miller(0,0,0,1,cs,'uvw') Miller(1,1,-2,0,cs,'uvw')]
% option 1:  use plotPDF
c_dir= ebsd('q').orientations*h(1)
a_dirs=ebsd('q').orientations*symmetrise(h(2))
mf=mtexFigure;
plot(c_dir,'antipodal','MarkerSize',1,'parent',mf.nextAxis)
mtexTitle('0001')
plot(a_dirs,'antipodal','MarkerSize',1,'MarkerFaceColor','r','parent',mf.nextAxis)
mtexTitle('11-20')
% just plot 1000 points
ind= discretesample(1:length(ebsd('q')),1000);
plot(c_dir(ind),'antipodal','MarkerSize',2)
% option 2: use plotPDF
plotPDF(ebsd('q').orientations,h,'points',1000)

% -------------------------------------------------------------------------
% c) plot orientaions as sigma sections
% hint: have a look at help plotSection
plotSection(ebsd('q').orientations,'points',4000,'sections',6,'sigma')
% -------------------------------------------------------------------------
% d) plot orientations in an inverse pole figure wrt to x,z and a direction 
% oriented NW-SE 
% option 1: learn that orientation\vector or inv(o)*vector gives a Miller
% and option 2, use plotIDPF (see help orientation/plotIPDF)

dir=[xvector zvector vector3d(1,1,0)]
for i=1:length(dir)
m(:,i) =inv(ebsd('q').orientations).*dir(i);
end
m=project2FundamentalRegion(m);
m=Miller(m,ebsd('q').CS);
% to plot just 1000 points e.g.
sel=discretesample(1:length(ebsd('q')),1000);
mf=mtexFigure
for i=1:length(dir)
plot(m(sel,i),'FundamentalRegion','MarkerSize',2,'parent',mf.nextAxis)
end
% option 2: use plotIPDF
plotIPDF(ebsd('q').orientations,dir,'points',1000,'MarkerSize',2)
% -------------------------------------------------------------------------
% e) do pixels indexed with just few bands (e.g. <=6) show a different 
% pf than those with many (e.g. >9) ? Color-code each point according to
% the mad

ebsdQ=ebsd('q')
ind6=ebsdQ.prop.bands <= 6;
ind9=ebsdQ.prop.bands > 9;
figure
plotPDF(ebsdQ(ind6).orientations,ebsdQ(ind6).prop.mad,h,'MarkerSize',3)
c = mtexColorbar; ylabel(c,'mad')
figure
plotPDF(ebsdQ(ind9).orientations,ebsdQ(ind9).prop.mad,h,'MarkerSize',3)
c = mtexColorbar; ylabel(c,'mad')

figure
plotPDF(ebsdQ(ind6).orientations,h,'MarkerSize',3)
hold on
plotPDF(ebsdQ(ind9).orientations,h,'MarkerSize',3)
hold off

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

oM = ipdfHSVOrientationMapping(ebsd('q').CS);
oM.inversePoleFigureDirection = zvector;
plot(oM)
figure;
plot(ebsd('q'),oM.orientation2color(ebsd('q').orientations))
hs= [... % planes
 Miller(0,0,0,1,cs,'hkl'),... %c
 Miller(1,1,-2,0,cs,'hkl'),... %a
 Miller(-1,-1,2,0,cs,'hkl'),... %-a
 Miller(1,0,-1,0,cs,'hkl'),... %m
 Miller(0,1,-1,1,cs,'hkl'),... %z
 Miller(1,0,-1,1,cs,'hkl')... %r
 ];
figure
plotPDF(ebsd('q').orientations,oM.orientation2color(ebsd('q').orientations),hs)

oM.inversePoleFigureDirection = vector3d(1,1,0);
figure;
plot(ebsd('q'),oM.orientation2color(ebsd('q').orientations))
figure
plotPDF(ebsd('q').orientations,oM.orientation2color(ebsd('q').orientations),hs)

% -------------------------------------------------------------------------
% 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 from and use it as inversePoleFigureDirection
% in the orientation mapping
c_dir = ebsd('q').orientations.*Miller(0,0,0,1,ebsd('q').CS);
oM = ipdfCenterOrientationMapping(ebsd('q').CS);
oM.inversePoleFigureDirection = mean(c_dir,'antipodal'); %specimendirection
oM.center = Miller(0,0,0,1,ebsd('q').CS); % fibre
oM.color = [1 0 1];
conewidth = 10*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))

% -------------------------------------------------------------------------
% 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 = TSLOrientationMapping(xvector)
plot(oMs,'zAxisDirection','intoPlane','xAxisDirection','east')
% you will be informed that this color-coding is not topologically continuous
c_dir = ebsd('q').orientations.*Miller(0,0,0,1,ebsd('q').CS)

figure %hold on % ploting on top of the oM gave some problems on a particular system, freezing Matlab to death
ind= discretesample(1:length(ebsd('q')),500);
plot(c_dir(ind),oMs.Miller2color(c_dir(ind)),'antipodal','MarkerSize',5)

figure
plot(ebsd('q'),oMs.Miller2color(c_dir))

% -------------------------------------------------------------------------
% 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
a_all = symmetrise(Miller(1,1,-2,0,ebsd('q').CS))
aang = zeros(length(ebsd('q')),3);
for i=1:length(a_all)
    aang(:,i) = angle(ebsd('q').orientations.* a_all(i),vector3d(1,1,0),'antipodal')./degree;
end
aang = min(aang,[],2);
    
figure
plot(ebsd('q'),aang)
c=mtexColorbar; ylabel(c,'min. angle between flow plane and <a>')
%plot only those within a specific range
figure
ebsdq = ebsd('q');
plot(ebsdq(aang < 30),aang(aang < 30))
c=mtexColorbar; ylabel(c,'angles between <a> and flow plane < 30 degree')


%% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% 
%              /$$$$$$$$ /$$$$$$$   /$$$$$$  /$$$$$$$                       
%             | $$_____/| $$__  $$ /$$__  $$| $$__  $$                      
%             | $$      | $$  \ $$| $$  \__/| $$  \ $$                      
%             | $$$$$   | $$$$$$$ |  $$$$$$ | $$  | $$                      
%             | $$__/   | $$__  $$ \____  $$| $$  | $$                      
%             | $$      | $$  \ $$ /$$  \ $$| $$  | $$                      
%             | $$$$$$$$| $$$$$$$/|  $$$$$$/| $$$$$$$/                      
%             |________/|_______/  \______/ |_______/         

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

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

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

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

% let's use cleaned data
% (cleanup procedure was: grains were calculated, non-indexed areas <= 2px have been filled, grains 
% calculated again and dauphine twins were merged, qtz-grains of 1 px have been deleted and filled
% been deleted and filled)
load ebsd_clean_workshop.mat

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


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

% -------------------------------------------------------------------------
% 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)
c = mtexColorbar
ylabel(c,'GOS')

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

% 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: How long is the total qtz-qtz boundary in mm
% hint: inspect grains.boundary('q','q').[TAB]
sum(grains('q').boundary('q','q').segLength)/1000

% -------------------------------------------------------------------------
% d) plot the misorientation angle distribution
% hint: inspect grains.boundary('q','q').[TAB]
% and also consider using plotAngleDistribution
% option 1)
hist(grains('q').boundary('q','q').misorientation.angle./degree,1:110)
% option 2)
plotAngleDistribution(grains('q').boundary('q','q').misorientation)

% in each case, we not that there is a huge peak at 60degree,
% this is due to dauphine twinning, a 60 degree rotation around the c-axis

% -------------------------------------------------------------------------
% e) get rid of dauphine twins
% hint: you can use merge, however it comes with a new meanorientation of
% grains which might not be any of the twinned portions, hence setting a
% higer symmetry might be an option

% option 1) merge the grains
% get the grain boundaries
gb_qtz = grains.boundary('q','q');
 %define a rotation assosciated with dauphine twinning
rot = rotation('axis',Miller(0,0,0,1,ebsd('q').CS),'angle',60*degree);
% 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 threshold
ind = ang_r < 6.*degree;
%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.3)
hold on
plot(grains_merged('q').boundary,'linewidth',2)
hold on
plot(twinBoundary,'linecolor','r','linewidth',2)

% option 2) temporarily switch to higher 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


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

% -------------------------------------------------------------------------
% 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) plot a map of grains which have an aspect ratio larger 2
%create a new variable of qtz-only ebsd (makes some things easier)
gmqtz = grainsD('q');

%define an index for a condition
condA = gmqtz.aspectRatio > 2;
figure
plot(gmqtz(condA),gmqtz(condA).aspectRatio)
c= mtexColorbar; ylabel(c,'aspect ratio')


% -------------------------------------------------------------------------
% g) plot only grains with an aspect ratio larger 2 AND with a long axis
% direction +/- 10 degree from NW-SE, plot ipdf standar colorcoding
%get long axis direction
[ev,ew]= gmqtz.principalComponents; %ev being the angle as complex, ew the length of axs
laxdir = atan2(ev(1,2,:),ev(1,1,:))./degree.*-1;
laxdir=laxdir(:);
%dir could go ccw with 0 in the east
laxdir(laxdir<0) = 180 + laxdir(laxdir<0);


plot(gmqtz,laxdir)
mtexColorMap('hsv')
c = mtexColorbar; ylabel(c,'trend of long axis')

% get the angle between vector3d(2,1,0) and laxdir within range of 10
% degree, in our case 2,1,0 is eqal to 153.4 (otherwise one can use polar, take care for the correct direction)
[~, azi]= polar(vector3d(2,1,0));
condB = abs(azi./degree+90-laxdir) <= 10;
figure
plot(gmqtz(condA & condB),gmqtz(condA & condB).meanOrientation)


% -------------------------------------------------------------------------
% h) compare the area weighted pole figures of large and small grains and
% compare only the pole figures obtained by grains.meanOrientation
gmqtz = grainsD('q');
condB = gmqtz.equivalentRadius.*2 < 5;

figure
plotPDF(ebsd(gmqtz(condB),'q').orientations,h,'contourf','halfwidth',10*degree,'minmax','colorrange','equal')
figure
plotPDF(ebsd(gmqtz(~condB),'q').orientations,h,'contourf','halfwidth',10*degree,'minmax','colorrange','equal')

figure
plotPDF(gmqtz(condB).meanOrientation,h,'contourf','halfwidth',10*degree,'minmax','colorrange','equal')
figure
plotPDF(gmqtz(~condB).meanOrientation ,h,'contourf','halfwidth',10*degree,'minmax','colorrange','equal')


%% 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(grains_merged('q'),grains_merged('q').meanOrientation,'FaceAlpha',0.2)
hold on
plot(grains_merged('q').boundary('q','q'),grains_merged('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(grains_merged('q').CS);
plot(grains_merged('q').boundary('q','q'),oMC.Miller2color(grains_merged('q').boundary('q','q').misorientation.axis), 'linewidth',2)
figure
plot(oMC)

% c) plot a map with only those boundaries which have misorientation axes 
% close to 10-10

gb= grains_merged('q').boundary('q','q');
misax = gb.misorientation.axis;
cond = angle(misax,Miller(1,0,-1,0,misax.CS),'antipodal') < 15*degree;
plot(gb(cond),oMC.Miller2color(misax(cond)), 'linewidth',2)

% 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

gbo = ebsd_merged(grains_merged.boundary('q','q').ebsdId,'q').orientations;    
% get the misorientation axis
misaxS = axis(gbo(:,1),gbo(:,2));
oMS = TSLOrientationMapping(misaxS);
    plot(oMS,'zAxisDirection','intoPlane','xAxisDirection','east')
    hold on
    plot(misaxS,'antipodal','contour','contours',0:1:3)
    hold off
figure    
plot(grains_merged('q').boundary('q','q'),oMS.Miller2color(misaxS), 'linewidth',2)

% plot those misorientation axes in a pf/ipf with misorientation angles between 10 and
% 20degree 
misAng = grains_merged('q').boundary('q','q').misorientation.angle./degree;
cond = misAng >=10 & misAng <=20;
plot(misax(cond),'contourf','FundamentalRegion','antipodal','minmax')
figure
plot(misaxS(cond),'antipodal','contourf','halfwidth',10*degree,'minmax')



% -------------------------------------------------------------------------
%i) Big-Exercise: Compare misorientations of low angle boundaries (2-10 degree)
% of two types of grains, those with a flat and those with a steep 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.5*degree);

gmqtz = grainsD('q');
c_dir =gmqtz.meanOrientation.*Miller(0,0,0,1,ebsd('q').CS);
ebsdq =ebsd('q');

% 1) with the c-axis close to the center
cond1 = angle(c_dir, zvector,'antipodal') < 20*degree;
% 2) with the c-axis close to the periphery
cond2 = angle(c_dir, zvector,'antipodal') > 70*degree;

% check selection
figure
plot(gmqtz(cond1),gmqtz(cond1).meanOrientation)
figure
plot(gmqtz(cond2),gmqtz(cond2).meanOrientation)

%get the ebsd of selected  grains
ebsd1 = ebsdq(gmqtz(cond1));
ebsd2 = ebsdq(gmqtz(cond2));

% make new grains with a very small threshold
[subgrains_1, ebsd1.grainId,ebsd1.mis2mean] = calcGrains(ebsd1,'angle',0,5*degree,'unitcell');
[subgrains_2, ebsd2.grainId,ebsd2.mis2mean] = calcGrains(ebsd2,'angle',0,5*degree,'unitcell');


gb1 = subgrains_1.boundary('q','q');
gb2 = subgrains_2.boundary('q','q');


% misorientation angle
misor_angle1 = gb1.misorientation.angle./degree;
misor_angle2 = gb2.misorientation.angle./degree;

% misorientation axis (crytal ref. frame)
misor_ax_crystal1 = gb1.misorientation.axis;
misor_ax_crystal2 = gb2.misorientation.axis;


% misor axes (sample ref. frame)
%get the boundary orientations
gbo1 = ebsd1(gb1.ebsdId).orientations;    
gbo2 = ebsd2(gb2.ebsdId).orientations;
% calc the axes
misor_ax_sample1 = axis(gbo1(:,1),gbo1(:,2));
misor_ax_sample2 = axis(gbo2(:,1),gbo2(:,2));

% we can all the above to a specific range of angles
min_angle = 2;
max_angle = 10;

conda1 = misor_angle1 > min_angle & misor_angle1 < max_angle;
conda2 = misor_angle2 > min_angle & misor_angle2 < max_angle;


% set 1
% make a pleasant colorcoding for misor axes
oMC = TSLOrientationMapping(misor_ax_crystal1.CS,misor_ax_crystal1.CS);
figure
plot(ebsd1,ebsd1.orientations,'faceAlpha',0.3)
hold on
plot(gb1(conda1),oMC.Miller2color(misor_ax_crystal1(conda1)),'linewidth',2)
hold off

figure
plot(misor_ax_crystal1(conda1),misor_angle1(conda1),'MarkerSize',3,'antipodal','FundamentalRegion')
c = mtexColorbar(gcm);
mtexColorMap('WhiteJet')

ylabel(c,'misor angle')

figure
plot(misor_ax_crystal1(conda1),'contourf','antipodal','FundamentalRegion','minmax')
mtexColorbar;
mtexColorMap('LaboTeX')

figure
plot(misor_ax_sample1(conda1),'antipodal','minmax','MarkerSize',2)
hold on
plot(misor_ax_sample1(conda1),'contour','antipodal','minmax','linewidth',2,'halfwidth',10*degree)
hold off



% set 2
figure
plot(ebsd2('q'),ebsd2('q').orientations,'faceAlpha',0.3)
hold on
plot(gb2(conda2),oMC.Miller2color(misor_ax_crystal2(conda2)),'linewidth',2)
hold off


figure
plot(misor_ax_crystal2(conda2),misor_angle2(conda2),'MarkerSize',3,'antipodal','FundamentalRegion')
c = mtexColorbar(gcm);
mtexColorMap('WhiteJet')
ylabel(c,'misor angle')

figure
plot(misor_ax_crystal2(conda2),'contourf','antipodal','FundamentalRegion','minmax')
mtexColorbar;
mtexColorMap('LaboTeX')

figure
plot(misor_ax_sample2(conda2),'antipodal','minmax','MarkerSize',2)
hold on
plot(misor_ax_sample2(conda2),'contour','antipodal','minmax','linewidth',2,'halfwidth',10*degree)
hold off


% 