%% Physical Properties using MTEX
%%
%% Luiz F. G. Morales - ScopeM / SGT group - ETH Zürich
% luiz.morales@scopem.ethz.ch
%
% based on MTEX development by David Mainprice (Geosciences Montpellier),
% Ralf Hielscher, Helmut Schaeben (Universität Freiberg) and Florian Bachmann
%
%
%% A potential plan for this lecture
% Let’s start easy – tensor basics, single crystal 2nd rank tensor
% 3rd-rank tensors - Piezoelectricity
% Single crystal elastic 4th rank tensors
% Averaging elastic properties for aggregates – Voigt, Reuss and Hill averaging schemes
%
%% And what is anisotropy?
% Anisotropy is a direction-dependent response of the crystal: the same stimulus applied in different directions (i.e. a tensile stress) leads to different responses (elastic strain);
% The tensors that describe the stimulus-response relation in crystals (elastic, thermal, plastic, electric) are usually anisotropic;
% Single crystals are always anisotropic, and the anisotropy reflects the symmetry of the crystal lattice. For example, while the number of independent components of the elastic tensor is two for an isotropic medium, it is:
% three for cubic symmetry;
%  five for hexagonal symmetry
% nine for orthogonal symmetry
%%
%
%% Optical anisotropy of calcite
%
% Calcite optical properties : 2nd Rank Tensor
%
% When light passes through a calcite crystal it splits into two rays, extra ordinary and an ordinary ray. Ordinary ray travels straight and extra ordinary ray is refracted.
%
%% Anisotropy
% If a property is anisotropic in the grains and if the grains have preferred orientation, than the aggregate will also be anisotropic;
% The symmetry of the property is the same as the symmetry of the texture:
% an aggregate with randomly distributed orientations will have isotropic properties
% Depending on the phases evolved and the type of deformation, anisotropy may have axial, orthorrhombic or monoclinic symmetry;
%
%% What is a tensor?
% Tensors are mathematical objects used to describe physical properties, just like scalars and vectors;
% Tensors are merely a generalisationof scalars and vectors; a scalar is a zero rank tensor, and a vector is a first rank tensor;
% The rank (or order) of a tensor is defined by the number of directions required to describe it;
% For example, properties that require one direction (first rank) can be fully described by a 3×1 column vector, and properties that require two directions (second rank tensors), can be described by 9 numbers, as a 3×3 matrix;
%
%% Tensor ranking of physical properties
%
%
%% Second-rank tensors
% Typically relates 2 vectors – for example thermal conductivity :  applied vector (negative) temperature gradient  and resulting vector heat flow density, exception thermal expansion relates temperature (0) and strain (2).
%
% The generic 2nd rank tensor T is the relation between an applied vector p and resultant vector q. We can write relation between p and qas a tensor equation
%                                     p= T q  or  pi = Tij qj  (i=1,2,3 ; j=1,2,3)
% In general the vectors p and q are not parallel.
%
%% Radius-normal property and melting wax experiment
% Consider a red hot arrow touching the second order prism plane of a single crystal of quartz and creating a point heat source. A (negative) thermal gradient will result in heat flow away from the hot heat source (arrow tip) to the colder regions of the crystal. As the thermal resistivity (reciprocal of thermal conductivity) varies with direction the heat flow will not be equal in all radial directions, i.e. heat flow is anisotropic.
%
%
%
%% How to define a tensor in MTEX – thermal conductivity

%% Define tensor reference frame

cs_tensor_ort = crystalSymmetry('12/m1',[8.561, 12.996, 7.192],...
    [90, 116.01, 90]*degree, 'mineral','orthoclase','Y||b','Z||c') 

%% Define thermal conductivity tensor coefficients (in W/m/K unit)

M = [[1.45 0.00 0.19];[0.00 2.11 0.00];[0.19 0.00 1.79]]

%%  Define the thermal conductivity tensor

k = tensor (M, 'name', 'thermal_conductivity','rank',2, 'unit', ...
    'W_1/m_1/K',cs_tensor_ort)

%
%% How to represent tensor properties in MTEX – thermal conductivity



