%*****************************************************************************************
% Demo : Shear modulus for a dislocation slip system
%*****************************************************************************************
%
% David Mainprice (2/05/2019)
%
% Reference Elastic constants
% alpha-Quartz McSkimin et al 1965 J.Appl.Phys.v.36p.1624-1632
% Reference Crystal Structure
% alpha quartz Trigonal SG P3221 25C d=2.648 g/cm3
% Elastic tensor Frame is IRE setting  : x= a-axis z= c-axis
% Define density (g/cm3)
 rho= 2.6480;
%
% Define Cartesian Tensor crystal crystalSymmetry & frame
cs_Tensor = crystalSymmetry('321',[  4.9134  4.9134  5.4052],...
 [  90.0000  90.0000 120.0000]*degree,'x||a','z||c',...
 'mineral','Alpha-Quartz 1965');
%
% elastic Cij stiffness tensor (GPa) as matrix M
 M =.... 
 [[   86.80    7.04   11.91  -18.04    0.00    0.00];...
 [     7.04   86.80   11.91   18.04    0.00    0.00];...
 [    11.91   11.91  105.75    0.00    0.00    0.00];...
 [   -18.04   18.04    0.00   58.20    0.00    0.00];...
 [     0.00    0.00    0.00    0.00   58.20  -18.04];...
 [     0.00    0.00    0.00    0.00  -18.04   39.88]];
%
% M as stiffness tensor C with MTEX tensor command
C = tensor(M,cs_Tensor);

% Shear Modulus of dislocation slip system
% b= slip direction    (in specimen coordinates - vector3d)
% n= slip plane normal (in shear plane normal  - vector3d)
%
a1_direction = Miller(2,-1,-1,0,cs_Tensor,'UVTW')
c_plane = Miller(0,0,0,1,cs_Tensor,'hkil')
% convert crystallographic to vector3d
b=vector3d(a1_direction)
n=vector3d(c_plane)
% shear modulus in GPa
G_a1_c_slip = shearModulus(C,n,b)
