%% Import Script for EBSD Data
%**************************************************************************
% Version only plots and prints the most commonly used Hill average
% also known as Voigt-Reuss-Hill average
%**************************************************************************
% Version using the harmonic method with ODF calculated using the Dirichlet
% kernel and a Bandwidth of 4 for 4th rank elastic tensor properties
%**************************************************************************
%
% Tested with MTEX 4.0.9
%
% Reference for EBSD dataset used in this example
%
% Andrea Tommasi, Alain Vauchez, Dmitri A. Ionov (2008)
% Deformation, static recrystallization, and reactive melt transport in shallow
% subcontinental mantle xenoliths (Tok Cenozoic volcanic field, SE Siberia)
% Earth and Planetary Science Letters 272 (2008) 65-77.
%
%
% Worked MTEX/MatLab examples for the seismic properties from EBSD
% 
% A) Import of Oxford/HKL *.ctf (ASCII) file
%
% B) Area fractions of all indexed phases
%
% C) Correction of area fractions
%
% D) Use pre-defined set of single crystal elastic tensors
%    of Olivine, Enstatite and Diopside
%    Analysis their anisotropy and plot velocity diagrams
%
% E) Calculate single phase aggregate elastic tensors of
%    Olivine, Enstatite and Diopside using Voigt,Hill and Ruess averages
%    Analysis their anisotropy and plot velocity diagrams
%
% F) Calculate multi-phase aggregate elastic tensors of
%    Olivine + Enstatite + Diopside using Voigt,Hill and Ruess averages
%    Analysis their anisotropy and plot velocity diagrams
% 
%****** functions *********************************************************
% Func_Print_Seismic_Cij_MTEX4.m
% Func_Seismic_multiplot_gfs_MTEX4.m
% ALL MUST BE IN YOUR MATLAB CURRENT DIRECTORY
% (CHECK top of Command Window for current directory)
%**************************************************************************
%
% Uses *.ctf data file at pname = '/MatLab_Programs/MTEX_Workshop_Distibution';
%**************************************************************************
% set font size
% N.B. You need to adjust this to your screen size,
% for example try 20, 30 or 40
% set here to global_fontSize = 20
% for Func_Seismic_multiplot_gfs_MTEX4
% **********************************************************************
%
% Revised 19/01/2015 by David Mainprice
%
% Reference publication
%
% Mainprice, D., Hielscher, R., Schaeben, H. (2011) 
% Calculating anisotropic physical properties from texture data using
% the MTEX open source package. In: Prior, D.J., Rutter, E.H., 
% Tatham, D. J. (eds) Deformation Mechanisms, Rheology and Tectonics:
% Microstructures, Mechanics and Anisotropy. 
% Geological Society, London, Special Publications, 360, 175-192.
%
%% clear memory and close plots - for new analysis of EBSD data
clear all
close all
%**************************************************************************
%% set global_fontSize (varies as function of your screen resolution)
% for Func_Seismic_multiplot_gfs_MTEX4 
global_fontSize = 20;
%% Specify Crystal and Specimen Symmetries
%
% crystal symmetry
CS = {... 
  'notIndexed',...
  crystalSymmetry('mmm', [4.756 10.207 5.98], 'mineral', 'Forsterite', 'color', 'dark green'),...
  crystalSymmetry('mmm', [18.2406 8.8302 5.1852], 'mineral', 'Enstatite', 'color', 'light red'),...
  crystalSymmetry('12/m1', [9.746 8.99 5.251], [90,105.63,90]*degree, 'X||a*', 'Y||b*', 'Z||c', 'mineral', 'Diopside', 'color', 'dark blue'),...
  crystalSymmetry('m-3m', [5.431 5.431 5.431], 'mineral', 'Silicon', 'color', 'cyan')};

% specimen symmetry : the MTEX4 default is triclinic so do not need this line
SS = specimenSymmetry('triclinic'); 
% plotting convention
setMTEXpref('xAxisDirection','east');
setMTEXpref('zAxisDirection','outOfPlane');

%% Specify File Names
% path to files
pname = '/MatLab_Programs/MTEX_Workshop_Distibution';
% which files to be imported
fname = [pname '/9510-16-map-stitch-pscleaned.ctf'];
%% Import the Data
% create an EBSD variable containing the data
ebsd = loadEBSD(fname,CS,'interface','ctf',...
  'convertSpatial2EulerReferenceFrame', 'wizard')