%% Plot thermal conductivity tensor

plot (k, 'complete','lower','contourf')
mtexColorbar

% annotate with single crystal pole to planes

hold on
text(Miller({1,0,0},{0,1,0},{0,0,1},cs_tensor_ort), {'(100)','(010)',...
    '(001)'},'backgroundColor','w','fontweight','bold','fontsize',18)
hold off

% don't forget to save your figures

saveFigure ('thermal conductivity tensor orthoclase.png')

%
%% A very important property of tensors
% One of the fundamental properties of tensors of physical properties is that the magnitude of the physical property is invariant with the rotation of the coordinate reference frame, which is attached to the crystal;
% A crystal may rotate in sample reference frame (which is an external reference frame), but the tensor reference frame is related to the crystal structure and rotates with the crystal, like a rigid rotation
% What happens if we rotate the tensor previously defined?
%
%% Rotate tensor and replot

%% Let's rotate the tensor and plot it again

% define a rotation

g = rotation.byEuler(10*degree,20*degree,30*degree)

% rotate tensor

k_rot = rotate (k,g)

% plot rotated tensor 

plot (k_rot, 'complete','lower','contourf','minmax')
mtexColorbar

% annotate with single crystal pole to planes

hold on
text(g*Miller({1,0,0},{0,1,0},{0,0,1},cs_tensor_ort),{'(100)','(010)',...
    '(001)'},'backgroundColor','w','fontweight','bold','fontsize',18)
hold off

% don't forget to save your figures

saveFigure ('thermal conductivity tensor orthoclase rotated.png')

%
%% The quadric representation for 2nd-rank tensors
%
% A common way to represent second-rank tensors is via a second degree surface (also named a "quadric");
%  The quadric may be an ellipsoid or a hyperboloid.
% Most common second-rank tensors are symmetric (Tij = Tji). In that case, and if the 3 main coefficients are positive, then the property can be represented by an ellipsoid with 1/√T1, 1/√T2 and 1/√T3,
% That is the case of thermal conductivity, electric polarization and optical properties
% In MTEX, this elipsoid can be plotted with the command 'surf'
%

%% Plot thermal conductivity as 3D surface

% plot quadric representation for thermal conductivity

newMtexFigure('layout',[1,2]) % figure layout: one line, two columns
surf(k.directionalMagnitude)

% plot section along X on the surface
nextAxis 
plotSection(k.directionalMagnitude,vector3d.X)
drawNow(gcm)
saveFigure ('thermal conductivity tensor orthoclase surface.png')

% In this example, the section is parallel to X vector of the elipsoid.
%  Note that the quadric representation can be also rotated on the fly for the visualization in different directions.
%
%% The hyperboloid representation for 2nd-rank tensors
% When one of the main diagonal coefficients of the tensor is negative, the surface is a hyperboloid of one sheet (e.g. thermal expansion of plagioclase);
% If two values are negative, then the surface will be a hyperboloid of two sheets (e.g. thermal expansion of calcite - contraction on the basal plane)
% If all three values are negative, as in the case of magnetic suceptibility of some minerals, the surface is a imaginary ellipsoid
%
%% The example of thermal expansion of calcite

%% Define tensor reference frame

cs_tensor_cal = crystalSymmetry('321', [4.99 4.99 17.064], 'X||a', ...
    'Y||b*', 'Z||c','mineral', 'Calcite', 'color', 'light blue');

%% Define thermal expansion coefficient

%  Coefficients from Fei 1995 AGU Bookshelf Series Mineral 
%  Physics vol.2 Chapter 6 Thermal expansion in x 10-6 1/K units

M_cal = [[-3.2E-6 0.00 0.00];[0.00 -3.2E-6 0.00];[0.00 0.00 13.3E-6]]
mte
%%  Define the thermal expansion tensor

k_cal = tensor (M_cal, 'name', 'thermal_expansion','rank',2, 'unit', ...
    '[K^-1]',cs_tensor_cal)

% plot the thermal expansion tensor 

