Plot 2D velocity surfaces - phase, slowness & wavefront For given plane normal direction, for example xvector, yvector or zvector Stishovite (SiO2) at high pressure
David Mainprice 1/05/2018
The Elastic Stiffness Tensor of Stishovite Reference: Elastic constants Stishovite, Weidner et al 1982 JGR
Reference: Crystal Structure Sinclair and Ringwood 78 Nature P422
d=4.2901
% Define density (g/cm3)
rho = 4.2901 ;
% elastic Cij stiffness tensor (GPa) as matrix M
M = ...
[[ 453.00 211.00 203.00 0.00 0.00 0.00 ]; ...
[ 211.00 453.00 203.00 0.00 0.00 0.00 ]; ...
[ 203.00 203.00 776.00 0.00 0.00 0.00 ]; ...
[ 0.00 0.00 0.00 252.00 0.00 0.00 ]; ...
[ 0.00 0.00 0.00 0.00 252.00 0.00 ]; ...
[ 0.00 0.00 0.00 0.00 0.00 302.00 ]];
% define cartesian tensor crystal crystalSymmetry & frame
cs_Tensor = crystalSymmetry ( '4/mmm' ,[ 4.9133 4.9133 5.4048 ], ...
[ 90.0000 90.0000 90.0000 ] * degree , 'x||a' , 'z||c' , ...
'mineral' , 'Stishovite 1982' );
%
% define elastic stiffness tensor Cijkl
C = stiffnessTensor ( M , cs_Tensor , 'density' , rho );
compute wave velocities and polarization directions
% the propagation direction is just the vector normal to the sphere
prop = S2AxisFieldHarmonic . normal ;
% the wave velocities and polarization directions as directional dependend
% functions
[ vp , vs1 , vs2 , pp , ps1 , ps2 ] = velocity ( C );
Plotting settings
% plane normal direction for 2d sections
planeNormal = vector3d . Z ;
% plotting convention - plot X-axis to east
plotx2east ;
% close all open graphics
close all
% some global options for the titles
%titleOpt = {'FontSize',getMTEXpref('FontSize'),'visible','on','color','k'};
titleOpt = { 'FontSize' , 25 };
% option for legend
legendOpt = { 'location' , 'best' };
% option for mtexColorbar
ColorbarOpt = { 'location' , 'southoutside' };
% options for sections
optSec = { 'color' , 'interp' , 'linewidth' , 5 , 'doNotDraw' };
% options for quiver
optQuiver = { 'linewidth' , 5 , 'autoScaleFactor' , 0.25 , 'doNotDraw' };
% options for prop
optQuiverProp = { 'color' , 'k' , 'linewidth' , 2 , 'autoScaleFactor' , 0.15 , 'doNotDraw' };
1: Phase velocity surface (km/s)
figure ( 1 )
% phase velocities
plotSection ( vp , planeNormal , optSec {:})
hold on
plotSection ( vs1 , planeNormal , optSec {:})
plotSection ( vs2 , planeNormal , optSec {:})
% polarization directions
quiverSection ( vp , pp , planeNormal , 'color' , 'c' , optQuiver {:})
quiverSection ( vs1 , ps1 , planeNormal , 'color' , 'g' , optQuiver {:})
quiverSection ( vs2 , ps2 , planeNormal , 'color' , 'm' , optQuiver {:})
% plot propagation directions as reference
quiverSection ( vp , prop , planeNormal , optQuiverProp {:})
quiverSection ( vs1 , prop , planeNormal , optQuiverProp {:})
quiverSection ( vs2 , prop , planeNormal , optQuiverProp {:})
axis off tight
legend ( 'Vp' , 'Vs1' , 'Vs2' , 'pp' , 'ps1' , 'ps2' , 'X' , legendOpt {:})
mtexTitle ( 'Phase velocity surfaces (km/s)' , titleOpt {:})
% seismic velocity slow = red 2 blue =fast
mtexColorMap red2blue
mtexColorbar ( 'Title' , '(km/s)' , ColorbarOpt {:})
% Az El for planeNormal = Z
%view(0,85)
hold off
drawNow ( gcm , 'figSize' , 'large' )
2: Slowness surface (s/km)
figure ( 2 )
% slowness
plotSection ( 1. / vp , planeNormal , optSec {:})
hold on
plotSection ( 1. / vs1 , planeNormal , optSec {:})
plotSection ( 1. / vs2 , planeNormal , optSec {:})
% polarization directions
quiverSection ( 1. / vp , pp , planeNormal , 'color' , 'c' , optQuiver {:})
quiverSection ( 1. / vs1 , ps1 , planeNormal , 'color' , 'g' , optQuiver {:})
quiverSection ( 1. / vs2 , ps2 , planeNormal , 'color' , 'm' , optQuiver {:})
% plot propagation directions as reference
quiverSection ( 1. / vp , prop , planeNormal , optQuiverProp {:})
quiverSection ( 1. / vs1 , prop , planeNormal , optQuiverProp {:})
quiverSection ( 1. / vs2 , prop , planeNormal , optQuiverProp {:})
axis off tight
legend ( 'Vp' , 'Vs1' , 'Vs2' , 'pp' , 'ps1' , 'ps2' , 'X' , legendOpt {:})
mtexTitle ( 'Slowness surfaces (km/s)' , titleOpt {:})
% seismic slowness slow = blue 2 red =fast
mtexColorMap blue2red
mtexColorbar ( 'Title' , '(s/km)' , ColorbarOpt {:})
% Az El for planeNormal = Z
%view(0,85)
hold off
drawNow ( gcm , 'figSize' , 'large' )
Select the two S-waves (Vs1 and Vs2 where Vs1>Vs2 in velocity) by the orientation of the polarization vectors ps1 and ps2 with respect to the specimen Z direction. So that sv (v=vertical) is S-wave with polarization closest to Z sh (h=horizontal) has the polarization closest to the plane normal to Z Both polarizations pairs (sv and sh) and (ps1 and ps2) are orthogonal
which values to switch this defines a function which is either one or zero
id = angle ( ps1 , vector3d . Z ) <= 89.9 * degree ;
vsv = id .* vs1 + ( 1 - id ) .* vs2 ;
vsh = id .* vs2 + ( 1 - id ) .* vs1 ;
psv = id .* ps1 + ( 1 - id ) .* ps2 ;
psh = id .* ps2 + ( 1 - id ) .* ps1 ;
1: Phase velocity surface (km/s) with sv1 & vs2
figure ( 1 )
% phase velocities
plotSection ( vp , planeNormal , optSec {:})
hold on
plotSection ( vs1 , planeNormal , optSec {:})
plotSection ( vs2 , planeNormal , optSec {:})
% polarization directions
quiverSection ( vp , pp , planeNormal , 'color' , 'c' , optQuiver {:})
quiverSection ( vs1 , ps1 , planeNormal , 'color' , 'g' , optQuiver {:})
quiverSection ( vs2 , ps2 , planeNormal , 'color' , 'm' , optQuiver {:})
% plot propagation directions as reference
quiverSection ( vp , prop , planeNormal , optQuiverProp {:})
quiverSection ( vs1 , prop , planeNormal , optQuiverProp {:})
quiverSection ( vs2 , prop , planeNormal , optQuiverProp {:})
axis off tight
legend ( 'Vp' , 'Vs1' , 'Vs2' , 'pp' , 'ps1' , 'ps2' , 'X' , legendOpt {:})
mtexTitle ( 'Phase velocity surfaces (km/s)' , titleOpt {:})
mtexColorMap blue2red
mtexColorbar ( 'Title' , '(km/s)' , ColorbarOpt {:})
% Az El for planeNormal = Z
%view(0,85)
hold off
drawNow ( gcm , 'figSize' , 'large' )
2: Phase velocity surface (km/s) with svs & vsh
figure ( 2 )
% phase velocities
plotSection ( vp , planeNormal , optSec {:})
hold on
plotSection ( vsv , planeNormal , optSec {:})
plotSection ( vsh , planeNormal , optSec {:})
% polarization directions
quiverSection ( vp , pp , planeNormal , 'color' , 'c' , optQuiver {:})
quiverSection ( vsv , ps1 , planeNormal , 'color' , 'g' , optQuiver {:})
quiverSection ( vsh , ps2 , planeNormal , 'color' , 'm' , optQuiver {:})
% plot propagation directions as reference
quiverSection ( vp , prop , planeNormal , optQuiverProp {:})
quiverSection ( vsv , prop , planeNormal , optQuiverProp {:})
quiverSection ( vsh , prop , planeNormal , optQuiverProp {:})
axis off tight
legend ( 'Vp' , 'Vsv' , 'Vsh' , 'pp' , 'ps1' , 'ps2' , 'X' , legendOpt {:})
mtexTitle ( 'Phase velocity surfaces (km/s)' , titleOpt {:})
mtexColorMap blue2red
mtexColorbar ( 'Title' , '(km/s)' , ColorbarOpt {:})
% Az El for planeNormal = Z
%view(0,85)
hold off
drawNow ( gcm , 'figSize' , 'large' )
3: plot slowness in plane normal Z
figure ( 3 )
plotSection ( 1. / vp , planeNormal , optSec {:})
hold on
plotSection ( 1. / vs1 , planeNormal , optSec {:})
plotSection ( 1. / vs2 , planeNormal , optSec {:})
% polarization vectors pp,ps1,ps2
quiverSection ( 1. / vp , pp , planeNormal , 'color' , 'c' , optQuiver {:})
quiverSection ( 1. / vs1 , ps1 , planeNormal , 'color' , 'g' , optQuiver {:})
quiverSection ( 1. / vs2 , ps2 , planeNormal , 'color' , 'm' , optQuiver {:})
% propagation vectors (prop)
quiverSection ( 1. / vp , prop , planeNormal , optQuiverProp {:})
quiverSection ( 1. / vs1 , prop , planeNormal , optQuiverProp {:})
quiverSection ( 1. / vs2 , prop , planeNormal , optQuiverProp {:})
axis off tight
legend ( 'Sp' , 'Ss1' , 'Ss2' , 'Sp' , 'Ss1' , 'Ss2' , 'X' , legendOpt {:})
mtexTitle ( 'Stishovite : Slowness surfaces (s/km)' , titleOpt {:})
% seismic slowness slow = blue 2 red =fast
mtexColorMap red2blue
mtexColorbar ( 'Title' , '(s/km)' , ColorbarOpt {:})
hold off
drawNow ( gcm , 'figSize' , 'large' )
Compute WaveFront as spherical functions: EnergyVectors Evp, Evs1 & Evs2
Evp = energyVector ( C ,[], vp , pp );
Evs1 = energyVector ( C ,[], vs1 , ps1 );
Evs2 = energyVector ( C ,[], vs2 , ps2 );
plot wavefront in plane normal Z
close all
% Vp,Vs1,Vs2 wavefronts with motif
%optSec = {'linewidth',5,'doNotDraw'};
optSec = { 'color' , 'interp' , 'linewidth' , 2 , 'doNotDraw' };
plotSection ( Evp , planeNormal , optSec {:})
hold on
plotSection ( Evs1 , planeNormal , optSec {:})
plotSection ( Evs2 , planeNormal , optSec {:})
%
% polarization vectors pp,ps1,ps2
quiverSection ( Evp , Evp , planeNormal , 'color' , 'c' , optQuiver {:})
quiverSection ( Evs1 , Evs1 , planeNormal , 'color' , 'g' , optQuiver {:})
quiverSection ( Evs2 , Evs2 , planeNormal , 'color' , 'm' , optQuiver {:})
% propagation vectors (prop)
quiverSection ( Evp , prop , planeNormal , optQuiverProp {:})
quiverSection ( Evs1 , prop , planeNormal , optQuiverProp {:})
quiverSection ( Evs2 , prop , planeNormal , optQuiverProp {:})
axis off tight
legend ( 'Evp' , 'Es1' , 'Es2' , 'Epv' , 'Eps1' , 'Eps2' , 'X' , legendOpt {:})
mtexTitle ( 'Stishovite : Wavefront surfaces (km/s)' , titleOpt {:})
mtexColorbar ( 'Title' , '(km/s)' , ColorbarOpt {:})
hold off
drawNow ( gcm , 'figSize' , 'large' )
% End of demo