%% Generate simple Crystal Symmetry variables for all phases (tested with MTEX4)
fprintf(' \n');
fprintf(' Crystal Symmetry (CS) : names of variables \n');
fprintf(' \n');
% phase names
Phase_names = ebsd.mineralList;
for i=1:length(CS)
% indexed point for mineral
     N_Points = length(ebsd(Phase_names(i)));
% print only indexed phases
    if((~strcmpi(Phase_names{i},'notIndexed')) && (N_Points > 0))
% if minus sign in name remplace it with underscore
      Phase_names{i} = strrep(Phase_names{i},'-','_');
% retain first part of mineral name
      P_Name = strtok(char(Phase_names{i}),' ');
% make variables for CS with mineral names
      myvariable = strcat(P_Name,'_CS');
      datavalues = CS{i};
% print variable names for CS
      eval([sprintf(myvariable) ' = CS{i};'])
      fprintf(' %s %s \n','Variable name =',myvariable);
    end
end
disp(' ')
disp(' Use mineral_CS to defined your CS in m-file ')
disp(' ')
%% Default MTEX colormap - jet colormap begin with white
setMTEXpref('defaultColorMap',WhiteJetColorMap);
%**************************************************************************
%% Set-up an export file with area fractions results
%**************************************************************************
% it will be in folder with the datafile called *._Area_fractions.txt
filename_export = '/MatLab_Programs/Forsterite_MTEX_Example_Area_fractions.txt';
% data file name
data_file = 'Forsterite_MTEX_Example_Area_fractions.txt';
% open export file to write 'w'
[fid2] = fopen(filename_export, 'w');
fprintf(fid2,'%s \n','Area fractions analysis');
fprintf(fid2,' \n');
fprintf(fid2,'%s %s \n','Data file:',data_file);
fprintf(fid2,' \n');
%***********************************************************************
%% STEP 1 : You must create an ODF for contour plots
%***********************************************************************
%  to calcluated Texture Index, to plot contoured ODF, Pole Figures,
% Inverse Pole figures ...
% Calculate an ODF for the mineral Olivine
% Here we define the half-width of radial bell-shaped kernel function
% The de vallee Poussin is the default kernel function for texture analyisis 
odf_olivine = calcODF(ebsd('Forsterite').orientations,'halfwidth',10*degree)
%***********************************************************************
%% STEP 2 : create a list of hkls or uvws for your pole figure plots
%***********************************************************************
% list two : [uvw] [directions]
pfs_olivine_uvw = (...
    [Miller(1,0,0,Forsterite_CS,'uvw'),...
     Miller(0,1,0,Forsterite_CS,'uvw'),...
     Miller(0,0,1,Forsterite_CS,'uvw')])
