% Crystallography_In_MTEX.m
%
% Crystallography in MTEX (version tested with MTEX 4.0.6)
%
% crystal symmetry (CS) and specimen symmetry (SS)
%
% By David Mainprice  last revised 21/11/2014
%
%% clear memory and close all plots
clear all
close all
%% List point groups in MTEX4
pG = crystalSymmetry.pointGroups
strvcat(pG.Inter)
% In International notation
%
% Triclinic
%  1    
% -1 Laue class
%
% Monoclinic
% 211  
% m11  
% 2/m11 Laue class two-fold a-axis
% 121  
% 1m1  
% 12/m1 Laue class two-fold b-axis
% 112  
% 11m  
% 112/m Laue class two-fold c-axis
%
% Orthorhombic
% 222  
% 2mm  
% m2m  
% mm2  
% mmm  Laue class
%
% Trigonal
%  3    
% -3   
% 321  
% 3m1  
%-3m1 
% 312  
% -31m
%
% Tetragonal
% 4    
% -4   
% 4/m  Laue class
% 422  
% 4mm  
% -42m 
% -4m2 
% 4/mmm Laue class
%
% Hexagonal
%  6    
% -6   
% 6/m   Laue class
% 622  
% 6mm  
% -62m 
% -6m2 
% 6/mmm Laue class
%
% Cubic
% 23   
% m-3   Laue class
% 432  
% -43m 
% m-3m  Laue class
%%
%***********************************************************************
% specify specimen symmetry (SS)
%***********************************************************************
disp(' ')
disp('                  Specimen Symmetry (SS)')
disp('      3 Proper rotational point groups in MTEX')
disp('     Crystal system in CAPITALS for Laue classes ')
disp('------:--------------:------------:------:--------:-----------')
disp(' Code : crystal      : Schoen-  : Point  : Laue   : Rotational')
disp('      : system       : flies    : group  : class  : point group')
disp('------:--------------:----------:--------:--------:-----------')
disp('   1  : TRICLINIC    :   Ci     :  -1    :   -1   :     1')
disp('------:--------------:----------:--------:--------:-----------')
disp('   2  : MONOCLINIC*  :   C2h    :  2/m   :  2/m   :     2')
disp('------:--------------:----------:--------:--------:-----------')
disp('   3  : ORTHORHOMBIC :   D2h    :  mmm   :  mmm   :    222')
disp('------:--------------:----------:--------:--------:-----------')
disp(' ')
disp(' * monoclinic setting 2-fold axis is Y-axis')
disp(' ')
%%
%***********************************************************************
% specify crystal symmetry (CS)
%***********************************************************************
disp(' ')
disp('                  Crystal Symmetry (CS)')
disp('      32 Proper rotational point groups in MTEX')
disp('     Crystal system in CAPITALS for Laue classes ')
disp('------:--------------:------------:------:--------:-----------')
disp(' Code : crystal      : Schoen-  : Point  : Laue   : Rotational')
disp('      : system       : flies    : group  : class  : point group')
disp('------:--------------:----------:--------:--------:-----------')
disp('   1  : triclinic    :   C1     :   1    :   -1   :     1')
disp('   2  : TRICLINIC    :   Ci     :  -1    :   -1   :     1')
disp('------:--------------:----------:--------:--------:-----------')
disp('   3  : monoclinic*  :   C2     :   2    :  2/m   :     2')
disp('   4  : monoclinic*  :   Cs     :   m    :  2/m   :     2')
disp('   5  : MONOCLINIC*  :   C2h    :  2/m   :  2/m   :     2')
disp('------:--------------:----------:--------:--------:-----------')
disp('   6  : orthorhombic :   D2     :  222   :  mmm   :    222')
disp('   7  : orthorhombic :   C2v    :  mm2   :  mmm   :    222')
disp('   8  : ORTHORHOMBIC :   D2h    :  mmm   :  mmm   :    222')
disp('------:--------------:----------:--------:--------:-----------')
disp('   9  : tetragonal   :   C4     :   4    :  4/m   :     4')
disp('  10  : tetragonal   :   S4     :  -4    :  4/m   :     4')
disp('  11  : TETRAGONAL   :   C4h    :  4/m   :  4/m   :     4')
disp('------:--------------:----------:--------:--------:-----------')
disp('  12  : tetragonal   :   D4     :  422   : 4/mmm  :    422')
disp('  13  : tetragonal   :   C4v    :  4mm   : 4/mmm  :    422')
disp('  14  : tetragonal   :   D2d    :  -42m  : 4/mmm  :    422')
disp('  15  : TETRAGONAL   :   D4h    :  4/mmm : 4/mmm  :    422')
disp('------:--------------:----------:--------:--------:-----------')
disp('  16  : trigonal     :   C3     :   3    :   -3   :     3')
disp('  17  : TRIGONAL     :   C3i    :  -3    :   -3   :     3')
disp('------:--------------:----------:--------:--------:-----------')
disp('  18  : trigonal     :   D3     :  32    :   -3m  :    32')
disp('  19  : trigonal     :   C3v    :  3m    :   -3m  :    32')
disp('  20  : TRIGONAL     :   D3d    : -3m    :   -3m  :    32')
disp('------:--------------:----------:--------:--------:-----------')
disp('  21  : hexagonal    :   C6     :   6    :   6/m  :     6')
disp('  22  : hexagonal    :   C3h    :  -6    :   6/m  :     6')
disp('  23  : HEXAGONAL    :   C6h    :  6/m   :   6/m  :     6')
disp('------:--------------:----------:--------:--------:-----------')
disp('  24  : hexagonal    :   D6     :  622   : 6/mmm  :    622')
disp('  25  : hexagonal    :   C6v    :  6mm   : 6/mmm  :    622')
disp('  26  : hexagonal+   :   D3h    : -6m2   : 6/mmm  :    622')
disp('  27  : HEXAGONAL    :   D6h    : 6/mmm  : 6/mmm  :    622')
disp('------:--------------:----------:--------:--------:-----------')
disp('  28  : cubic        :   T      :  23    :  m-3   :     23')
disp('  29  : CUBIC        :   Th     :  m-3   :  m-3   :     23')
disp('------:--------------:----------:--------:--------:-----------')
disp('  30  : cubic        :   O      :  432   :  m-3m  :    432')
disp('  31  : cubic        :   Td     :  -43m  :  m-3m  :    432')
disp('  32  : CUBIC        :   Oh     :  m-3m  :  m-3m  :    432')
disp('------:--------------:----------:--------:--------:-----------')
disp(' ')
disp(' * monoclinic setting 2-fold axis is b-axis')
disp(' ')
disp(' + hexagonal point group -6m2 can also written as -62m')
disp(' ')
%%
%***********************************************************************
% specify specimen symmetry (SS)
%***********************************************************************
% N.B.
% You do not need to specify triclinic ('-1') specimen symmetry in MTEX 4
% as it is the DEFAULT specimen symmetry
% specimen symmetry , syntax is 
SS = specimenSymmetry('-1');
%%
%***********************************************************************
% specify crystal symmetry (CS) using a *.cif
%***********************************************************************
% Defining a Crystal Symmetry by a CIF File 
% MTEX allows to define a crystal symmetry by
% importing a crystallographic information file (*.cif)
% these files are stored in the MTEX folder
% MatLab_Programs/mtex-4.0.6/data/cif
% The exact address will depend on your PC
%
%% Alpha-quartz : file quartz.cif
CS = loadCIF('quartz')
%% Forsterite (Setting Pbnm used in Mineralogy) : file Forsterite1973.cif
CS = loadCIF('Forsterite1973')
%% Fe rich Olivine (Setting Pnma used in Physics) : file olivin.cif
CS = loadCIF('olivin')
%%
%**************************************************************************
% Crystal Symmetry (CS)
% Example of Triclinic calcium rich plagioclase Anorthite (An100)
%**************************************************************************
% 1. crystal symmetry system
% 2. [a,b,c] - unit cell lengths (typically in angstroms 10^-10m = 0.1 nm)
% 3. [alpha,beta,gamma]*degree -  unit cell angles in degrees
% 'X||a*','Z||c' - Two orthogonal reference directions from X,Y or Z,
% which are  parallel to two orthogonal crystallographic directions
% from a,a*,b,b* and c,c*
% For example in the Oxford/HKL EBSD system for Euler angles the
% reference frame is a* and c
% 4. mineral name 'An100'
%
% plagioclase feldspar anorthite CaAl2Si2O8 An100
% Wainwright and Starkey 1971
CS = crystalSymmetry('-1',[8.173,12.869,14.165],[93.11,115.91,91.26]*degree,...
  'X||a*','Z||c','mineral','An100')