figure (1)
plot (k_cal, 'complete','lower','contourf')
mtexColorbar
h = Miller({1,0,-1,0},{1,1,-2,0},{1,0,-1,4},{0,1,-1,0},{0,0,0,1},cs_tensor_cal);
h_s=symmetrise(h)
hold on
annotate(h_s,'labeled','backgroundColor','w','grid','upper','doNotDraw')
hold off
drawNow(gcm,'figSize','large')
saveFigure('calcite thermal expansion.png')

% plot the thermal expansion tensor surface (hyperboloid)

figure (2)
surf(k_cal.directionalMagnitude)
saveFigure('calcite thermal expansion hyperboloid.png')

%
%
%% True hyperboloid


clear all;
[X,Y,Z] = meshgrid(-10:0.5:10,-10:0.5:10,-10:0.5:10);
a=1;
b=1;
c=1;
V = X.^2/a^2 + Y.^2/b^2 - Z.^2/c^2;
p=patch(isosurface(X,Y,Z,V,1)); % This is the key step. It involves getting the part of the volume corresponding to the surface defined by the equation
set(p,'FaceColor','red','EdgeColor','none');
daspect([1 1 1])
view(3);
camlight

%
%% The effect of symmetry on 2nd rank tensors
%
%
%
%% Neumann's principle - effect of symmetry in physical properties of crystals
%
% In 1885, F.E. Neumann stated that that the “symmetry elements of any physical property of a crystal must include ALL the symmetry elements of the point group of the crystal”
% This implies that certain physical properties of a given material may have higher symmetry than the material itself, but will never be less symmetric than it.
% For example: optical properties in cubic materials are isotropic
% Many physical properties are inherently centrosymmetric (e.g. all symmetric second order tensors and elasticity)
% This "adds" a center of symmetry in many minerals (e.g. quartz) and result in a higher symmetry than the possessed by the crystal.
%
%
%
%% Piezoelectricity
%
% Eletrical charge that accumulates in certain solids in response to applied stress
% Reported by the first time in 1880 by Jacques and Pierre Currie, while doing experiments in a series of minerals, including tourmaline, spharelite, alpha-quartz (between others)
% Some geological materials, some ceramics and even biological matter (e.g. bones, some proteins) have piezoeletric properties
% Not only a scientific "curiosity", but very important in the industry, particularly transducers and resonators stems
% Piezoeletricitiy is reversible - removing the stress will remove the induced eletrical field
% The semi-conducting elements tellurium (Te) and selenium (Se), along with the mineral pyrolusite (MnO2), have the greatest piezoelectric effect of naturally occurring compounds.
% Although piezoelectric minerals are generally not very common in the Earth’s crust, one is very common: quartz is the third-most common mineral at about 12% of the Earth’s crust, according to Taylor & McLennan (1985).
%
%% Direct effect - applying stress to generate electric current
%
%
%
% Figure source: https://www.autodesk.com/products/eagle/blog/piezoelectricity/
%
%% Converse effect - apply current to produce mechanical energy (release in form of sound waves)
%
%
%
%
% Figure source: https://www.autodesk.com/products/eagle/blog/piezoelectricity/
%
%% Piezoeletricity in quartz
%
% Figure source: https://www.autodesk.com/products/eagle/blog/piezoelectricity/
%
%% Some basics on piezoeletricity 3rd-rank tensor
% The direct effect can be written as the relationship between the 2nd rank stress tensor  rank electric polarization vector Pi, linked by the piezoelectric tensor dijk as follows
%
% This can be written in a more compact tensor notation:
%
% In this equation Pi is the electric polarization, dijk is the piezoeletric tensor and  sigma jk is the stress tensor
% The piezoelectric tensor will always have three indices and Voigt matrix notation two indices. The direct and converse effects can also be written in Voigt matrix notation as
%
%
%% The piezoelectric tensor
% 1st rank electric polarization (Pi) and electric field (Ek) vectors have the index i=1,2 or 3
% 2nd rank stress \sigma \textrm{ij} and elastic strain \varepsilon \textrm{jk}also have indices j=1,2,3 and k=1,2,3, sp the piezoelectric tensor dijk has 3 x 3 x 3 (=27) coefficients
% The stress and strain 2nd rank tensors  are symmetric for linear elasticity, resultingin interchangeability of jk and kj indices of piezoelectric tensor dijk
% That reduces the number of independent components to 18, im the case of triclinic crystals
% This number is considerably reduced in higher crystalline symmetries. In the case of monoclinic crystals that have a 2-fold rotation of 180° around [b]or-[c]-axis, the number of independent components is reduced to 8, and can be 1 in the case of hexagonal -62 and cubic -43m
%
%% Symmetry properties related to piezoelectricity
% Crystals that belong to centrosymmetric point groups (i.e. the 11 Laue classes) are piezoelectrically inactive
% All crystals belonging to the 21 noncentrosymmetric point groups are piezoelectrically active, with the exception of cubic 432
% In the 432 point group, the presence of four-fold axes parallel to [100], [010] and [001] makes all directions perpendicular to these axes non-polar
% Because of that, all the tensor oefficients that are non-zero in two other piezoeletric cubic point groups 23 and -43m (d14=d25=d36), are zero in 432
% That results in 20 non-centrosymmetric point groups for piezoelectric active crystals, where 10 enantiomorphic space group pairs do not have improper rotations (i.e. no mirror planes).
% Such crystals occur in right-handed and left-handed forms (e.g. a-quartz) and are distributed in 6 point groups (4, 422, 3, 32, 6 and 622)
%
%% Loading quartz piezoelectricity in MTEX
%

