% Calculating and plotting elastic velocities (km/s) from elasticity stiffness
%
% Illustrates types of plot
%
% Cijkl tensor (GPa) and density (g/cm3)
%
% Tested with MTEX 4.0.9
%
% By David Mainprice (9/12/2014)
%****************************************
% Requires
% Func_Seismic_multiplot_gfs_MTEX4
% In Current Folder
%****************************************
%
%% clear memory and close all plots
clear all
close all
%% Crystal Symmetry
% 
% crystal symmetry - Orthorhombic mmm
% Olivine structure 
% (4.7646 10.2296 5.9942 90.00 90.00 90.00) - Orthorhombic

CS_tensor = crystalSymmetry('mmm',[4.7646,10.2296,5.9942],...
    [90.00,90.00,90.00]*degree,'x||a','z||c','mineral','Olivine');

%% Import 4th rank tensor as 6 by 6 matrix
% 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 (g/cm3)
rho=3.355;

% Define tenor object in MTEX
% Cij -> Cijkl - elastic stiffness tensor
C = tensor(M,CS_tensor,'name','elastic stiffness','unit','GPa','density',rho)
%% Func_Seismic_multiplot_gfs_MTEX4(C,rho,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)
%
% C - MTEX object elastic stiffness tensor
% rho - density (gm/cm3)
% plot_name - plot name (can include path '/MatLab_Programs/Test_plot.pdf')
% plot_x : 1= X to north, 2= X to East, 3= X to South, 4= X to West
% ref_axis_1,ref_axis_2,ref_axis_3 = text reference axes labels e.g.;
%                                    'X','Y','Z' or '[100]','[010]','[001]'
plot_name = '/MatLab_Programs/Test_olivine_single_crystal_complete_seismics_plot_x_north.pdf';
plot_x = 1;     % plot X to North
Vp_polar = 0;   % plot Vp polarizations
Vs1_polar = 1;  % plot Vs1 polarizations on AVs pole figure
Vs2_polar = 0;  % plot Vs2 polarizations on Avs pole figure
% reference axes - text labels
ref_axis_1 ='[100]';
ref_axis_2 ='[010]';
ref_axis_3 ='[001]';
% define directions for text labels
ref_vec_1 = xvector;
ref_vec_2 = yvector;
ref_vec_3 = zvector;
% plot hemisphere 'upper' or 'lower'
HS ='upper';
% global_fontSize linked to plot size
global_fontSize = 20;
% Seismic multiplot
Func_Seismic_multiplot_gfs_MTEX4(C,rho,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)
%% Import 4th rank tensor as 6 by 6 matrix
% Reference Elastic constants
% Alexsandrov et al(1974) An57-60 XIIa YII(Z x X) ZIIc*
% Elastic constants measured as monoclinic called "pseudo-monoclinic"
% 
% Reference Crystal Structure
% LABRADORITE  An 66 Wenk and Kroll(1985)
%
% N.B. Crystal structure is triclinic
% Define Cartesian Tensor crystal symmetry & frame
CS_tensor = crystalSymmetry('-1',[  8.1747 12.8708 14.2030],...
 [  93.4610 116.0860  90.5140]*degree,'x||a','z||c*',...
 'mineral','An57-60 1974')
%
% elastic Cij stiffness tensor (GPa) as matrix M
 M =.... 
 [[   99.40   62.80   48.70    0.00   -0.25    0.00];...
 [    62.80  158.40   26.70    0.00  -10.70    0.00];...
 [    48.70   26.70  150.30    0.00  -12.40    0.00];...
 [     0.00    0.00    0.00   21.70    0.00   -5.40];...
 [    -0.25  -10.70  -12.40    0.00   34.50    0.00];...
 [     0.00    0.00    0.00   -5.40    0.00   37.10]];
%
% Define density (g/cm3)
rho= 2.68;
% Define tenor object in MTEX
% Cij -> Cijkl - elastic stiffness tensor
C = tensor(M,CS_tensor,'name','elastic stiffness','unit','GPa','density',rho)
%% Func_Seismic_multiplot_gfs_MTEX4(C,rho,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)
%
% C - MTEX object elastic stiffness tensor
% rho - density (gm/cm3)
% plot_name - plot name (can include path '/MatLab_Programs/Test_plot.pdf')
% plot_x : 1= X to north, 2= X to East, 3= X to South, 4= X to West
% ref_axis_1,ref_axis_2,ref_axis_3 = text reference axes labels e.g.;
%                                    'X','Y','Z' or '[100]','[010]','[001]'
plot_name = '/MatLab_Programs/Test_An57_single_crystal_complete_seismics_plot_x_north.pdf';
plot_x = 1;     % plot X to North
Vp_polar = 0;   % plot Vp polarizations
Vs1_polar = 1;  % plot Vs1 polarizations on AVs pole figure
Vs2_polar = 0;  % plot Vs2 polarizations on Avs pole figure
% reference axes - text labels
ref_axis_1 ='a[100]';
ref_axis_2 ='b[010]';
ref_axis_3 ='c*(001)';
% define directions for text labels
ref_vec_1 = xvector;
ref_vec_2 = yvector;
ref_vec_3 = zvector;
% plot hemisphere 'upper' or 'lower'
HS ='upper';
% global_fontSize linked to plot size
global_fontSize = 20;
% Seismic multiplot
Func_Seismic_multiplot_gfs_MTEX4(C,rho,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)
%%
disp(' ')
disp(' Script has successfully terminated')
disp(' ')
%
% End of script
%
