% Calculating and plotting elastic velocities (km/s) from elasticity stiffness
%
% Illustrates 3 types of plot
%
% A) Func_Seismic_multiplot_MTEX4 with Vp,AVs,Vs1 polarization, Vs1,Vs2,dVs,Vs1/Vp,Vs2/Vp 
%
% B) Func_Seismic_multiplot_Vp_AVs_MTEX4 with Vp,AVS with polarization options
%
% C) Func_Seismic_multiplot_Vp_dVs_MTEX4 with Vp,dVS with polarization options
%
% Cijkl tensor (GPa) and density (g/cm3)
%
% Tested with MTEX 4.0.9
%
% By David Mainprice (18/12/2014)
%****************************************
% Requires
% Func_Seismic_multiplot_MTEX4
% Func_Seismic_multiplot_Vp_AVs_MTEX4
% Func_Seismic_multiplot_Vp_dVs_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_MTEX4(C,rho,plot_name,plot_x,Vp_polar,Vs1_polar,Vs2_polar,ref_axis_1,ref_axis_2,ref_axis_3)
%
% 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.png';
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]';
%
Func_Seismic_multiplot_MTEX4(C,rho,plot_name,plot_x,Vp_polar,Vs1_polar,Vs2_polar,ref_axis_1,ref_axis_2,ref_axis_3)
%
%% Func_Seismic_multiplot_Vp_AVs_MTEX4(C,rho,plot_name,plot_x,Vp_polar,Vs1_polar,Vs2_polar,ref_axis_1,ref_axis_2,ref_axis_3)
% Seismic_multiplot_Vp_AVs
%
% 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_seismics_Vp_AVs_plot_x_north.png';
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]';
%
Func_Seismic_multiplot_Vp_AVs_MTEX4(C,rho,plot_name,plot_x,Vp_polar,Vs1_polar,Vs2_polar,ref_axis_1,ref_axis_2,ref_axis_3)

%% Func_Seismic_multiplot_Vp_dVs_MTEX4(C,rho,plot_name,plot_x,Vp_polar,Vs1_polar,Vs2_polar,ref_axis_1,ref_axis_2,ref_axis_3)
% Seismic_multiplot_Vp_AVs
%
% 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_seismics_Vp_dVs_plot_x_north.png';
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]';
%
Func_Seismic_multiplot_Vp_dVs_MTEX4(C,rho,plot_name,plot_x,Vp_polar,Vs1_polar,Vs2_polar,ref_axis_1,ref_axis_2,ref_axis_3)
%
disp(' ')
disp(' Script has successfully terminated')
disp(' ')
%
% End of script
%