% Define quartz tensor reference frame

cs_tensor_quartz = crystalSymmetry('32', [4.916 4.916 5.4054],...
    'X||a*', 'Z||c', 'mineral', 'Quartz');

% Load the tensor (this one can be loaded automatically from MTEX)

fname = fullfile(mtexDataPath,'tensor', 'Single_RH_quartz_poly.P');

P = tensor.load(fname,cs_tensor_quartz,'propertyname',...
    'piecoelectricity','unit','C/N','DoubleConvention') 

%
%% Making some colorful plots

% Let's change the colorcode for the plots

setMTEXpref('defaultColorMap',blue2redColorMap);

% Now let's plot

plot(P)
mtexColorbar
saveFigure('quartz piezoelectricity 120 sector.png')

% This plot shows thew magnitude polarization variation, from red (strongest) to blue (weakest). As a standard, we use the unique region implied by crystal symmetry, that in the case of quartz (trigonal) is represented by a 120° sector of the full pole figure
%
% One naturally can also plot the complete crystal behavior

% Plot quartz piezoelectricity (full circle)

plot(P,'complete','smooth','upper')
h = Miller({0,0,0,1},{1,1,-2,0},{2,-1,-1,0},{1,0,-1,0},{1,0,-1,1},...
    {0,1,-1,1},{1,0,-1,2},{0,1,-1,2},cs_tensor_quartz);
h_s=symmetrise(h)
hold on
annotate(h_s,'labeled','backgroundColor','w','grid','upper','doNotDraw')
hold off
drawNow(gcm,'figSize','large')
mtexColorbar
saveFigure('quartz piezoelectricity complete.png')

% As in the other examples, one can plot single crystal piezoelectricity as a surface

% Plot quartz piezoelectricity as a surface

surf(P.directionalMagnitude)
saveFigure('quartz piezoelectricity surface.png')

% Plot quartz piezoelectricity section along basal plane
nextAxis 
plotSection(P.directionalMagnitude,vector3d.Z)
axis on
drawNow(gcm)
saveFigure ('section quartz piezoelectricity Z.png')