%***********************************************************************
%% STEP 3 : contoured pole figures
%***********************************************************************
close all
% Options:
%
% 'antipodal', also known as non-polar, the +hkl and -hkl directions
% are rotated into the same hemisphere (default upper),(e.g. 001 and 00-1),
% resulting all data is plotted in one hemisphere. This typically the
% default setting used by most people.
%
% Alternative options
% 'complete' = whole hemisphere, plots 'upper' and 'lower' hemispheres.
% 'upper' = plots data in the upper hemisphere.
% 'lower' = plots data in the lower hemisphere.
%
% Contouring options:
% 'resolution' - resolution in degrees
% 'contour' - colour contour lines
% 'contourf' - filled color between contours
% 'logarithmic' - log scale
% 'gray' - gray colour scale
%
% figure position vectors - [left bottom width height]
%       I------------I   I 
%       I            I   I
%       I<--width--->I height
%       I            I   I
%       I------------I   I
%--left--      
%       I
%    bottom
%       I
% N.B. plot to EAST
plotx2east
figure
plotPDF(odf_olivine,pfs_olivine_uvw,'lower','contourf')
colorbar
% annotate with text labels X,Y,Z at positions in the pole figure
% specimen coordinates are given by MTEX vector3d X=xvector ,Y=yvector, Z=zvector
% structural geology coordinates X=lineation Y = normal to X and Y is  Z=pole
% to foliation, so that MTEX xvector = Zs MTEX yvector = X
% and  MTEX -zvector (lower hemisphere) =  Y
% In this example the structural directions are approximate
annotate([xvector,yvector,-zvector],...
'label',{'Z','Y','X'},'backgroundcolor','w','MarkerSize',20,'FontSize',20)
% use colorbar if you all pole figure contoured on the SAME SCALE
% remove colorbar if want the pole figures contoured INDIVIDUALLY
%***********************************************************************
% save graphics file
%***********************************************************************
% file name  plus .pdf, .eps. , .png
saveFigure('/MatLab_Programs/PFs_LH_Forsterite_MTEX_Example.pdf')
%%
%***********************************************************************
% Elminate poor EBSD measurements with low fit parameter called MAD
%***********************************************************************
%
% Get the Mean Angular Deviation (MAD) for this EBSD data set using
% an Oxford-HKL system
%
% Low MAD means good fit for the indexing and diffraction pattern
%
close all
% histogram of mad
hist(ebsd.mad)
title('Mean Angular Deviation (MAD)','FontSize',20)
% define maximum acceptable MAD value, typically 1.3 degrees
MAD_MAXIMUM = 1.3;
% Eliminate all measurements with MAD larger than MAD_MAXIMUM=1.3
% create new ebsd object with corrected data set
ebsd_corrected = ebsd(ebsd.mad<MAD_MAXIMUM)
%%
close all
%%
%***********************************************************************
% Set up volume fractions                            (tested with MTEX4)
%***********************************************************************
% list of phase names
Phase_names = ebsd.mineralList;
% zero area fractions
Area_Fractions = zeros(numel(Phase_names));
% number of measurement points in the map
Map_points = length(ebsd_corrected);
% default first phase is index
N_first_indexed_phase = 1;
% remove nonIndexed points from total
% 'notIndexed' is present
if(strcmp(Phase_names(1),'notIndexed'))
    N_first_indexed_phase = 2;
end
%**************************************************************************
% Area fractions for all phases in the EBSD object including non-indexed
%**************************************************************************
% print results to screen
fprintf(' \n')
fprintf(' Area fractions for all phases in the EBSD object \n')
fprintf(' \n')
fprintf('   #      Phase       Points  Area fraction \n')
for i=1:numel(Phase_names);
% indexed point for mineral
N_Points = length(ebsd_corrected(Phase_names(i)));
  if(N_Points > 0)
% area fraction for mineral
    Area_Fractions(i) = N_Points/Map_points;
% retain first part of mineral name
    P_Name = strtok(char(Phase_names{i}),' ');
    fprintf(' %3i %12s %10i %11.4f \n',...
    i,P_Name,N_Points,Area_Fractions(i))
  end
end
AF_total = sum(sum(Area_Fractions));
fprintf(' \n')
fprintf(' Total area fraction = %8.4f \n',AF_total)
fprintf(' \n')
% print results to file
fprintf(fid2,' \n');
fprintf(fid2,' Area fractions for all phases in the EBSD object \n');
fprintf(fid2,' \n');
fprintf(fid2,'   #      Phase       Points  Area fraction \n');
for i=1:numel(Phase_names);
% indexed point for mineral
N_Points = length(ebsd_corrected(Phase_names(i)));
  if(N_Points > 0)
% retain first part of mineral name
    P_Name = strtok(char(Phase_names{i}),' ');
    fprintf(fid2,' %3i %12s %10i %11.4f \n',...
    i,P_Name,N_Points,Area_Fractions(i));
  end
end
fprintf(fid2,' \n');
fprintf(fid2,' Total area fraction = %8.4f \n',AF_total);
fprintf(fid2,' \n');
% total number of indexed points in map without 'notIndexed' phase
Indexed_map_points = sum(length(ebsd_corrected(Phase_names(N_first_indexed_phase:numel(Phase_names)))));
% percentage of indexed point in map
Percent_indexed = 100*Indexed_map_points/Map_points;
% print to screen
fprintf(' \n')
fprintf(' Percentage indexed points in map = %5.1f \n',Percent_indexed)
fprintf(' \n')
% print to file
fprintf(fid2,' \n');
fprintf(fid2,' Percentage indexed points in map = %5.1f \n',Percent_indexed);
fprintf(fid2,' \n');
%**************************************************************************
% Area fractions based on minerals with EBSD indexed points only
%**************************************************************************
% print results to screen
fprintf(' \n')
fprintf(' Area fractions based on minerals with EBSD indexed points only \n')
fprintf(' \n')
fprintf('   #      Phase       Points  Area fraction \n')
for i=N_first_indexed_phase:numel(Phase_names);
% indexed point for mineral
N_Points = length(ebsd_corrected(Phase_names(i)));
  if(N_Points > 0)