%% plot convention x-axis to north
plotx2north
%% Plotting symmetry elements - 3 two fold axes at 90 degrees - > mirrors
CS = crystalSymmetry('mmm')
plot(CS)
%% Plotting symmetry elements
CS = crystalSymmetry('12/m1')
plot(CS)
%% Plotting symmetry elements
CS = crystalSymmetry('2/m')
plot(CS)
%% Plotting symmetry elements - three fold parallel to c-axis
CS = crystalSymmetry('3')
plot(CS)
%% Plotting symmetry elements - add an inversion
CS = crystalSymmetry('-3')
plot(CS)
%% Plotting symmetry elements - add a mirror plane
CS = crystalSymmetry('-3m')
plot(CS)
%% Plotting equal angle stereograms 
CS = crystalSymmetry('m-3m');
% reset default fontsize for all following plots !
%setMTEXpref('FontSize',20)
plotHKL(CS,'projection','eangle','upper','grid_res',10*degree,'BackGroundColor','w','figSize','large','FontSize',20)
%% Plotting equal area stereograms 
CS = crystalSymmetry('m-3m');
plotHKL(CS,'projection','earea','upper','grid_res',10*degree,'BackGroundColor','w','figSize','large','FontSize',20)
%% plagioclase feldspar anorthite CaAl2Si2O8 An100 (Wainwright and Starkey
% (1971)
CS = crystalSymmetry('-1',[8.173,12.869,14.165],[93.11,115.91,91.26]*degree,...
  'X||a*','Z||c','mineral','An100')
% Plotting equal angle stereograms 
plotHKL(CS,'projection','eangle','upper','grid_res',10*degree,'BackGroundColor','w','figSize','large','FontSize',20)
%%
%***********************************************************************
% Define a crystal directions
%***********************************************************************
% pole to plane hkl h =  h.a* + k.b* + l.c*
h = Miller(1,1,1,CS,'An100','hkl')
% direction u = u.a + v.b + w.c
u = Miller(1,1,1,CS,'An100','uvw') 
%% magnitudes of h and u
% magnitude of hkl
mag_h = norm(h)
% d-spacing of hkl  d_hkl = 1/|ha* + kb* + lc*|
dspace = 1/norm(h)
% alternative syntax in MTEX
dspace_h = dspacing(h)
% magnitude of a direction
mag_u = norm(u)
%% angle between h and u
angle_between_h_and_u = angle(h,u)/degree
% plot h - hkls and u - uvws on separate plots
plot(h,'labeled')
hold on
plot(u,'labeled')
hold off
%% volume of unit cell
% a,b,c axes in Euclidian frame
%a = get(CS,'axis') ** check results **
a = CS.axes
% cell volume V = a.(b x c)
V  = dot(a(1),cross(a(2),a(3)))
% cell volume V = b.(c x a)
V2  = dot(a(2),cross(a(3),a(1)))
% cell volume V = c.(a x b)
V3  = dot(a(3),cross(a(1),a(2)))
%% reciprocal lattice vectors
% are defined such that 
% a* is normal to b and c
% b* is normal to c and a
% c* is normal to a and b
% a* = (b x c)/V
a_star = cross(a(2),a(3)) ./V
% b* = (c x a)/V
b_star = cross(a(3),a(1)) ./V
% c* = (a x b)/V
c_star = cross(a(1),a(2)) ./V
%% recoprical lattice volume
% cell volume V* = a*.(b* x c*)
V_star  = dot(a_star,cross(b_star,c_star))
% cell volume V* = b*.(c* x a*)
V2_star  = dot(b_star,cross(c_star,a_star))
% cell volume V* = c*.(a* x b*)
V3_star  = dot(c_star,cross(a_star,b_star))
%% direct lattice vectors
% a = (b* x c*)/V*
a_direct = cross(b_star,c_star) ./V_star
% b = (c* x a*)/V*
b_direct = cross(c_star,a_star) ./V_star
% c = (a* x b*)/V*
c_direct = cross(a_star,b_star) ./V_star
%% relationship between direction and reciprocal lattices
% like
% a*.a  where a* =  ha* + kb* + lc* = ha* = a* with h=1
%         and a  =  ua  + vb  + wc  = ua  = a  with u=1
% a*.a = |a*|.|a| cos(theta) = 1
% fixes the modulus of a*.a = 1 so that |a*| = 1/|a|  or |a| = 1/|a*| etc
% a*.a
dot(a_star,a_direct)
% b*.b
dot(b_star,b_direct)
% c*.c
dot(c_star,c_direct)
%% relationship between direction and reciprical lattices
% unlike hint cos(90) = 0
% a*.b = (ha* + kb* + lc*).(ua + vb + wc) = ha*.vb = 0
% a* is normal to the plane containing b and c
% a*.b
dot(a_star,b_direct)
% a*.c
dot(a_star,c_direct)
% b* is normal to the plane containing a and c
% b*.a
dot(b_star,a_direct)
% b*.c
dot(b_star,c_direct)
% c* is normal to the plane containing a and b
% c*.a
dot(c_star,a_direct)
% c*.b
dot(c_star,b_direct)
%% Direct lattice metric tensor (G)
%
%  Define direct lattice metric tensor
%
% dot product of direct lattice vectors a,b,c
%
%     |a.a a.b a.c|
% G = |b.a b.b b.c|
%     |c.a c.b c.c|
%
G = ...
[[dot(a_direct,a_direct) dot(a_direct,b_direct) dot(a_direct,c_direct)];...
[ dot(b_direct,a_direct) dot(b_direct,b_direct) dot(b_direct,c_direct)];...
[ dot(c_direct,a_direct) dot(c_direct,b_direct) dot(c_direct,c_direct)]];

G

%% Reciprocal lattice metric tensor (G*)
%
%  Define reciprocal lattice metric tensor
%
% dot product of reciprocal lattice vectors a*,b*,c*
%
%      |a*.a* a*.b* a*.c*|
% G* = |b*.a* b*.b* b*.c*|
%      |c*.a* c*.b* c*.c*|
%
G_star = ...
[[dot(a_star,a_star) dot(a_star,b_star) dot(a_star,c_star)];...
[ dot(b_star,a_star) dot(b_star,b_star) dot(b_star,c_star)];...
[ dot(c_star,a_star) dot(c_star,b_star) dot(c_star,c_star)]];

G_star

%% Applications of the Direct and Reciprocal lattice Metric tensors
%
% 1. magnitude of a direction
% indices of direction as column vector
uvw_100 = [1; 0; 0]
% dot = (uvw_100)t * G * uvw_100
dot_product = uvw_100' * G * uvw_100
% |a| Direct lattice Metric tensor
mag_u100_m = sqrt(dot_product)
% |a| MTEX
u100 = Miller(1,0,0,CS,'An100','uvw')
mag_u100 = norm(u100)
% 2. magnitude of a hkl
% indices of pole to plane as column vector
hkl_100 = [1; 0; 0]
% dot = (hkl_100)t * (G*) * hkl_100
dot_product = hkl_100' * G_star * hkl_100
% |a*| Reciprocal lattice Metric tensor
mag_h100_m = sqrt(dot_product)
% |a*| MTEX
h100 = Miller(1,0,0,CS,'An100','hkl')
mag_h100 = norm(h100)
%%
% 3. angle between directions
% define directions
% indices of direction as column vector
uvw_100 = [1; 0; 0]
% MTEX direction u100 = u.a + v.b + w.c
u100 = Miller(1,0,0,CS,'An100','uvw')
% indices of direction as column vector
uvw_001 = [0; 0; 1]
% MTEX direction u001 = u.a + v.b + w.c
u001 = Miller(0,0,1,CS,'An100','uvw') 
% dot = (uvw_100)t * G * uvw_001
dot_product = uvw_100' * G * uvw_001 
% angle = arccos(dot_product/|u100|/|u001|)
angle_degrees = acosd(dot_product/norm(u100)/norm(u001))
% check using MTEX
angle_between_u100_and_u001 = angle(u100,u001)/degree
%%
% 4. angle between plane normals
% define plane normals
% indices of pole to plane as column vector
hkl_100 = [1; 0; 0]
% MTEX direction h100 = h.a* + k.b* + l.c*
h100 = Miller(1,0,0,CS,'An100','hkl')
% indices of pole to plane as column vector
hkl_001 = [0; 0; 1]
% MTEX direction h001 = h.a* + k.b* + l.c*
h001 = Miller(0,0,1,CS,'An100','hkl') 
% dot = (hkl_100)t * (G*) * hkl_001
dot_product = hkl_100' * G_star * hkl_001 
% angle = arccos(dot_product/|h100|/|h001|))
angle_degrees = acosd(dot_product/norm(h100)/norm(h001))
% check using MTEX
angle_between_h100_and_h001 = angle(h100,h001)/degree

