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 ''' Abstract traction field
39 Fields of this type a re multiplied in the
40 :py:class:`~pyrocko.gf.tractions.TractionComposition`
41 '''
42 operation = 'mult'
44 def get_tractions(self, nx, ny, patches):
45 raise NotImplementedError
48class TractionField(AbstractTractionField):
49 ''' Traction field
51 Fields of this type are added in the
52 :py:class:`~pyrocko.gf.tractions.TractionComposition`
53 '''
54 operation = 'add'
56 def get_tractions(self, nx, ny, patches):
57 raise NotImplementedError
60class TractionComposition(TractionField):
61 ''' Composition of traction fields
63 :py:class:`~pyrocko.gf.tractions.TractionField` and
64 :py:class:`~pyrocko.gf.tractions.AbstractTractionField` can be combined
65 to realize a combination of different fields.
66 '''
67 components = List.T(
68 AbstractTractionField.T(),
69 default=[],
70 help='Ordered list of tractions')
72 def get_tractions(self, nx, ny, patches=None):
73 npatches = nx * ny
74 tractions = num.zeros((npatches, 3))
76 for comp in self.components:
77 if comp.operation == 'add':
78 tractions += comp.get_tractions(nx, ny, patches)
79 elif comp.operation == 'mult':
80 tractions *= comp.get_tractions(nx, ny, patches)
81 else:
82 raise AttributeError(
83 'Component %s has an invalid operation %s.' %
84 (comp, comp.operation))
86 return tractions
88 def add_component(self, field):
89 logger.debug('adding traction component')
90 self.components.append(field)
93class UniformTractions(TractionField):
94 ''' Uniform traction field
96 The traction field is uniform in strike, dip and normal direction.
97 This realisation is not only simple but also unrealistic.
98 '''
99 traction = Float.T(
100 default=1.,
101 help='Uniform traction in strike, dip and normal direction [Pa]')
103 def get_tractions(self, nx, ny, patches=None):
104 npatches = nx * ny
105 return num.full((npatches, 3), self.traction)
108class HomogeneousTractions(TractionField):
109 ''' Homogeneous traction field
111 The traction vectors in strike, dip and normal direction are acting
112 homogeneously on the rupture plane.
113 '''
115 strike = Float.T(
116 default=1.,
117 help='Tractions in strike direction [Pa]')
118 dip = Float.T(
119 default=1.,
120 help='Traction in dip direction (up) [Pa]')
121 normal = Float.T(
122 default=1.,
123 help='Traction in normal direction [Pa]')
125 def get_tractions(self, nx, ny, patches=None):
126 npatches = nx * ny
128 return num.tile(
129 (self.strike, self.dip, self.normal), npatches) \
130 .reshape(-1, 3)
133class DirectedTractions(TractionField):
134 ''' Directed traction field
136 The traction vectors are following a uniform ``rake``.
137 '''
139 rake = Float.T(
140 default=0.,
141 help='rake angle in [deg], '
142 'measured counter-clockwise from right-horizontal '
143 'in on-plane view. Rake is translated into homogenous tractions '
144 'in strike and up-dip direction.')
145 traction = Float.T(
146 default=1.,
147 help='Traction in rake direction [Pa]')
149 def get_tractions(self, nx, ny, patches=None):
150 npatches = nx * ny
152 strike = num.cos(self.rake*d2r) * self.traction
153 dip = num.sin(self.rake*d2r) * self.traction
154 normal = 0.
156 return num.tile((strike, dip, normal), npatches).reshape(-1, 3)
159class SelfSimilarTractions(TractionField):
160 ''' Traction model following Power & Tullis (1991).
162 The traction vectors are calculated as a sum of 2D-cosines with a constant
163 amplitude / wavelength ratio. The wavenumber kx and ky are constant for
164 each cosine function. The rank defines the maximum wavenumber used for
165 summation. So, e.g. a rank of 3 will lead to a summation of cosines with
166 ``kx = ky`` in (1, 2, 3).
167 Each cosine has an associated phases, which defines both the phase shift
168 and also the shift from the rupture plane centre.
169 Finally the summed cosines are translated into shear tractions based on the
170 rake and normalized with ``traction_max``.
172 '''
173 rank = Int.T(
174 default=1,
175 help='maximum summed cosine wavenumber/spatial frequency.')
177 rake = Float.T(
178 default=0.,
179 help='rake angle in [deg], '
180 'measured counter-clockwise from right-horizontal '
181 'in on-plane view. Rake is translated into homogenous tractions '
182 'in strike and up-dip direction.')
184 traction_max = Float.T(
185 default=1.,
186 help='maximum traction vector length [Pa]')
188 phases = Array.T(
189 optional=True,
190 dtype=num.float,
191 shape=(None,),
192 help='phase shift of the cosines in [rad].')
194 def get_phases(self):
195 if self.phases is not None:
196 if self.phases.shape[0] == self.rank:
197 return self.phases
199 return (num.random.random(self.rank) * 2. - 1.) * num.pi
201 def get_tractions(self, nx, ny, patches=None):
202 z = num.zeros((ny, nx))
203 phases = self.get_phases()
205 for i in range(1, self.rank+1):
206 x = num.linspace(-i*num.pi, i*num.pi, nx) + i*phases[i-1]
207 y = num.linspace(-i*num.pi, i*num.pi, ny) + i*phases[i-1]
208 x, y = num.meshgrid(x, y)
209 r = num.sqrt(x**2 + y**2)
210 z += 1. / i * num.cos(r + phases[i-1])
212 t = num.zeros((nx*ny, 3))
213 t[:, 0] = num.cos(self.rake*d2r) * z.ravel(order='F')
214 t[:, 1] = num.sin(self.rake*d2r) * z.ravel(order='F')
216 t *= self.traction_max / num.max(num.linalg.norm(t, axis=1))
218 return t
221class FractalTractions(TractionField):
222 ''' Fractal traction field
223 '''
225 rseed = Int.T(
226 default=None,
227 optional=True,
228 help='Seed for :py:class:`~numpy.random.RandomState`.'
229 'If ``None``, an random seed will be initialized.')
231 rake = Float.T(
232 default=0.,
233 help='rake angle in [deg], '
234 'measured counter-clockwise from right-horizontal '
235 'in on-plane view. Rake is translated into homogenous tractions '
236 'in strike and up-dip direction.')
238 traction_max = Float.T(
239 default=1.,
240 help='maximum traction vector length [Pa]')
242 def __init__(self, *args, **kwargs):
243 super().__init__(*args, **kwargs)
244 if self.rseed is None:
245 self.rseed = num.random.randint(0, 2**32-1)
246 self._data = None
248 def _get_data(self, nx, ny):
249 if self._data is None:
250 rstate = num.random.RandomState(self.rseed)
251 self._data = rstate.rand(nx, ny)
253 return self._data
255 def get_tractions(self, nx, ny, patches=None):
256 if patches is None:
257 raise AttributeError(
258 'patches needs to be given for this traction field')
259 npatches = nx * ny
260 dx = -patches[0].al1 + patches[0].al2
261 dy = -patches[0].aw1 + patches[0].aw2
263 # Create random data and get spectrum and power spectrum
264 data = self._get_data(nx, ny)
265 spec = num.fft.fftshift(num.fft.fft2(data))
266 power_spec = (num.abs(spec)/spec.size)**2
268 # Get 0-centered wavenumbers (k_rad == 0.) is in the centre
269 kx = num.fft.fftshift(num.fft.fftfreq(nx, d=dx))
270 ky = num.fft.fftshift(num.fft.fftfreq(ny, d=dy))
271 k_rad = num.sqrt(ky[:, num.newaxis]**2 + kx[num.newaxis, :]**2)
273 # Define wavenumber bins
274 k_bins = num.arange(0, num.max(k_rad), num.max(k_rad)/10.)
276 # Set amplitudes within wavenumber bins to power_spec * 1 / k_max
277 amps = num.zeros(k_rad.shape)
278 amps[k_rad == 0.] = 1.
280 for i in range(k_bins.size-1):
281 k_min = k_bins[i]
282 k_max = k_bins[i+1]
283 r = num.logical_and(k_rad > k_min, k_rad <= k_max)
284 amps[r] = power_spec.T[r]
285 amps = num.sqrt(amps * data.size * num.pi * 4)
287 amps[k_rad > k_bins.max()] = power_spec.ravel()[num.argmax(power_spec)]
289 # Multiply spectrum by amplitudes and inverse fft into demeaned noise
290 spec *= amps.T
292 tractions = num.abs(num.fft.ifft2(spec))
293 tractions -= num.mean(tractions)
294 tractions *= self.traction_max / num.abs(tractions).max()
296 t = num.zeros((npatches, 3))
297 t[:, 0] = num.cos(self.rake*d2r) * tractions.ravel(order='C')
298 t[:, 1] = num.sin(self.rake*d2r) * tractions.ravel(order='C')
300 return t
303class RectangularTaper(AbstractTractionField):
304 width = Float.T(
305 default=.2,
306 help='Width of the taper as a fraction of the plane.')
308 type = StringChoice.T(
309 choices=('tukey', ),
310 default='tukey',
311 help='Type of the taper, default "tukey"')
313 def get_tractions(self, nx, ny, patches=None):
314 if self.type == 'tukey':
315 x = tukey_window(nx, self.width)
316 y = tukey_window(ny, self.width)
317 return (x[:, num.newaxis] * y).ravel()[:, num.newaxis]
319 raise AttributeError('unknown type %s' % self.type)
322class DepthTaper(AbstractTractionField):
323 depth_start = Float.T(
324 help='Depth where the taper begins [km]')
326 depth_stop = Float.T(
327 help='Depth where taper stops, and drops to 0. [km]')
329 type = StringChoice.T(
330 choices=('linear', ),
331 default='linear',
332 help='Type of the taper, default "linear"')
334 def get_tractions(self, nx, ny, patches):
335 assert self.depth_stop > self.depth_start
336 depths = num.array([p.depth for p in patches])
338 if self.type == 'linear':
339 slope = self.depth_stop - self.depth_start
340 depths -= self.depth_stop
341 depths /= -slope
342 depths[depths > 1.] = 1.
343 depths[depths < 0.] = 0.
344 return depths[:, num.newaxis]
347def plot_tractions(tractions, nx=15, ny=12, depth=10*km, component='strike'):
348 '''Plot choosen traction model for quick inspection
350 :param tractions: traction field or traction composition to be displayed
351 :type tractions: :py:class:`pyrocko.gf.tractions.TractionField`
352 :param nx: number of patches along strike
353 :type nx: optional, int
354 :param ny: number of patches down dip
355 :type ny: optional, int
356 :param depth: depth of the rupture plane center in [m]
357 :type depth: optional, float
358 :param component: choice, which component of the traction is displayed.
359 Choose between:
360 - "tx": along strike tractions
361 - "ty": up dip tractions
362 - "tz": normal tractions
363 - "absolut": length of the traction vector
364 :type component: optional, str
365 '''
366 import matplotlib.pyplot as plt
367 from pyrocko.modelling.okada import OkadaSource
369 comp2idx = dict(
370 tx=0, ty=1, tz=2)
372 source = OkadaSource(
373 lat=0.,
374 lon=0.,
375 depth=depth,
376 al1=-20*km, al2=20*km,
377 aw1=-15*km, aw2=15*km,
378 strike=120., dip=90., rake=90.,
379 slip=5.)
381 patches, _ = source.discretize(nx, ny)
382 tractions = tractions.get_tractions(nx, ny, patches)
384 if component in comp2idx:
385 tractions = tractions[:, comp2idx[component]].reshape(nx, ny)
386 elif component == 'absolut':
387 tractions = num.linalg.norm(tractions, axis=1).reshape(nx, ny)
388 else:
389 raise ValueError('given component is not valid.')
391 fig = plt.figure()
392 ax = fig.gca()
394 ax.imshow(tractions)
396 plt.show()
399if __name__ == '__main__':
400 tractions = TractionComposition(
401 components=[
402 UniformTractions(traction=45e3),
403 RectangularTaper(),
404 DepthTaper(depth_start=10.*km, depth_stop=30.*km)
405 ])
407 plot_tractions(tractions)
410__all__ = [
411 'AbstractTractionField',
412 'TractionField',
413 'TractionComposition',
414 'UniformTractions',
415 'HomogeneousTractions',
416 'DirectedTractions',
417 'FractalTractions',
418 'SelfSimilarTractions',
419 'RectangularTaper',
420 'DepthTaper',
421 'plot_tractions']