%% Import Script for EBSD Data
%
% This script was automatically created by the import wizard. You should
% run the whoole script or parts of it in order to import your data. There
% is no problem in making any changes to this script.

%% Specify Crystal and Specimen Symmetries
% Corrected your code here using the blue command Import Wizard in the Command Window
% or typing import_wizard('EBSD') in the Command Window
% 
% crystal symmetry
CS = {... 
  'notIndexed',...
  symmetry('2/m', [9.794 8.906 5.319], [90,105.9,90]*degree, 'X||a*', 'Y||b*', 'Z||c', 'mineral', 'Clinopyroxene', 'color', 'light blue'),...
  symmetry('mmm', [4.756 10.207 5.98], 'mineral', 'Forsterite', 'color', 'light green'),...
  symmetry('2/m', [43.3 9.23 7.27], [90,91.6,90]*degree, 'X||a*', 'Y||b*', 'Z||c', 'mineral', 'Antigorite superstruc', 'color', 'light red')};

% specimen symmetry
SS = symmetry('triclinic');

% plotting convention your may need to change this depending on your SEM
% configuration
setMTEXpref('xAxisDirection','north');
setMTEXpref('zAxisDirection','outOfPlane');

%% Specify File Names

% path to files
pname = '/MatLab_Programs';

% which files to be imported
fname = [pname '/area1.ctf'];

%% Import the Data

% create an EBSD variable containing the data
ebsd = loadEBSD(fname,CS,SS,'interface','ctf',...
  'convertEuler2SpatialReferenceFrame')

%%
% phases use for subsequent calculations: 
% phase 8:olivine;phase 9:
% antigorite
phases=[8,9];

% Alternatively enter user defined volume fractions (olivine: 0.8; antigorite:0.2)
volFrac=[0.8 0.2];

% define tensor frame
cs_Tensor_Olivine = symmetry('mmm',[4.7646,10.2296,5.9942],...
    [90.00,90.00,90.00]*degree,'x||a','z||c','mineral','Olivine')

% Olivine Density : 3.355 gm/cm3
% TYPING ERROR
% *** rho_olivine to rho_Olivine ***
rho_Olivine=3.355;
% The input Olivine tensor as matrix M_Epidote
M_Olivine = ....
    [[320.5  68.15  71.6     0     0     0];...
    [ 68.15  196.5  76.8     0     0     0];...
    [  71.6   76.8 233.5     0     0     0];...
    [   0      0      0     64     0     0];...
    [   0      0      0      0    77     0];...
    [   0      0      0      0     0  78.7]];
% Define tensor Olivine stiffness tensor with MTEX command tensor
C_Olivine=tensor(M_Olivine,cs_Tensor_Olivine)

% define tensor frame
cs_Tensor_Antigorite= symmetry('monoclinic',[43.5852,9.2624,7.246],...
[90.0,91.160,90.0]*degree,'X||a*','Z||c','mineral','Antigorite')

% Antigorite Density : 2.62gm/cm3
rho_Antigorite=2.62;
% The input Antigorite tensor as matrix M_Epidote
M_Antigorite =....
    [[208.1   66.4   16.0   0.0   5.5   0.0];...
    [  66.4  201.6    4.9   0.0  -3.1   0.0];...
    [  16.0    4.9   96.9   0.0   1.6   0.0];...
    [   0.0    0.0    0.0  16.9   0.0 -12.1];...
    [   5.5   -3.1    1.6   0.0  18.4   0.0];...
    [   0.0    0.0    0.0 -12.1   0.0  65.5]];
% Define tensor Antigorite stiffness tensor with MTEX command tensor
C_Antigorite=tensor(M_Antigorite,cs_Tensor_Antigorite)
%%

rho_aggregate = 0.8*rho_Olivine + 0.2*rho_Antigorite;

%%
% Olivine and Antigorite aggregate: Voigt,Reuss, Hill averages
%
% This only works if ALL phases in EBSD object are used 
% The example in the on-line help is a two-phase aggregate 
% The Voigt, Reuss, and Hill averages for all phases are computed by
% [CVoigt,CReuss,CHill] =  calcTensor(ebsd,CGlaucophane,CEpidote)
%
%[C_Voigt_Ag_EBSD,C_Reuss_Ag_EBSD,C_Hill_Ag_EBSD]=...
%    calcTensor(ebsd,C_Olivine,C_Antigorite)

%%
% Alternative method - this works
%
[C_Voigt_Ol_EBSD,C_Reuss_Ol_EBSD,C_Hill_OL_EBSD]=...
    calcTensor(ebsd('Forsterite'),C_Olivine)

[C_Voigt_Atg_EBSD,C_Reuss_Atg_EBSD,C_Hill_Atg_EBSD]=...
    calcTensor(ebsd('Antigorite'),C_Antigorite)

C_Hill_Ag_EBSD = 0.8*C_Hill_OL_EBSD + 0.2*C_Hill_Atg_EBSD

%%
% calculate Vp in yz-plane
%
% Define the 180 degress rotation about the x-specimen direction  = vector3d(1,0,0)
o = rotation('axis',vector3d(1,0,0),'angle',(0:180)*degree);
% Define a plane normal to x-specimen direction (yz plane)
% by the rotation of the starting vector in the yz plane, for example z or y
starting_vector = vector3d(0,0,1);
% generate an array of vectors in 1 degree intervals in yz-plane (rotation o about x)
yz_plane = o*starting_vector;
% generate velocities in yz-plane
[vp,vs1,vs2,pp,ps1,ps2] = velocity(C_Hill_Ag_EBSD,yz_plane,rho_aggregate);