%%
%
%% CPO-derived seismic properties - elasticity and 4th-rank tensor
%% Why it is important?
%
% •The calculation of the seismic properties from microstructural information (crystal orientation, volume fraction, grain shape etc.) is important for rocks because it gives insight into the role of microstructure in determining the bulk properties;
% •Can be calculated for the conditions where the microstructures were developed via using the appropriate single crystal derivatives;
% •Additional features not necessarily preserved in the recovered microstructure, such as the presence of fluids (e.g. magma) can be modeled;
% •The effect of phase change on the seismic properties can also be modeled using these methods;
% •Modeling is essential for understanding anisotropic properties as experimental measurements in the many directions necessary to fully characterize anisotropy is not currently feasible for the majority of the temperature and pressure conditions found in the deep Earth.
%
%
%% What causes seismic anisotropy in rocks?
%
%
%% Shear wave splitting
% https://en.wikipedia.org/wiki/Shear_wave_splitting#/media/File:Garnero_Splits.gif
%
%
%% Elasticity - a 4th-rank tensor
%
% In any crystalline material there is balance between Coulomb attractive forces between oppositely charge ions and Born repulsive forces due to the overlap of electron shells. At any given thermodynamic state the crystal will tend toward an equilibrium structure;
%
% For a change in hydrostatic or non-hydrostatic stress, the crystal structure will adjust at the atomic level to the new thermodynamic state.
% The fundamental nature of atomic forces in the determination of elastic properties has been illustrated by the emergence of first principles atomic modeling to predict single crystal elastic tensors of geophysical  importance at lower mantle PT conditions.
%
%% If the stress and strain tensor were not symmetric
%
%
%% Stress and strain tensor
%
% That results in 36 independent components
%
%% Crystal symmetry plays a important role,
%% In the case of less symmetric materials (triclinic) - there will be...
%%
%
%
%% Tensor in a cubic material
%
% There are a maximum of 21 elastic constants for a crystalline body, but for cubic crystals the elastic constants, cij may be reduced to just three independent elastic constants:
% c11= c22 = c33 ⇒ modulus for axial compression, i.e., a stress σ11results in a strain ε11along an axis;
% c44 = c55 = c66 ⇒ shear modulus, i.e., a shear stress σ23results in a shear strain ε23 across a face;
% c12 = c13 = c23 ⇒ modulus for dilation on compression, i.e., an axial stress σ11results in a strain ε22 along a perpendicular axis.
% All other cij = 0.
%
%% Tensors and the crystalline symmetry
%
%
% from Almqvist & Mainprice 2017
%
%% What can you do with elastic constants?
% Volume compressibility
% Linear compressibility
% Young's modulus
% Shear Modulus
% Poisson's ratio
% Christoffel tensor - seismic velocities
% etc...
%
%% Let's upload the olivine single crystal elastic tensor
%%
%% First we define the single crystal tensor reference frame (olivine)

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

%% Then we define the density of the material (g/cm3)

rho_olivine = 3.3550;

%% And the coefficients Cij for olivine in Voigt notation (Abramson et al., 1997 - JGR)

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

%% To define the stiffness tensor, you will use the command stiffnessTensor

C_olivine = stiffnessTensor(Cij,cs_tensor_olivine,'density',...
    rho_olivine);

%
%% Now you can plot a certain number of properties derived from the stiffness tensor.
%% For instance, Young's modulus

 setMTEXpref('defaultColorMap',blue2redColorMap);
 figure(1)
 plot(C_olivine.YoungsModulus,'complete','upper','contourf');
 mtexColorbar
 h = Miller({1,0,0},{0,1,0},{0,0,1},cs_tensor_olivine,'uvw');
 hold on
 annotate(h,'labeled','backgroundColor','w','grid','upper',...
     'doNotDraw')
 hold off
 saveFigure ('olivine single crystal youngs modulus.png')

%% Or the linear compressibility

figure(2)
 plot(C_olivine.linearCompressibility,'complete','upper','contourf');
 mtexColorbar
 h = Miller({1,0,0},{0,1,0},{0,0,1},cs_tensor_olivine,'uvw');
 hold on
 annotate(h,'labeled','backgroundColor','w','grid','upper',...
     'doNotDraw')
 hold off
 saveFigure ('olivine single crystal linear compressibility.png')
 

%
%% The single crystal seismic velocities are computed with command seismicVelocities

plotSeismicVelocities(C_olivine)

