1import logging
2import numpy as num
3from pyrocko.guts import Object, Float, List, StringChoice, Int
4from pyrocko.guts_array import Array
6logger = logging.getLogger('pyrocko.gf.tractions')
7km = 1e3
8d2r = num.pi/180.
9r2d = 180./num.pi
12def tukey_window(N, alpha):
13 assert alpha <= 1.
14 window = num.ones(N)
15 n = num.arange(N)
17 N_f = int((alpha * N)//2)
18 window[:N_f] = .5 * (1. - num.cos((2*num.pi * n[:N_f])/(alpha * N)))
19 window[(N-N_f):] = window[:N_f][::-1]
20 return window
23def planck_window(N, epsilon):
24 assert epsilon <= 1.
25 window = num.ones(N)
26 n = num.arange(N)
28 N_f = int((epsilon * N))
29 window[:N_f] = \
30 (1. + num.exp((epsilon * N) / n[:N_f] -
31 (epsilon * N) / ((epsilon * N - n[:N_f]))))**-1.
32 window[(N-N_f):] = window[:N_f][::-1]
33 return window
36class AbstractTractionField(Object):
37 '''
38 Base class for multiplicative traction fields (tapers).
40 Fields of this type a re multiplied in the
41 :py:class:`~pyrocko.gf.tractions.TractionComposition`
42 '''
43 operation = 'mult'
45 def get_tractions(self, nx, ny, patches):
46 raise NotImplementedError
49class TractionField(AbstractTractionField):
50 '''
51 Base class for additive traction fields.
53 Fields of this type are added in the
54 :py:class:`~pyrocko.gf.tractions.TractionComposition`
55 '''
56 operation = 'add'
58 def get_tractions(self, nx, ny, patches):
59 raise NotImplementedError
62class TractionComposition(TractionField):
63 '''
64 Composition of traction fields.
66 :py:class:`~pyrocko.gf.tractions.TractionField` and
67 :py:class:`~pyrocko.gf.tractions.AbstractTractionField` can be combined
68 to realize a combination of different fields.
69 '''
70 components = List.T(
71 AbstractTractionField.T(),
72 default=[],
73 help='Ordered list of tractions.')
75 def get_tractions(self, nx, ny, patches=None):
76 npatches = nx * ny
77 tractions = num.zeros((npatches, 3))
79 for comp in self.components:
80 if comp.operation == 'add':
81 tractions += comp.get_tractions(nx, ny, patches)
82 elif comp.operation == 'mult':
83 tractions *= comp.get_tractions(nx, ny, patches)
84 else:
85 raise AttributeError(
86 'Component %s has an invalid operation %s.' %
87 (comp, comp.operation))
89 return tractions
91 def add_component(self, field):
92 logger.debug('Adding traction component.')
93 self.components.append(field)
96class HomogeneousTractions(TractionField):
97 '''
98 Homogeneous traction field.
100 The traction vectors in strike, dip and normal direction are acting
101 homogeneously on the rupture plane.
102 '''
104 strike = Float.T(
105 default=1.,
106 help='Tractions in strike direction [Pa].')
107 dip = Float.T(
108 default=1.,
109 help='Traction in dip direction (up) [Pa].')
110 normal = Float.T(
111 default=1.,
112 help='Traction in normal direction [Pa].')
114 def get_tractions(self, nx, ny, patches=None):
115 npatches = nx * ny
117 return num.tile(
118 (self.strike, self.dip, self.normal), npatches) \
119 .reshape(-1, 3)
122class DirectedTractions(TractionField):
123 '''
124 Directed traction field.
126 The traction vectors are following a uniform ``rake``.
127 '''
129 rake = Float.T(
130 default=0.,
131 help='Rake angle in [deg], '
132 'measured counter-clockwise from right-horizontal '
133 'in on-plane view. Rake is translated into homogenous tractions '
134 'in strike and up-dip direction.')
135 traction = Float.T(
136 default=1.,
137 help='Traction in rake direction [Pa].')
139 def get_tractions(self, nx, ny, patches=None):
140 npatches = nx * ny
142 strike = num.cos(self.rake*d2r) * self.traction
143 dip = num.sin(self.rake*d2r) * self.traction
144 normal = 0.
146 return num.tile((strike, dip, normal), npatches).reshape(-1, 3)
149class SelfSimilarTractions(TractionField):
150 '''
151 Traction model following Power & Tullis (1991).
153 The traction vectors are calculated as a sum of 2D-cosines with a constant
154 amplitude / wavelength ratio. The wavenumber kx and ky are constant for
155 each cosine function. The rank defines the maximum wavenumber used for
156 summation. So, e.g. a rank of 3 will lead to a summation of cosines with
157 ``kx = ky`` in (1, 2, 3).
158 Each cosine has an associated phases, which defines both the phase shift
159 and also the shift from the rupture plane centre.
160 Finally the summed cosines are translated into shear tractions based on the
161 rake and normalized with ``traction_max``.
163 '''
164 rank = Int.T(
165 default=1,
166 help='Maximum summed cosine wavenumber/spatial frequency.')
168 rake = Float.T(
169 default=0.,
170 help='Rake angle in [deg], '
171 'measured counter-clockwise from right-horizontal '
172 'in on-plane view. Rake is translated into homogenous tractions '
173 'in strike and up-dip direction.')
175 traction_max = Float.T(
176 default=1.,
177 help='Maximum traction vector length [Pa].')
179 phases = Array.T(
180 optional=True,
181 dtype=num.float,
182 shape=(None,),
183 help='Phase shift of the cosines in [rad].')
185 def get_phases(self):
186 if self.phases is not None:
187 if self.phases.shape[0] == self.rank:
188 return self.phases
190 return (num.random.random(self.rank) * 2. - 1.) * num.pi
192 def get_tractions(self, nx, ny, patches=None):
193 z = num.zeros((ny, nx))
194 phases = self.get_phases()
196 for i in range(1, self.rank+1):
197 x = num.linspace(-i*num.pi, i*num.pi, nx) + i*phases[i-1]
198 y = num.linspace(-i*num.pi, i*num.pi, ny) + i*phases[i-1]
199 x, y = num.meshgrid(x, y)
200 r = num.sqrt(x**2 + y**2)
201 z += 1. / i * num.cos(r + phases[i-1])
203 t = num.zeros((nx*ny, 3))
204 t[:, 0] = num.cos(self.rake*d2r) * z.ravel(order='F')
205 t[:, 1] = num.sin(self.rake*d2r) * z.ravel(order='F')
207 t *= self.traction_max / num.max(num.linalg.norm(t, axis=1))
209 return t
212class FractalTractions(TractionField):
213 '''
214 Fractal traction field.
215 '''
217 rseed = Int.T(
218 default=None,
219 optional=True,
220 help='Seed for :py:class:`~numpy.random.RandomState`.'
221 'If ``None``, an random seed will be initialized.')
223 rake = Float.T(
224 default=0.,
225 help='Rake angle in [deg], '
226 'measured counter-clockwise from right-horizontal '
227 'in on-plane view. Rake is translated into homogenous tractions '
228 'in strike and up-dip direction.')
230 traction_max = Float.T(
231 default=1.,
232 help='Maximum traction vector length [Pa].')
234 def __init__(self, *args, **kwargs):
235 super().__init__(*args, **kwargs)
236 if self.rseed is None:
237 self.rseed = num.random.randint(0, 2**32-1)
238 self._data = None
240 def _get_data(self, nx, ny):
241 if self._data is None:
242 rstate = num.random.RandomState(self.rseed)
243 self._data = rstate.rand(nx, ny)
245 return self._data
247 def get_tractions(self, nx, ny, patches=None):
248 if patches is None:
249 raise AttributeError(
250 'Patches needs to be given for this traction field.')
251 npatches = nx * ny
252 dx = -patches[0].al1 + patches[0].al2
253 dy = -patches[0].aw1 + patches[0].aw2
255 # Create random data and get spectrum and power spectrum
256 data = self._get_data(nx, ny)
257 spec = num.fft.fftshift(num.fft.fft2(data))
258 power_spec = (num.abs(spec)/spec.size)**2
260 # Get 0-centered wavenumbers (k_rad == 0.) is in the centre
261 kx = num.fft.fftshift(num.fft.fftfreq(nx, d=dx))
262 ky = num.fft.fftshift(num.fft.fftfreq(ny, d=dy))
263 k_rad = num.sqrt(ky[:, num.newaxis]**2 + kx[num.newaxis, :]**2)
265 # Define wavenumber bins
266 k_bins = num.arange(0, num.max(k_rad), num.max(k_rad)/10.)
268 # Set amplitudes within wavenumber bins to power_spec * 1 / k_max
269 amps = num.zeros(k_rad.shape)
270 amps[k_rad == 0.] = 1.
272 for i in range(k_bins.size-1):
273 k_min = k_bins[i]
274 k_max = k_bins[i+1]
275 r = num.logical_and(k_rad > k_min, k_rad <= k_max)
276 amps[r] = power_spec.T[r]
277 amps = num.sqrt(amps * data.size * num.pi * 4)
279 amps[k_rad > k_bins.max()] = power_spec.ravel()[num.argmax(power_spec)]
281 # Multiply spectrum by amplitudes and inverse fft into demeaned noise
282 spec *= amps.T
284 tractions = num.abs(num.fft.ifft2(spec))
285 tractions -= num.mean(tractions)
286 tractions *= self.traction_max / num.abs(tractions).max()
288 t = num.zeros((npatches, 3))
289 t[:, 0] = num.cos(self.rake*d2r) * tractions.ravel(order='C')
290 t[:, 1] = num.sin(self.rake*d2r) * tractions.ravel(order='C')
292 return t
295class RectangularTaper(AbstractTractionField):
296 width = Float.T(
297 default=.2,
298 help='Width of the taper as a fraction of the plane.')
300 type = StringChoice.T(
301 choices=('tukey', ),
302 default='tukey',
303 help='Type of the taper, default: "tukey".')
305 def get_tractions(self, nx, ny, patches=None):
306 if self.type == 'tukey':
307 x = tukey_window(nx, self.width)
308 y = tukey_window(ny, self.width)
309 return (x[:, num.newaxis] * y).ravel()[:, num.newaxis]
311 raise AttributeError('Unknown type: %s' % self.type)
314class DepthTaper(AbstractTractionField):
315 depth_start = Float.T(
316 help='Depth where the taper begins [m].')
318 depth_stop = Float.T(
319 help='Depth where taper ends and drops to zero [m].')
321 type = StringChoice.T(
322 choices=('linear', ),
323 default='linear',
324 help='Type of the taper, default: "linear".')
326 def get_tractions(self, nx, ny, patches):
327 assert self.depth_stop > self.depth_start
328 depths = num.array([p.depth for p in patches])
330 if self.type == 'linear':
331 slope = self.depth_stop - self.depth_start
332 depths -= self.depth_stop
333 depths /= -slope
334 depths[depths > 1.] = 1.
335 depths[depths < 0.] = 0.
336 return depths[:, num.newaxis]
339def plot_tractions(tractions, nx=15, ny=12, depth=10*km, component='strike'):
340 '''
341 Plot choosen traction model for quick inspection.
343 :param tractions: traction field or traction composition to be displayed
344 :type tractions: :py:class:`pyrocko.gf.tractions.TractionField`
345 :param nx: number of patches along strike
346 :type nx: optional, int
347 :param ny: number of patches down dip
348 :type ny: optional, int
349 :param depth: depth of the rupture plane center in [m]
350 :type depth: optional, float
351 :param component: choice, which component of the traction is displayed.
352 Choose between:
353 - "tx": along strike tractions
354 - "ty": up dip tractions
355 - "tz": normal tractions
356 - "absolut": length of the traction vector
357 :type component: optional, str
358 '''
359 import matplotlib.pyplot as plt
360 from pyrocko.modelling.okada import OkadaSource
362 comp2idx = dict(
363 tx=0, ty=1, tz=2)
365 source = OkadaSource(
366 lat=0.,
367 lon=0.,
368 depth=depth,
369 al1=-20*km, al2=20*km,
370 aw1=-15*km, aw2=15*km,
371 strike=120., dip=90., rake=90.,
372 slip=5.)
374 patches, _ = source.discretize(nx, ny)
375 tractions = tractions.get_tractions(nx, ny, patches)
377 if component in comp2idx:
378 tractions = tractions[:, comp2idx[component]].reshape(nx, ny)
379 elif component == 'absolut':
380 tractions = num.linalg.norm(tractions, axis=1).reshape(nx, ny)
381 else:
382 raise ValueError('Given component is not valid.')
384 fig = plt.figure()
385 ax = fig.gca()
387 ax.imshow(tractions)
389 plt.show()
392__all__ = [
393 'AbstractTractionField',
394 'TractionField',
395 'TractionComposition',
396 'HomogeneousTractions',
397 'DirectedTractions',
398 'FractalTractions',
399 'SelfSimilarTractions',
400 'RectangularTaper',
401 'DepthTaper',
402 'plot_tractions']