%%
% 5. angle between direction and plane normal
% define plane normals
% indices direction as column vector
uvw_100 = [1; 0; 0]
% MTEX direction uvw_100 = u.a + v.b + w.c
u100 = Miller(1,0,0,CS,'An100','uvw')
% indices of pole to plane as column vector
hkl_001 = [0; 0; 1]
% MTEX direction h001 = h.a* + k.b* + l.c*
h001 = Miller(0,0,1,CS,'An100','hkl') 
% dot = (hkl_100)t * hkl_001
dot_product = uvw_100' * hkl_001 
% angle = arccos(dot_product/|u100|/|h001|))
angle_degrees = acosd(dot_product/norm(u100)/norm(h001))
% check using MTEX
angle_between_u100_and_h001 = angle(u100,h001)/degree
%% 6. Convert hkl to uvw : direction (uvw) normal to plane (hkl)
% define plane normals
% indices of pole to plane as column vector
hkl_100 = [1; 0; 0]
% uvw = (G*) * hkl_100
uvw_h100_m = G_star * hkl_100
% MTEX
h100 = Miller(1,0,0,CS,'An100','hkl')
% h100 -> uvw_u100
uvw_u100 = h100.uvw
%% 7. convert uvw to hkl : normal to plane (hkl) parallel to a direction (uvw)
% define direction
% indices
uvw_100 = [1; 0; 0]
% hkl = G * uvw_100
hkl_u100_m = G * uvw_100
% MTEX
u100 = Miller(1,0,0,CS,'An100','uvw')
% u100  -> hkl_u100
hkl_u100 = u100.hkl
%% Forsterite (101) plane
% define crystal symmetry as orthorhombic
CS = crystalSymmetry('mmm', [4.756 10.207 5.98], 'mineral', 'Forsterite', 'color', 'dark green')
% plane(101)
m = Miller(1,0,1,CS,'Forsterite','hkl')
% plot
plot(m,'labeled')
% save file
saveFigure('/MatLab_Programs/Forsterite_plane_101.pdf')
%% Symmetry equivalent planes to (101)
m_equivalent = symmetrise(m)
% plot
plot(m_equivalent,'labeled')
% save file
saveFigure('/MatLab_Programs/Forsterite_planes_101_equivalents.pdf')
%% Ice Ih direction a1 [2-1-10]
% crystal symmetry
CS = crystalSymmetry('6/mmm', [4.5181 4.5181 7.3560],[90.0 90.0 120.0]*degree,...
    'X||a*','Z||c', 'mineral', 'Ice Ih', 'color', 'light blue')