%
%% Calculating CPO-derived seismic properties
%% First, we have to calculate an ODF
% In texture analysis, the orientation of a crystal is defined historically by the letter g, where=(φ1,ϕ, φ2) - an Euler angle triplet;
% The ODF f(g) is the volume fraction of orientations within an orientation range between g and g+dg, in a space containing all possible orientations given by
%
% Region of integration in orientation space – symmetrically controlled;
% The ODF is normalized and ODF densities are expressed in multiples of the uniform density (m.u.d.). The normalization for the ODF is given by
%
% For a uniform distribution the ODF density is 1 for all orientations in the ODF, defined by f(g)uniform=1
% The strength of preferred orientation of an ODF may be expressed by the ODF texture index
%
%% Cartesian Reference Frame
% Cartesian frame defined by orthogonal crystallographic directions used to report single crystal elastic properties may not be the same as those used for Euler angle reference frame (texture analyses);
% Crystallographic reference frame of tensor needs to be rotated to bring it into coincidence with Euler angle frame,
%
% Here Cijkl(gE) is the elastic property in the Euler reference and Cpqrt(gT) is the elastic property in the original tensor reference frame, BOTH  frames are in crystal  co-ordinates.
% To calculate the seismic properties of a polycrystal, one must evaluate the elastic properties of the aggregate;
% In an aggregate with CPO, the anisotropy of single crystal elastic properties need to be considered;
% For each orientation g, the single crystal properties have to be rotated into the specimen co-ordinate
%
% Here, Cijkl(g) is the elastic property in sample co-ordinates, gij is the crystal orientation in sample co-ordinates and Cpqrt(gE) is the elastic property in crystal co-ordinates of the Euler frame (E).
%
%% Elastic constant averaging schemes - Voigt and Reuss
%% Advantages
% easy to calculate
% widely separated bounds for highly anisotropic phases
%% Disadvantages
% cannot introduce microstructural effect (e.g. grain shape)
% very poor bounds for mixtures of materials with very different stiffness (e.g. solids + liquids + voids)
% In this cases, better to use differential effective media or self-consistent approaches (to be incorporated on MTEX)
%
%% Elastic properties of an aggregate
% The elastic properties of the polycrystal may be calculated by integration over all possible orientations of the ODF.
%
% where <Cijkl>m is the elastic properties of the aggregate of mineral m.
% Alternatively it may be determined by simple summation of individual orientation measurements (e.g. U-stage or EBSD),
%
% in the equation above, V(g) is the volume fraction of the grain in orientation g
%
%% Christoffel equation & Tensor
%
% The final step is the calculation of the three seismic phase velocities by the solution of the Christoffel equation,
%
% Here, V are the phase velocities, ρ is the density of the rock, nj is the propagation unit vector (propagation direction), Pk  is the polarization unit vectors (polarization direction), and Cijkl is the aggregate elastic constant;
% The polatization unit vectors Pk are obtained as eigenvectors and eigenvalues of the roots of the equation
%
%
%
%% Vp anisotropy of some minerals
%
%
%% Vp anisotropy of some mantle minerals
%
%
%% Calculating CPO-derived seismic properties
%% Let's import some data

mtexdata forsterite
plot (ebsd)

%
%% We need now to rotate the data, to have the foliation here in the horizontal position

ebsd = rotate(ebsd,rotation('axis',-zvector,'angle',90*degree));
plot(ebsd)

%
%% Now we import the single crystal elastic tensors we need
% The elastic stiffness tensor of Olivine was reported in Abramson et al., 1997 (Journal of Geophysical Research) with respect to the crystal reference frame

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

% we also need the density of this phase

rho_olivine = 3.3550;

% and the coefficients Cij, in Voigt matrix notation

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]];

%
%% In order to define the stiffness tensor as an MTEX variable we use the command stiffnessTensor.

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

%
%% ODF and polycrystalline elastic stiffness
% For large datasets, the calcTensor might be very slow, so it is better to computer the ODF of the phases before combining them with the single crystal elastic stiffness

odf_ol = calcDensity(ebsd('f').orientations,'halfwidth',...
    10*degree);

% then you combine with the single crystal elastic tensor

[CVoigt_ol, CReuss_ol, CHill_ol]    = mean(C_olivine,odf_ol)

%
%% And finally you can plot the seismic velocities of your aggregate
%

plotSeismicVelocities(CHill_ol)