% area fraction for mineral
    Area_Fractions(i) = N_Points/Indexed_map_points;
% retain first part of mineral name
    P_Name = strtok(char(Phase_names{i}),' ');
    fprintf(' %3i %12s %10i %11.4f \n',...
    i,P_Name,N_Points,Area_Fractions(i))
  end
end
AF_total = sum((Area_Fractions(N_first_indexed_phase:numel(Phase_names))));
fprintf(' \n')
fprintf(' Total area fraction = %8.4f \n',AF_total)
fprintf(' \n')
% print results to file
fprintf(fid2,' \n');
fprintf(fid2,' Area fractions based on minerals with EBSD indexed points only \n');
fprintf(fid2,' \n');
fprintf(fid2,'   #      Phase       Points  Area fraction \n');
for i=N_first_indexed_phase:numel(Phase_names);
% indexed point for mineral
N_Points = length(ebsd_corrected(Phase_names(i)));
  if(N_Points > 0)
% area fraction for mineral
    Area_Fractions(i) = N_Points/Indexed_map_points;
% retain first part of mineral name
    P_Name = strtok(char(Phase_names{i}),' ');
    fprintf(fid2,' %3i %12s %10i %11.4f \n',...
    i,P_Name,N_Points,Area_Fractions(i));
  end
end
AF_total = sum((Area_Fractions(N_first_indexed_phase:numel(Phase_names))));
fprintf(fid2,' \n');
fprintf(fid2,' Total area fraction = %8.4f \n',AF_total);
fprintf(fid2,' \n');
%% 
%***********************************************************************
% Revised Area fraction calculation
%***********************************************************************
% set area fraction for 'notIndexed' to zero
if(strcmp(Phase_names(1),'notIndexed'))
    Area_Fractions(1)=0.0;
end
fprintf(' \n')
fprintf(' ********************************* \n')
fprintf(' Revised Area fraction calculation \n')
fprintf(' ********************************* \n')
fprintf(' \n')
fprintf(' You may use the volume fractions determined by the \n')
fprintf(' script from your *.ctf EBSD file. \n')
fprintf(' Alternatively you can enter area fractions \n')
fprintf(' determined by another method or remove some phases \n')
fprintf(' (e.g. optical microscope, electron probe etc.) \n')
fprintf(' MAKE SURE YOU DO NOT HAVE THE Si CALIBRATION STANDARD \n')
fprintf(' IN YOUR LIST OF INDEXED PHASES - Silicon Points > 0 \n')
fprintf(' IF SO REMOVE IT BY ANSWERING NO below and \n')
fprintf(' set Si area fraction to 0.0 \n')
fprintf(' \n')
fprintf(' \n')
fprintf(' Use the volume fractions given \n')
ANS = input('by the script from your EBSD file(Y/N) ?','s');
if(strcmp(ANS,'N')||strcmp(ANS,'n'))
   fprintf(' \n')
   fprintf(' You have chosen to enter the volume fractions \n')
   fprintf(' determined by the another method or remove some phases \n')
   fprintf(' \n')
   fprintf(' The total area fraction will be automatically re-adjusted to 1.0 \n')
   fprintf(' (phases with EBSD zero indexed points are not \n')
   fprintf(' NORMALLY considered, for example Si standard !) \n')
   AF_total=0.0;
   for i=N_first_indexed_phase:numel(Phase_names);
       if(Area_Fractions(i) > 0)
          fprintf(' \n')
          P_Name = strtok(char(Phase_names{i}),' ');
          fprintf(' Phase No. %2i %s  \n',i,P_Name)
          fprintf('  OLD area fraction = %8.4f \n',Area_Fractions(i))
          fprintf(' -1 to keep original value \n')
          fprintf('  0 to remove phase \n')
          AF_NEW = input('NEW area fraction ?');
% record new area fraction
          if(AF_NEW >= 0)
            Area_Fractions(i) = AF_NEW;
          end
          fprintf('  *Revised* area fraction = %8.4f \n',Area_Fractions(i))
          AF_total = AF_total + Area_Fractions(i);
        end
   end