% direction a1 = [2-1-10]
m = Miller(2,-1,-1,0,CS,'Ice Ih','uvtw')
% plot
plot(m,'labeled')
% save file
saveFigure('/MatLab_Programs/Ice_direction_a1.pdf')
%% Ice Ih symmetry equivalent directions <a>  = <2-1-10>
m_equivalent = symmetrise(m)
% plot
plot(m_equivalent,'labeled')
% save file
saveFigure('/MatLab_Programs/Ice_direction_a1_equivalents.pdf')
%% symmetically equivalent planes and directions
% define crystal symmetry as orthorhombic
CS = crystalSymmetry('mmm', [4.756 10.207 5.98], 'mineral', 'Forsterite', 'color', 'dark green')
% pole to plane hkl h =  h.a* + k.b* + l.c*
h = Miller(1,1,1,CS,'Forsterite','hkl')
% direction u = u.a + v.b + w.c
u = Miller(1,1,1,CS,'Forsterite','uvw')
% angles between h and v
angle_between_h_and_u = angle(h,u)/degree
% symmetically equivalent hkls
h_sym = symmetrise(h)
% symmetically equivalent uvws
u_sym = symmetrise(u)
plot(h_sym,'all','labeled')
hold on
plot(u_sym,'all','labeled')
hold off
%% antipodal (non-polar) option
% symmetically equivalent hkls
h_sym = symmetrise(h,'antipodal')
% symmetically equivalent uvws
u_sym = symmetrise(u,'antipodal')
plot(h_sym,'all','labeled')
hold on
plot(u_sym,'all','labeled')
hold off
%% Define hkl or uvw
% pole to plane hkl h =  h.a* + k.b* + l.c*
h = Miller(1,1,1,CS,'Forsterite','hkl')
% direction u = u.a + v.b + w.c
u = Miller(1,1,1,CS,'Forsterite','uvw')
% angle between h and u
angle_between_h_and_u = angle(h,u)/degree
%% Conversion between hkl and uvw
% N.B. converted uvw and hkl not always integer values
h = Miller(1,1,1,CS,'Forsterite','hkl')
% change hkl (h) to uvw
uvw = h.uvw
u = Miller(1,1,1,CS,'Forsterite','uvw')
% change uvw (u) to hkl
hkl = u.hkl
%% conversion of hkl and uvw to specimen coordinates
% N.B. converted uvw and hkl not always integer values
% change hkl (h) to vector3d r1
% Miller indices h parallel to specimen direction r1
r1 = vector3d(h)
% change uvw (u) to vector3d r2
% Miller indices u parallel to specimen direction r2
r2 = vector3d(u)
%% conversion of specimen coordinates to hkl and uvw
% N.B. converted uvw and hkl not always integer values
% change vector3d(1,1,1) to hkl
hkl_111 = Miller(vector3d(1,1,1),CS,'Forsterite','hkl')
% change vector3d(1,1,1) to uvw
uvw_111 = Miller(vector3d(1,1,1),CS,'Forsterite','uvw')
%% Example of plotting a crystallographic plane
% define crystal symmetry
CS = crystalSymmetry('mmm', [4.756 10.207 5.98], 'mineral', 'Forsterite','color', 'dark green');
% define some planes (hkls)
olivine_hkls = (...
    [Miller(1,0,0,CS,'Forsterite','hkl'),...
     Miller(0,1,0,CS,'Forsterite','hkl'),...
     Miller(0,0,1,CS,'Forsterite','hkl'),...
     Miller(1,1,1,CS,'Forsterite','hkl')])
