% %
% Example of mapping  Olivine grain long and short axes
% and Olivine [100],[010] and [001] directions
% using mean grain orientation
%
% Using only MTEX quiver with arrows at grain centroid
%
% Using MTEX data set forsterite
%
%
% Magenta upper hemisphere
%   Green lower hemiphere
%
% David Mainprice 9/02/2017
%
%**************************************************************************
%% Clear memory and close all plots
clear all
close all
%% import MTEX example data set
mtexdata forsterite
%
% correct EBSD spatial coordinates
%
% rotate only the spatial data about the z-axis
ebsd = rotate(ebsd,180*degree,'keepEuler')
% rotate only the spatial data about the y-axis
ebsd = rotate(ebsd,rotation('axis',yvector,'angle',180*degree),'keepEuler')
% plot X to East
plotx2east
%**************************************************************************
%% Grain modelling prototcole - using indexed points and segmentation angle
% and re-calculate grain model to cleanup grain boundaries
%**************************************************************************
% segmentation angle typically 10 to 15 degrees
seg_angle = 10;
% minimum indexed points per grain between 5 and 10
min_points = 10;
% Restrict to indexed only points
[grains,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd('indexed'),'angle',seg_angle*degree);
% Remove small grains with less than min_points indexed points 
grains = grains(grains.grainSize > min_points);
% Re-calculate grain model to cleanup grain boundaries with less than minimum index points 
% used ebsd points within grains having the minium indexed number of points (e.g. 10 points)
ebsd = ebsd(grains);
[grains,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd('indexed'),'angle',seg_angle*degree);
plot(grains)
%**************************************************************************
%% Principle components using SVD analysis
%**************************************************************************
% Eigen-vectors (ev) Eigen-values (ew)
[ev,ew]= principalComponents(grains('fo'));
% Angle from X-axis using Eigen-vectors (ev)
Theta = atan2(ev(1,2,:),ev(1,1,:))/degree;
% Remove unwanted singleton dimensions
Theta=squeeze(Theta);
% keep angles in the range 0-180 degrees from X-axis
% Theta(Theta<0) = 180 + Theta;
for i=1:length(Theta)
    if(Theta(i)<0)
       Theta(i) = 180 + Theta(i);
    end
end
% theta for all grains zero azimuth from specimen X-axis 
% vector3d for all long and short axes
%**************************************************************************
% Create long and short axis vectors in specimen co-ordinates
%**************************************************************************
polar_angle_array(1:grains('fo').length) = 90.0;
% Remove unwanted singleton dimensions
polar_angle_array = polar_angle_array(1,:)';
azimuth_angle_array = Theta;
% create long-axis specimen 3dvector arrays 
long_axis_vectors = vector3d('polar',polar_angle_array*degree,azimuth_angle_array*degree);
% create short-axis specimen 3dvector arrays 
short_azimuth_angle_array = Theta + 90.0;
short_axis_vectors  = vector3d('polar',polar_angle_array*degree,short_azimuth_angle_array*degree);
%**************************************************************************
%% plot long(magenta) and short(green) axes in map coordinates at grain centroid
%**************************************************************************
Scale = 0.25;
close all
figure
% plot grains
plot(grains('fo'))
hold on
plot(grains('fo').boundary,'lineColor','k','LineWidth',2)
hold on
% plot grain short-axes in specimen coordinates
quiver(grains('fo'),Scale*grains('fo').diameter.*short_axis_vectors,'Color','g','LineWidth',2,'autoScale','off')
hold on
quiver(grains('fo'),Scale*grains('fo').diameter.*-short_axis_vectors,'Color','g','LineWidth',2,'autoScale','off','displayName','short-axis')
hold on
% plot grain long-axes in specimen coordinates
quiver(grains('fo'),Scale*grains('fo').diameter.*long_axis_vectors,'Color','m','LineWidth',2,'autoScale','off')
hold on
quiver(grains('fo'),Scale*grains('fo').diameter.*-long_axis_vectors,'Color','m','LineWidth',2,'autoScale','off','displayName','long-axis')
hold off
%**************************************************************************
%% Create crystal directions in specimen coordinates
%**************************************************************************
% use grain mean orientaton
r_100_mO = grains('fo').meanOrientation*Miller(1,0,0,ebsd('fo').CS,'uvw');
r_010_mO = grains('fo').meanOrientation*Miller(0,1,0,ebsd('fo').CS,'uvw');
r_001_mO = grains('fo').meanOrientation*Miller(0,0,1,ebsd('fo').CS,'uvw');
% Convert all to upper hemisphere
r_100_mO(r_100_mO.z<0)=-r_100_mO(r_100_mO.z<0);
r_010_mO(r_010_mO.z<0)=-r_010_mO(r_010_mO.z<0);
r_001_mO(r_001_mO.z<0)=-r_010_mO(r_001_mO.z<0);
%**************************************************************************
%% plot crystallograpic axes in map coordinates at grain centroid
% a-axis red b-axis green c-axis blue upper hemisphere directions
%**************************************************************************
% autoScale Factor default = 0.25
Scale_a = 0.05;
Scale_b = 0.02;
Scale_c = 0.02;
close all
figure
% plot grains
plot(grains('fo'))
hold on
% plot grain boundary outline
plot(grains('fo').boundary,'lineColor','k','LineWidth',1)
hold on
% plot crystallograpic a-axes in map coordinates at grain centroid
quiver(grains('fo'),Scale_a*grains('fo').diameter.*r_100_mO,'Color','r','LineWidth',2,'autoScale','off','displayName','a-axis')
hold on
% plot crystallograpic b-axes in map coordinates at grain centroid
quiver(grains('fo'),Scale_b*grains('fo').diameter.*r_010_mO,'Color','g','LineWidth',2,'autoScale','off','displayName','b-axes')
hold on
% plot crystallograpic c-axes in map coordinates at grain centroid
quiver(grains('fo'),Scale_c*grains('fo').diameter.*r_001_mO,'Color','b','LineWidth',2,'autoScale','off','displayName','c-axes')
hold off
drawNow(gcm,'figSize','large')