end
% re-adjust total area fraction to 1.00
Area_Fractions = Area_Fractions ./ AF_total;
% print results to screen
fprintf(' \n');
if(strcmp(ANS,'N')||strcmp(ANS,'n'))
  fprintf(' *Revised* Area fractions based on minerals with EBSD indexed points only \n');
end
if(strcmp(ANS,'Y')||strcmp(ANS,'y'))
  fprintf(' Area fractions based on minerals with EBSD indexed points only \n');
end
fprintf(' \n');
fprintf('   #      Phase       Points  Area fraction \n');
AF_total=0.0;
for i=N_first_indexed_phase:numel(Phase_names);
% indexed point for mineral
N_Points = length(ebsd_corrected(Phase_names(i)));
    if(Area_Fractions(i) > 0)
% retain first part of mineral name
      P_Name = strtok(char(Phase_names{i}),' ');
      fprintf(' %3i %12s %10i %11.4f \n',...
      i,P_Name,N_Points,Area_Fractions(i));
      AF_total = AF_total + Area_Fractions(i);
    end
end
AF_total = sum(sum(Area_Fractions));
fprintf(' \n');
fprintf(' Total area fraction = %8.4f \n',AF_total);
fprintf(' \n');
%
% print results to file
fprintf(fid2,' \n');
if(strcmp(ANS,'N')||strcmp(ANS,'n'))
  fprintf(fid2,' *Revised* Area fractions based on minerals with EBSD indexed points only \n');
end
if(strcmp(ANS,'Y')||strcmp(ANS,'y'))
  fprintf(fid2,' Area fractions based on minerals with EBSD indexed points only \n');
end
fprintf(fid2,' \n');
fprintf(fid2,'   #      Phase       Points  Area fraction \n');
for i=N_first_indexed_phase:numel(Phase_names);
% indexed point for mineral
N_Points = length(ebsd_corrected(Phase_names(i)));
    if(Area_Fractions(i) > 0)
% retain first part of mineral name
      P_Name = strtok(char(Phase_names{i}),' ');
      fprintf(fid2,' %3i %12s %10i %11.4f \n',...
      i,P_Name,N_Points,Area_Fractions(i));
    end
end
AF_total = sum(sum(Area_Fractions));
fprintf(fid2,' \n');
fprintf(fid2,' Total area fraction = %8.4f \n',AF_total);
fprintf(fid2,' \n');

%%
% Assign area fractions to variables 'af_Mineral_Name'
%
fprintf(' \n');
fprintf(' Area fractions : names of variables \n');
fprintf(' \n');
for i=N_first_indexed_phase:numel(Phase_names);
% indexed point for mineral
N_Points = length(ebsd_corrected(Phase_names(i)));
    if(Area_Fractions(i) > 0)
% if minus sign in name remplace it with underscore
      Phase_names{i} = strrep(Phase_names{i},'-','_');
% retain first part of mineral name
      P_Name = strtok(char(Phase_names{i}),' ');
% make variables for area fractions with mineral names
      myvariable = strcat('af_',P_Name);
      datavalues = Area_Fractions(i);
% print valiable names and area fraction 
      eval([sprintf(myvariable) '=Area_Fractions(i);'])
      fprintf(' %s %s \n','Variable name =',myvariable);
      fprintf(' %s %8.4f \n',' Area Fraction  =',datavalues);
    end
