1import logging 

2import numpy as num 

3from pyrocko.guts import Object, Float, List, StringChoice, Int 

4from pyrocko.guts_array import Array 

5 

6logger = logging.getLogger('pyrocko.gf.tractions') 

7km = 1e3 

8d2r = num.pi/180. 

9r2d = 180./num.pi 

10 

11 

12def tukey_window(N, alpha): 

13 assert alpha <= 1. 

14 window = num.ones(N) 

15 n = num.arange(N) 

16 

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 

21 

22 

23def planck_window(N, epsilon): 

24 assert epsilon <= 1. 

25 window = num.ones(N) 

26 n = num.arange(N) 

27 

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 

34 

35 

36class AbstractTractionField(Object): 

37 ''' 

38 Base class for multiplicative traction fields (tapers). 

39 

40 Fields of this type a re multiplied in the 

41 :py:class:`~pyrocko.gf.tractions.TractionComposition` 

42 ''' 

43 operation = 'mult' 

44 

45 def get_tractions(self, nx, ny, patches): 

46 raise NotImplementedError 

47 

48 

49class TractionField(AbstractTractionField): 

50 ''' 

51 Base class for additive traction fields. 

52 

53 Fields of this type are added in the 

54 :py:class:`~pyrocko.gf.tractions.TractionComposition` 

55 ''' 

56 operation = 'add' 

57 

58 def get_tractions(self, nx, ny, patches): 

59 raise NotImplementedError 

60 

61 

62class TractionComposition(TractionField): 

63 ''' 

64 Composition of traction fields. 

65 

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.') 

74 

75 def get_tractions(self, nx, ny, patches=None): 

76 npatches = nx * ny 

77 tractions = num.zeros((npatches, 3)) 

78 

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)) 

88 

89 return tractions 

90 

91 def add_component(self, field): 

92 logger.debug('Adding traction component.') 

93 self.components.append(field) 

94 

95 

96class HomogeneousTractions(TractionField): 

97 ''' 

98 Homogeneous traction field. 

99 

100 The traction vectors in strike, dip and normal direction are acting 

101 homogeneously on the rupture plane. 

102 ''' 

103 

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].') 

113 

114 def get_tractions(self, nx, ny, patches=None): 

115 npatches = nx * ny 

116 

117 return num.tile( 

118 (self.strike, self.dip, self.normal), npatches) \ 

119 .reshape(-1, 3) 

120 

121 

122class DirectedTractions(TractionField): 

123 ''' 

124 Directed traction field. 

125 

126 The traction vectors are following a uniform ``rake``. 

127 ''' 

128 

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].') 

138 

139 def get_tractions(self, nx, ny, patches=None): 

140 npatches = nx * ny 

141 

142 strike = num.cos(self.rake*d2r) * self.traction 

143 dip = num.sin(self.rake*d2r) * self.traction 

144 normal = 0. 

145 

146 return num.tile((strike, dip, normal), npatches).reshape(-1, 3) 

147 

148 

149class SelfSimilarTractions(TractionField): 

150 ''' 

151 Traction model following Power & Tullis (1991). 

152 

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``. 

162 

163 ''' 

164 rank = Int.T( 

165 default=1, 

166 help='Maximum summed cosine wavenumber/spatial frequency.') 

167 

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.') 

174 

175 traction_max = Float.T( 

176 default=1., 

177 help='Maximum traction vector length [Pa].') 

178 

179 phases = Array.T( 

180 optional=True, 

181 dtype=num.float, 

182 shape=(None,), 

183 help='Phase shift of the cosines in [rad].') 

184 

185 def get_phases(self): 

186 if self.phases is not None: 

187 if self.phases.shape[0] == self.rank: 

188 return self.phases 

189 

190 return (num.random.random(self.rank) * 2. - 1.) * num.pi 

191 

192 def get_tractions(self, nx, ny, patches=None): 

193 z = num.zeros((ny, nx)) 

194 phases = self.get_phases() 

195 

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]) 

202 

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') 

206 

207 t *= self.traction_max / num.max(num.linalg.norm(t, axis=1)) 

208 

209 return t 

210 

211 

212class FractalTractions(TractionField): 

213 ''' 

214 Fractal traction field. 