% plot hkls with their Miller indices
plot(olivine_hkls,'labeled')
hold on
% plane normal to pole to plane (111)
plot(Miller(1,1,1,CS,'Forsterite','hkl'),'plane','color','r','linewidth',3)
% alterative is 
% plane normal to pole to plane (111) with circle with 90*degree option
% circle(Miller(1,1,1,CS,'Forsterite','hkl'),90*degree,'linewidth',3)
% alternative is
% plot(Miller(1,1,1,CS,'Forsterite','hkl'),'plane','color','r','linewidth',3)
hold off
set(gcf,'renderer','zbuffer')
% save plot file as *.png
saveFigure('/MatLab_Programs/Forsterite_plane_plot.png')
%% plot specimen plane and its pole (e.g. 1,1,1)
plot(vector3d(1,1,1),'plane','color','r','linewidth',2)
hold on
plot(vector3d(1,1,1),'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',10)
hold off
%%
angle(Miller(0,0,1,CS,'Forsterite','hkl'),Miller(1,0,3,CS,'Forsterite','hkl'))/degree
angle(Miller(0,0,1,CS,'Forsterite','hkl'),Miller(1,0,2,CS,'Forsterite','uvw'))/degree
angle(vector3d(Miller(-2,0,1,CS,'Forsterite','hkl')),vector3d(Miller(1,0,2,CS,'Forsterite','uvw')))/degree
%% d_hkl inter-planar spacing
d_spacing_100 = dspacing(Miller(1,0,0,CS,'Forsterite','hkl'))
d_spacing_010 = dspacing(Miller(0,1,0,CS,'Forsterite','hkl'))
d_spacing_001 = dspacing(Miller(0,0,1,CS,'Forsterite','hkl'))
% The law of Bravais, Freidel, Donnay and Harker (BFDH)law (1937)
% Donnay J. D. H., Harker D. (1937)
% A new law of crystal morphology extending the law of Bravais.
% Am.Miner. 22, 446-467.
%
% It states that the larger a face's inter-planar spacing,
% the slower its growth-rate and the greater its surface area.
% So d_010>d_001>d_100 then
% (010) should have slowest growth-rate and the greatest surface area.
% (100) should have fastest growth-rate and the smallest surface area.
%% OPX
CS = crystalSymmetry('mmm', [18.2406 8.8302 5.1852], 'mineral','Enstatite')
d_spacing_100 = dspacing(Miller(1,0,0,CS,'Enstatite','hkl'))
d_spacing_010 = dspacing(Miller(0,1,0,CS,'Enstatite','hkl'))
d_spacing_001 = dspacing(Miller(0,0,1,CS,'Enstatite','hkl'))
%% CPX
CS = crystalSymmetry('2/m', [9.746 8.99 5.251], [90,105.63,90]*degree, 'X||a*', 'Y||b*', 'Z||c',...
    'mineral', 'Diopside')