end
disp(' ')
disp(' The af_mineral_name above will be used for area fractions in this m-file ')
disp(' ')
%***********************************************************************
% close export text file
%***********************************************************************
fclose(fid2);
%
%%
%***********************************************************************
% Single crystal tensors
%***********************************************************************
%
%              NORTH
%                X
%                |       
%         --------------------
%         |Single Crystal Vp |
%    Y ---|PHYSICAL       Vs1|
%         |PROPERTIES     Vs2|
%         --------------------
%             
%**************************************************************************
% plotting convention for sample direction X to North for
% UPPER HEMISPHERE single crystal plots
% set colour map to seismic color map : blue2redColorMap
setMTEXpref('defaultColorMap',blue2redColorMap);
%%
%***********************************************************************
% Define elastic stiffness tensor of Olivine
%***********************************************************************
% Olivine elastic stiffness (Cij) Tensor in GPa
% Abramson E.H., Brown J.M., Slutsky L.J., and Zaug J.(1997)
% The elastic constants of San Carlos olivine to 17 GPa.
% Journal of Geophysical Research 102: 12253-12263.
%
% Enter tensor as 6 by 6 matrix,M line by line.
M = [[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 density : San Carlos OLIVINE Fo90 room pressure d=3.355 g/cm3
rho_olivine = 3.3550;
% Define Cartesian Tensor crystal symmetry & frame
CS_Tensor_olivine = crystalSymmetry('mmm',[4.7646 10.2296 5.9942],...
 [90.0000  90.0000  90.0000]*degree,'x||a','z||c', 'mineral','Olivine')
% Define elastic stiffness tensor Cijkl with MTEX load tensor command
C_olivine = tensor(M,CS_Tensor_olivine,'name','elastic stiffness','unit','GPa','density',rho_olivine)
%%
% Plot and Print out results
%
% print
elastic_tensor = 'Single crystal Olivine';
print_file = '/MatLab_Programs/Single_Crystal_Olivine_Seismic_Analysis.txt';
Func_Print_Seismic_Cij_MTEX4(C_olivine,rho_olivine,elastic_tensor,print_file)
% plot
close all
plot_x = 1; % plot X to North
Vp_polar = 1; % plot Vp polarizations
HS = 'upper'; % plot upper hemisphere
Vp_polar = 1; % plot Vp polarizations
% define Vector3d for reference directions as for single crystal
ref_vec_1 = xvector;
ref_vec_2 = yvector;
ref_vec_3 = zvector;
% reference axes - text labels
ref_axis_1 ='a';
ref_axis_2 ='b';
ref_axis_3 ='c';
plot_name = '/MatLab_Programs/Single_Crystal_Olivine_Seismic_plot.png';
Func_Seismic_multiplot_gfs_MTEX4(C_olivine,rho_olivine,plot_name,plot_x,...
Vp_polar,ref_axis_1,ref_axis_2,ref_axis_3,ref_vec_1,ref_vec_2,ref_vec_3,HS,global_fontSize)
%%
%***********************************************************************
% Define elastic stiffness tensor of Enstatite
%***********************************************************************
%
% Reference Elastic constants
% M.Chai, J.M.Brown and L.J.Slutsky and J.Zang JGR 1997 v.102 pp14779-14785.      
%
% Reference Crystal Structure
% Enstatite (En 0.75 Fs 0.08 Ts 0.12) density 3.306 g/cm3 room pressure           
%
% Define density (g/cm3)
rho_opx = 3.3060;
%
% Define Cartesian Tensor crystal symmetry & frame
cs_Tensor_opx = crystalSymmetry('mmm',[ 18.2457  8.7984  5.1959],...
 [  90.0000  90.0000  90.0000]*degree,'x||a','z||c',...
 'mineral','Enstatite');
%
% elastic Cij stiffness tensor (GPa) as matrix M
%
 M =.... 
 [[  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]];
%
% M as stiffness tensor C with MTEX tensor command
C_opx = tensor(M,cs_Tensor_opx);
%
% Plot and Print out results
%
% print
elastic_tensor = 'Single crystal Enstatite';
print_file = '/MatLab_Programs/Single_Crystal_Enstatite_Seismic_Analysis.txt';
Func_Print_Seismic_Cij_MTEX4(C_opx,rho_opx,elastic_tensor,print_file)
% plot
close all
plot_x = 1; % plot X to North
Vp_polar = 1; % plot Vp polarizations
% reference axes - text labels
ref_axis_1 ='a';
ref_axis_2 ='b';
ref_axis_3 ='c';
plot_name = '/MatLab_Programs/Single_Crystal_Enstatite_Seismic_plot.png';
Func_Seismic_multiplot_gfs_MTEX4(C_opx,rho_opx,plot_name,plot_x,...
Vp_polar,ref_axis_1,ref_axis_2,ref_axis_3,ref_vec_1,ref_vec_2,ref_vec_3,HS,global_fontSize)
%%
%************************************************
% Define elastic stiffness tensor of Diopside
%************************************************
%
% Reference Elastic constants
% Diopside Isaak,Ohno and Lee 2005 Phys.Chem.MIn.
%
% Reference Crystal Structure
% Diopside Di93 He3 Ur2 X2  C2/c XIIa*,YIIb,ZIIc
%
% Define density (g/cm3)
rho_cpx = 3.2860;
%
% Define Cartesian Tensor crystal symmetry & frame
% N.B. Changed new point group representation '12/m1' replaces '2/m'
cs_Tensor_cpx = crystalSymmetry('12/m1',[  9.7456  8.9198  5.2516],...
 [  90.0000 105.8600  90.0000]*degree,'x||a*','z||c',...
 'mineral','Diopside');
%
% elastic Cij stiffness tensor (GPa) as matrix M
 M =.... 
 [[  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]];
%
% M as stiffness tensor C with MTEX tensor command
C_cpx = tensor(M,cs_Tensor_cpx);
%
% Plot and Print out results
%
% print
elastic_tensor = 'Single crystal Diopside';
print_file = '/MatLab_Programs/Single_Crystal_Diopside_Seismic_Analysis.txt';
Func_Print_Seismic_Cij_MTEX4(C_cpx,rho_cpx,elastic_tensor,print_file)
% plot
close all

plot_x = 1; % plot X to North
Vp_polar = 1; % plot Vp polarizations
% reference axes - text labels
ref_axis_1 ='a*';
ref_axis_2 ='b';
ref_axis_3 ='c';
plot_name = '/MatLab_Programs/Single_Crystal_Diopside_Seismic_plot.png';
Func_Seismic_multiplot_gfs_MTEX4(C_cpx,rho_cpx,plot_name,plot_x,...
Vp_polar,ref_axis_1,ref_axis_2,ref_axis_3,ref_vec_1,ref_vec_2,ref_vec_3,HS,global_fontSize)
%%
%***********************************************************************
% Average tensors for single phases from corrected EBSD data
%***********************************************************************
%
%              NORTH
%                Y
%                |       
%         ----------------
%         |AGGREGATE  Vp |
%   WEST  |PHYSICAL   Vs1|--> X  
%         |PROPERTIES Vs2|
%         ----------------
%             SOUTH
%%
%**************************************************************************
% ploting convention for sample direction X to EAST
plot_x = 2; % plot X to EAST
% lower hemisphere plot MUST compatable with Pole figure plots
HS = 'lower';
Vp_polar = 1; % plot Vp polarizations
% Structural reference axes - text labels should be same as pole figure plots
% define Vector3d for reference directions as for pole figures
ref_vec_1 = xvector
ref_vec_2 = yvector
ref_vec_3 = -zvector
% labels
ref_axis_1 ='Z';
ref_axis_2 ='Y';
ref_axis_3 ='X';
% set colour map to seismic color map : blue2redColorMap
setMTEXpref('defaultColorMap',blue2redColorMap);
%**************************************************************************
%
% 100% Olivine aggregate : Voigt,Hill and Reuss averages
%
fprintf(' \n')
fprintf(' 100%% Olivine : Voigt, Hill and Reuss averages \n')
fprintf(' \n')
% calculate ODF using Dirichlet Kernel Band-width = Lmax = 4 for elastic properties (tensor 4th rank)
K_BW4 = DirichletKernel(4);
olivine_ODF_Dirichlet = calcODF(ebsd('Forsterite').orientations,'kernel',K_BW4);
%
[C_Voigt_Ol,C_Reuss_Ol,C_Hill_Ol] = calcTensor(olivine_ODF_Dirichlet,C_olivine)
%
% Plot and Print out results
%
% Hill average
% print
elastic_tensor = '100%% Olivine : Hill average';
print_file = '/MatLab_Programs/100_Olivine_Hill_average_Seismic_Analysis.txt';
Func_Print_Seismic_Cij_MTEX4(C_Hill_Ol,rho_olivine,elastic_tensor,print_file)
% plot
close all
plot_name = '/MatLab_Programs/100_Olivine_Hill_average_Seismic_plot.png';
Func_Seismic_multiplot_gfs_MTEX4(C_Hill_Ol,rho_olivine,plot_name,plot_x,...
Vp_polar,ref_axis_1,ref_axis_2,ref_axis_3,ref_vec_1,ref_vec_2,ref_vec_3,HS,global_fontSize)
%%
%
% 100% Enstatite : Voigt,Reuss, Hill averages
%
fprintf(' \n')
fprintf(' 100%% Enstatite : Voigt, Hill and Reuss averages \n')
fprintf(' \n')
% calculate ODF using Dirichlet Kernel Band-width = Lmax = 4 for elastic properties (tensor 4th rank)
K_BW4 = DirichletKernel(4);
enstatite_ODF_Dirichlet = calcODF(ebsd('Enstatite').orientations,'kernel',K_BW4);
%
[C_Voigt_Opx,C_Reuss_Opx,C_Hill_Opx] = calcTensor(enstatite_ODF_Dirichlet,C_opx)
%
% Plot and Print out results
%
% Hill average
% print
elastic_tensor = '100%% Enstatite : Hill average';
print_file = '/MatLab_Programs/100_Enstatite_Hill_average_Seismic_Analysis.txt';
Func_Print_Seismic_Cij_MTEX4(C_Hill_Opx,rho_opx,elastic_tensor,print_file)
% plot
close all
plot_name = '/MatLab_Programs/100_Enstatite_Hill_average_Seismic_plot.png';
Func_Seismic_multiplot_gfs_MTEX4(C_Hill_Opx,rho_opx,plot_name,plot_x,...
Vp_polar,ref_axis_1,ref_axis_2,ref_axis_3,ref_vec_1,ref_vec_2,ref_vec_3,HS,global_fontSize)
%%
% 100% Diopside : Voigt,Reuss, Hill averages
%  
fprintf(' \n')
fprintf(' 100%% Diopside : Voigt, Hill and Reuss averages \n')
fprintf(' \n')
% calculate ODF using Dirichlet Kernel Band-width = Lmax = 4 for elastic properties (tensor 4th rank)
K_BW4 = DirichletKernel(4)
diopside_ODF_Dirichlet = calcODF(ebsd('Diopside').orientations,'kernel',K_BW4)
%
[C_Voigt_Cpx,C_Reuss_Cpx,C_Hill_Cpx] = calcTensor(diopside_ODF_Dirichlet,C_cpx)
%
% Plot and Print out results
% Hill average
% print
elastic_tensor = '100%% Diopside : Hill average';
print_file = '/MatLab_Programs/100_Diopside_Hill_average_Seismic_Analysis.txt';
Func_Print_Seismic_Cij_MTEX4(C_Hill_Cpx,rho_cpx,elastic_tensor,print_file)
% plot
close all
plot_name = '/MatLab_Programs/100_Diopside_Hill_average_Seismic_plot.png';
Func_Seismic_multiplot_gfs_MTEX4(C_Hill_Cpx,rho_cpx,plot_name,plot_x,...
Vp_polar,ref_axis_1,ref_axis_2,ref_axis_3,ref_vec_1,ref_vec_2,ref_vec_3,HS,global_fontSize)
%
%***********************************************************************
% Compute the aggregate density
%***********************************************************************
% Define density (g/cm3) of aggregate
rho_aggregate = af_Forsterite*rho_olivine + af_Enstatite*rho_opx + af_Diopside*rho_cpx;
fprintf(' \n')
fprintf(' Aggregate density (g/cm3) =  %8.4f \n', rho_aggregate)
fprintf(' \n')
%%
%***********************************************************************
% Average tensors for all phases from corrected EBSD data
%***********************************************************************
% Olivine, Enstatite and Diopside aggregate: Hill average
fprintf(' \n')
fprintf(' Aggregate (Olivine,Enstatite,Diopside) :   Hill average \n')
fprintf(' \n')
% Hill average
C_Hill_Ag = af_Forsterite*C_Hill_Ol + af_Enstatite*C_Hill_Opx ...
               + af_Diopside*C_Hill_Cpx
% Hill average
% print
elastic_tensor = 'Aggregate (Olivine,Enstatite,Diopside) : Hill average';
print_file = '/MatLab_Programs/Aggregate_Hill_average_Seismic_Analysis.txt';
Func_Print_Seismic_Cij_MTEX4(C_Hill_Ag,rho_aggregate,elastic_tensor,print_file)
% plot
close all
plot_name = '/MatLab_Programs/Aggregate_Hill_average_Seismic_plot.png';
Func_Seismic_multiplot_gfs_MTEX4(C_Hill_Ag,rho_aggregate,plot_name,plot_x,...
Vp_polar,ref_axis_1,ref_axis_2,ref_axis_3,ref_vec_1,ref_vec_2,ref_vec_3,HS,global_fontSize)
%% Default MTEX colormap - jet colormap begin with white
setMTEXpref('defaultColorMap',WhiteJetColorMap);
%
disp(' ')
disp(' Script has successfully terminated')
disp(' ')
%
% End of script
%