215 ''' 

216 

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.') 

222 

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.') 

229 

230 traction_max = Float.T( 

231 default=1., 

232 help='Maximum traction vector length [Pa].') 

233 

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 

239 

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) 

244 

245 return self._data 

246 

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 

254 

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 

259 

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) 

264 

265 # Define wavenumber bins 

266 k_bins = num.arange(0, num.max(k_rad), num.max(k_rad)/10.) 

267 

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. 

271 

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) 

278 

279 amps[k_rad > k_bins.max()] = power_spec.ravel()[num.argmax(power_spec)] 

280 

281 # Multiply spectrum by amplitudes and inverse fft into demeaned noise 

282 spec *= amps.T 

283 

284 tractions = num.abs(num.fft.ifft2(spec)) 

285 tractions -= num.mean(tractions) 

286 tractions *= self.traction_max / num.abs(tractions).max() 

287 

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') 

291 

292 return t 

293 

294 

295class RectangularTaper(AbstractTractionField): 

296 width = Float.T( 

297 default=.2, 

298 help='Width of the taper as a fraction of the plane.') 

299 

300 type = StringChoice.T( 

301 choices=('tukey', ), 

302 default='tukey', 

303 help='Type of the taper, default: "tukey".') 

304 

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] 

310 

311 raise AttributeError('Unknown type: %s' % self.type) 

312 

313 

314class DepthTaper(AbstractTractionField): 

315 depth_start = Float.T( 

316 help='Depth where the taper begins [m].') 

317 

318 depth_stop = Float.T( 

319 help='Depth where taper ends and drops to zero [m].') 

320 

321 type = StringChoice.T( 

322 choices=('linear', ), 

323 default='linear', 

324 help='Type of the taper, default: "linear".') 

325 

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]) 

329 

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] 

337 

338 

339def plot_tractions(tractions, nx=15, ny=12, depth=10*km, component='strike'): 

340 ''' 

341 Plot traction model for quick inspection. 

342 

343 :param tractions: 

344 Traction field or traction composition to be displayed. 

345 :type tractions: 

346 :py:class:`pyrocko.gf.tractions.TractionField` 

347 

348 :param nx: 

349 Number of patches along strike. 

350 :type nx: 

351 optional, int 

352 

353 :param ny: 

354 Number of patches down dip. 

355 :type ny: 

356 optional, int 

357 

358 :param depth: 

359 Depth of the rupture plane center in [m]. 

360 :type depth: 

361 optional, float 

362 

363 :param component: 

364 Choice of traction component to be shown. Available: ``'tx'`` (along 

365 strike), ``'ty'`` (up dip), ``'tz'`` (normal), ``'absolute'`` (vector 

366 length). 

367 :type component: 

368 optional, str 

369 ''' 

370 import matplotlib.pyplot as plt 

371 from pyrocko.modelling.okada import OkadaSource 

372 

373 comp2idx = dict( 

374 tx=0, ty=1, tz=2) 

375 

376 source = OkadaSource( 

377 lat=0., 

378 lon=0., 

379 depth=depth, 

380 al1=-20*km, al2=20*km, 

381 aw1=-15*km, aw2=15*km, 

382 strike=120., dip=90., rake=90., 

383 slip=5.) 

384 

385 patches, _ = source.discretize(nx, ny) 

386 tractions = tractions.get_tractions(nx, ny, patches) 

387 

388 if component in comp2idx: 

389 tractions = tractions[:, comp2idx[component]].reshape(nx, ny) 

390 elif component == 'absolute': 

391 tractions = num.linalg.norm(tractions, axis=1).reshape(nx, ny) 

392 else: 

393 raise ValueError('Given component is not valid.') 

394 

395 fig = plt.figure() 

396 ax = fig.gca() 

397 

398 ax.imshow(tractions) 

399 

400 plt.show() 

401 

402 

403__all__ = [ 

404 'AbstractTractionField', 

405 'TractionField', 

406 'TractionComposition', 

407 'HomogeneousTractions', 

408 'DirectedTractions', 

409 'FractalTractions', 

410 'SelfSimilarTractions', 

411 'RectangularTaper', 

412 'DepthTaper', 

413 'plot_tractions']