d_spacing_100 = dspacing(Miller(1,0,0,CS,'Diopside','hkl'))
d_spacing_010 = dspacing(Miller(0,1,0,CS,'Diopside','hkl'))
d_spacing_001 = dspacing(Miller(0,0,1,CS,'Diopside','hkl'))
d_spacing_110 = dspacing(Miller(1,1,0,CS,'Diopside','hkl'))
%% AMPHIBOLE
CS = crystalSymmetry('2/m', [9.889 17.978 5.298], [90,105.41,90]*degree, 'X||a*', 'Y||b*', 'Z||c', 'mineral', 'Pargasite')
d_spacing_100 = dspacing(Miller(1,0,0,CS,'Pargasite','hkl'))
d_spacing_010 = dspacing(Miller(0,1,0,CS,'Pargasite','hkl'))
d_spacing_001 = dspacing(Miller(0,0,1,CS,'Pargasite','hkl'))
d_spacing_110 = dspacing(Miller(1,1,0,CS,'Pargasite','hkl'))
%% LAYERED SILICATE
CS= crystalSymmetry('2/m', [5.237 9.07 14.285], [90,97.03,90]*degree, 'X||a*', 'Y||b*', 'Z||c', 'mineral', 'Chlorite')
d_spacing_100 = dspacing(Miller(1,0,0,CS,'Chlorite','hkl'))
d_spacing_010 = dspacing(Miller(0,1,0,CS,'Chlorite','hkl'))
d_spacing_001 = dspacing(Miller(0,0,1,CS,'Chlorite','hkl'))
d_spacing_110 = dspacing(Miller(1,1,0,CS,'Chlorite','hkl'))
%%
CS = crystalSymmetry('4/mmm', [6.5592 6.5592 5.9553], 'mineral', 'Zircon', 'color', 'magenta')
Slip_Direction_111 = Miller(-1,-1,1,CS,'Zircon','uvw')
Slip_Plane_112 = Miller(1,1,2,CS,'Zircon','hkl')
% angles between Slip Direction and Slip Plane normal
angle_between_b_and_n = angle(Slip_Direction_111,Slip_Plane_112)/degree
%%
b = vector3d(Slip_Direction_111)
n = vector3d(Slip_Plane_112)
% angles between Slip Direction (b) and Slip Plane normal (n)
angle_between_b_and_n = angle(b,n)/degree
% rotation axis for pure edge dislocation
rotation_axis =  cross(b,n)
rotation_axis_hkl = Miller(rotation_axis,CS,'Zircon','hkl')
rotation_axis_uvw = Miller(rotation_axis,CS,'Zircon','uvw')
symmetrise(rotation_axis_uvw)
%%
symmetrise(Slip_Direction_111)
%%
symmetrise(Slip_Plane_112)
%% AMPHIBOLE - hornblende
CS = crystalSymmetry('2/m', [9.803 18.046 5.313], [90,105.05,90]*degree, 'X||a*', 'Y||b*', 'Z||c', 'mineral', 'Hornblende')
d_spacing_100 = dspacing(Miller(1,0,0,CS,'Hornblende','hkl'))
d_spacing_010 = dspacing(Miller(0,1,0,CS,'Hornblende','hkl'))
d_spacing_001 = dspacing(Miller(0,0,1,CS,'Hornblende','hkl'))
d_spacing_110 = dspacing(Miller(1,1,0,CS,'Hornblende','hkl'))
% direction u = u.a + v.b + w.c
a_direction = norm(Miller(1,0,0,CS,'Hornblende','uvw'))
b_direction = norm(Miller(0,1,0,CS,'Hornblende','uvw'))
c_direction = norm(Miller(0,0,1,CS,'Hornblende','uvw'))
direction_110 = norm(Miller(1,1,0,CS,'Hornblende','uvw'))
direction_1m10 = norm(Miller(1,-1,0,CS,'Hornblende','uvw'))
direction_m110 = norm(Miller(-1,1,0,CS,'Hornblende','uvw'))
% [001](100) d-space (100) 9.4668  norm [001] = 5.3130
% [010] (100) norm [010] = 18.0460 1/2[010] = 9.0230
% 1/2<110>{1-10} norm<110> = 20.5367  1/2<110> = 10.2683
angle_cl_a_star = angle(Miller(1,1,0,CS,'Hornblende','hkl'),Miller(1,0,0,CS,'Hornblende','hkl'))/degree
angle_cl_b_star = angle(Miller(1,1,0,CS,'Hornblende','hkl'),Miller(0,1,0,CS,'Hornblende','hkl'))/degree
angle_cl_c_star = angle(Miller(1,1,0,CS,'Hornblende','hkl'),Miller(0,0,1,CS,'Hornblende','hkl'))/degree
angle_cl_b = angle(Miller(1,1,0,CS,'Hornblende','hkl'),Miller(0,1,0,CS,'Hornblende','uvw'))/degree
angle_cl_c = angle(Miller(1,1,0,CS,'Hornblende','hkl'),Miller(0,0,1,CS,'Hornblende','uvw'))/degree