%% MTEX workshop 2021
% Exercise of physical properties
% Seismic properties of mantle rock

%% 1. Import orientation map from sample SAM-1-X.ang

% crystal symmetry
CS = {... 
  'notIndexed',...
  crystalSymmetry('121', [9.6 8.8 5.3], [90,106.85,90]*degree, 'X||a*', 'Y||b*', 'Z||c', 'mineral', 'diopside', 'color', [0.53 0.81 0.98]),...
  crystalSymmetry('222', [18 8.8 5.2], 'mineral', 'enstatite', 'color', [0.56 0.74 0.56]),...
  crystalSymmetry('222', [4.8 10 6], 'mineral', 'olivine', 'color', [0.85 0.65 0.13]),...
  crystalSymmetry('432', [8.1 8.1 8.1], 'mineral', 'spinel', 'color', [0.94 0.5 0.5])};

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

% path to files
pname = '/Users/lfgmorales/OneDrive - ETH Zurich/Dokumente/MTEX_workshop_2021_physical_properties/ang files';

% which files to be imported
fname = [pname '/SAM-1-X-clean.ang'];


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

%%

plot(ebsd('o'),ebsd('o').orientations)

%% 2. Section the map to remove the non-mapped [0 0 35624 15900]

region = [0 0 35624 15900]
close all
figure
oM2 = ipfHSVKey(ebsd('olivine'))
color = oM2.orientation2color(ebsd('olivine').orientations);
plot(ebsd('olivine'),color,'figSize','medium')
rectangle('position',region,'edgecolor','r','linewidth',2)
condition = inpolygon(ebsd,region);
ebsd = ebsd(condition);

%% 3. Rotate the data to bring foliation to the horizontal
%
% first rotation should be -45 degrees in zvector, second should be 90 
% degrees in yvector

ebsd_rot1= rotate(ebsd,rotation('axis',zvector,'angle',-45*degree));
ebsd_rot2= rotate(ebsd_rot1,rotation('axis',yvector,'angle',90*degree));
ebsd=ebsd_rot2

plot(ebsd('olivine'),color,'figSize','medium')

%% 4. Import the olivine, enstatite and diopside elastic constants

%% diopside elastic tensor

cs_Tensor_cpx = crystalSymmetry('121',[9.585  8.776  5.26],...
  [90.0000 105.8600  90.0000]*degree,'x||a*','z||c',...
  'mineral','Diopside');

% the density
rho_cpx = 3.2860;

% the tensor coefficients
Cij_cpx =....
  [[  228.10   78.80   70.20    0.00    7.90    0.00];...
  [    78.80  181.10   61.10    0.00    5.90    0.00];...
  [    70.20   61.10  245.40    0.00   39.70    0.00];...
  [     0.00    0.00    0.00   78.90    0.00    6.40];...
  [     7.90    5.90   39.70    0.00   68.20    0.00];...
  [     0.00    0.00    0.00    6.40    0.00   78.10]];

% define the tensor
C_cpx = stiffnessTensor(Cij_cpx,cs_Tensor_cpx,'density',rho_cpx);

%% olivine elastic tensor

cs_Tensor_olivine = crystalSymmetry('222', [4.762 10.225 5.994],...
    'mineral', 'olivine', 'color', 'light red');

% we need the density of this phase
rho_olivine = 3.3550;

% and we need the coefficients Cij for olivine in Voigt notation 
% Abramson et al. 1997 Journal of Geophysical Research)

Cij_ol = [[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]];

% the elastic stiffness tensor is defined as

C_olivine = stiffnessTensor(Cij_ol,cs_Tensor_olivine,'density',...
    rho_olivine);

%% enstatite elastic tensor

cs_Tensor_opx = crystalSymmetry('mmm',[ 18.2457  8.7984  5.1959],...
  [  90.0000  90.0000  90.0000]*degree,'x||a','z||c',...
  'mineral','Enstatite');

% the density
rho_opx = 3.3060;

% the tensor coefficients
Cij_opx =....
  [[  236.90   79.60   63.20    0.00    0.00    0.00];...
  [    79.60  180.50   56.80    0.00    0.00    0.00];...
  [    63.20   56.80  230.40    0.00    0.00    0.00];...
  [     0.00    0.00    0.00   84.30    0.00    0.00];...
  [     0.00    0.00    0.00    0.00   79.40    0.00];...
  [     0.00    0.00    0.00    0.00    0.00   80.10]];

% define the tensor
C_opx = stiffnessTensor(Cij_opx,cs_Tensor_opx,'density',rho_opx);

%% 5. Calculate the ODFs of the individual phases

odf_ol = calcDensity(ebsd('o').orientations,'halfwidth',10*degree);
odf_opx = calcDensity(ebsd('e').orientations,'halfwidth',10*degree);
odf_cpx = calcDensity(ebsd('d').orientations,'halfwidth',10*degree);

%% 6.  Define pole figures for the phases

PFs_olivine_uvw = [ ...
    Miller(1,0,0,ebsd('o').CS,'uvw'), ...
    Miller(0,1,0,ebsd('o').CS,'uvw'), ...
    Miller(0,0,1,ebsd('o').CS,'uvw'), ...
    ]

PFs_enstatite_uvw = [ ...
    Miller(1,0,0,ebsd('e').CS,'uvw'), ...
    Miller(0,1,0,ebsd('e').CS,'uvw'), ...
    Miller(0,0,1,ebsd('e').CS,'uvw'), ...
    ]

PFs_diopside_hkl = [ ...
    Miller(1,0,0,ebsd('d').CS,'hkl'), ...
    Miller(0,1,0,ebsd('d').CS,'hkl'), ...
    Miller(0,0,1,ebsd('d').CS,'hkl'), ...
    ]

%% 7. Plot pole figures

figure
plotPDF(odf_olivine,PFs_olivine_uvw,'upper','resolution',5*degree,'contourf','colorrange','equal','minmax')
colorbar
annotate([xvector,yvector,zvector],...
'label',{'X','Z','Y'},'backgroundcolor','w','MarkerSize',20,'FontSize',20)
mtexColorbar
drawNow(gcm,'figSize','large')

%% 8. Compute the average stiffness tensor for each phase individually

[CVoigt_ol, CReuss_ol, CHill_ol]    = mean(C_olivine,odf_ol);
[CVoigt_opx, CReuss_opx, CHill_opx] = mean(C_opx,odf_opx);
[CVoigt_cpx, CReuss_cpx, CHill_cpx] = mean(C_cpx,odf_cpx);


%% 9. Calculate the volume of each phase in this map

vol_ol  = length(ebsd('o')) ./ length(ebsd('indexed'));
vol_opx = length(ebsd('e')) ./ length(ebsd('indexed'));
vol_cpx = length(ebsd('d')) ./ length(ebsd('indexed'));

%% 10. Compute the polycrystalline stiffness tensor for the EBSD map

CHill = vol_ol * CHill_ol + vol_opx * CHill_opx...
    + vol_cpx * CHill_cpx;

%% 11. Plot the seismic velocities

close all
plotSeismicVelocities(CHill)
saveFigure ('seismic velocities aggregate.png')



