/
ExSeismicVelocitySingleCrystalDemo2d.m
280 lines (220 loc) · 8.62 KB
/
ExSeismicVelocitySingleCrystalDemo2d.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
%% Seismic velocities and anisotropy - single crystal plots
%
% 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