|
assert alpha <= 1. window = num.ones(N) n = num.arange(N)
N_f = int((alpha * N)//2) window[:N_f] = .5 * (1. - num.cos((2*num.pi * n[:N_f])/(alpha * N))) window[(N-N_f):] = window[:N_f][::-1] return window
assert epsilon <= 1. window = num.ones(N) n = num.arange(N)
N_f = int((epsilon * N)) window[:N_f] = \ (1. + num.exp((epsilon * N) / n[:N_f] - (epsilon * N) / ((epsilon * N - n[:N_f]))))**-1. window[(N-N_f):] = window[:N_f][::-1] return window
''' Abstract traction field
Fields of this type a re multiplied in the :py:class:`~pyrocko.gf.tractions.TractionComposition` '''
raise NotImplementedError
''' Traction field
Fields of this type are added in the :py:class:`~pyrocko.gf.tractions.TractionComposition` '''
raise NotImplementedError
''' Composition of traction fields
:py:class:`~pyrocko.gf.tractions.TractionField` and :py:class:`~pyrocko.gf.tractions.AbstractTractionField` can be combined to realize a combination of different fields. ''' AbstractTractionField.T(), default=[], help='Ordered list of tractions')
npatches = nx * ny tractions = num.zeros((npatches, 3))
for comp in self.components: if comp.operation == 'add': tractions += comp.get_tractions(nx, ny, patches) elif comp.operation == 'mult': tractions *= comp.get_tractions(nx, ny, patches) else: raise AttributeError( 'Component %s has an invalid operation %s.' % (comp, comp.operation))
return tractions
logger.debug('adding traction component') self.components.append(field)
''' Uniform traction field
The traction field is uniform in strike, dip and normal direction. This realisation is not only simple but also unrealistic. ''' default=1., help='Uniform traction in strike, dip and normal direction [Pa]')
npatches = nx * ny return num.full((npatches, 3), self.traction)
''' Homogeneous traction field
The traction vectors in strike, dip and normal direction are acting homogeneously on the rupture plane. '''
default=1., help='Tractions in strike direction [Pa]') default=1., help='Traction in dip direction (up) [Pa]') default=1., help='Traction in normal direction [Pa]')
(self.strike, self.dip, self.normal), npatches) \ .reshape(-1, 3)
''' Directed traction field
The traction vectors are following a uniform ``rake``. '''
default=0., help='rake angle in [deg], ' 'measured counter-clockwise from right-horizontal ' 'in on-plane view. Rake is translated into homogenous tractions ' 'in strike and up-dip direction.') default=1., help='Traction in rake direction [Pa]')
''' Traction model following Power & Tullis (1991).
The traction vectors are calculated as a sum of 2D-cosines with a constant amplitude / wavelength ratio. The wavenumber kx and ky are constant for each cosine function. The rank defines the maximum wavenumber used for summation. So, e.g. a rank of 3 will lead to a summation of cosines with ``kx = ky`` in (1, 2, 3). Each cosine has an associated phases, which defines both the phase shift and also the shift from the rupture plane centre. Finally the summed cosines are translated into shear tractions based on the rake and normalized with ``traction_max``.
''' default=1, help='maximum summed cosine wavenumber/spatial frequency.')
default=0., help='rake angle in [deg], ' 'measured counter-clockwise from right-horizontal ' 'in on-plane view. Rake is translated into homogenous tractions ' 'in strike and up-dip direction.')
default=1., help='maximum traction vector length [Pa]')
optional=True, dtype=num.float, shape=(None,), help='phase shift of the cosines in [rad].')
if self.phases is not None: if self.phases.shape[0] == self.rank: return self.phases
return (num.random.random(self.rank) * 2. - 1.) * num.pi
z = num.zeros((ny, nx)) phases = self.get_phases()
for i in range(1, self.rank+1): x = num.linspace(-i*num.pi, i*num.pi, nx) + i*phases[i-1] y = num.linspace(-i*num.pi, i*num.pi, ny) + i*phases[i-1] x, y = num.meshgrid(x, y) r = num.sqrt(x**2 + y**2) z += 1. / i * num.cos(r + phases[i-1])
t = num.zeros((nx*ny, 3)) t[:, 0] = num.cos(self.rake*d2r) * z.ravel(order='F') t[:, 1] = num.sin(self.rake*d2r) * z.ravel(order='F')
t *= self.traction_max / num.max(num.linalg.norm(t, axis=1))
return t
''' Fractal traction field '''
default=None, optional=True, help='Seed for :py:class:`~numpy.random.RandomState`.' 'If ``None``, an random seed will be initialized.')
default=0., help='rake angle in [deg], ' 'measured counter-clockwise from right-horizontal ' 'in on-plane view. Rake is translated into homogenous tractions ' 'in strike and up-dip direction.')
default=1., help='maximum traction vector length [Pa]')
super().__init__(*args, **kwargs) if self.rseed is None: self.rseed = num.random.randint(0, 2**32-1) self._data = None
if self._data is None: rstate = num.random.RandomState(self.rseed) self._data = rstate.rand(nx, ny)
return self._data
if patches is None: raise AttributeError( 'patches needs to be given for this traction field') npatches = nx * ny dx = -patches[0].al1 + patches[0].al2 dy = -patches[0].aw1 + patches[0].aw2
# Create random data and get spectrum and power spectrum data = self._get_data(nx, ny) spec = num.fft.fftshift(num.fft.fft2(data)) power_spec = (num.abs(spec)/spec.size)**2
# Get 0-centered wavenumbers (k_rad == 0.) is in the centre kx = num.fft.fftshift(num.fft.fftfreq(nx, d=dx)) ky = num.fft.fftshift(num.fft.fftfreq(ny, d=dy)) k_rad = num.sqrt(ky[:, num.newaxis]**2 + kx[num.newaxis, :]**2)
# Define wavenumber bins k_bins = num.arange(0, num.max(k_rad), num.max(k_rad)/10.)
# Set amplitudes within wavenumber bins to power_spec * 1 / k_max amps = num.zeros(k_rad.shape) amps[k_rad == 0.] = 1.
for i in range(k_bins.size-1): k_min = k_bins[i] k_max = k_bins[i+1] r = num.logical_and(k_rad > k_min, k_rad <= k_max) amps[r] = power_spec.T[r] amps = num.sqrt(amps * data.size * num.pi * 4)
amps[k_rad > k_bins.max()] = power_spec.ravel()[num.argmax(power_spec)]
# Multiply spectrum by amplitudes and inverse fft into demeaned noise spec *= amps.T
tractions = num.abs(num.fft.ifft2(spec)) tractions -= num.mean(tractions) tractions *= self.traction_max / num.abs(tractions).max()
t = num.zeros((npatches, 3)) t[:, 0] = num.cos(self.rake*d2r) * tractions.ravel(order='C') t[:, 1] = num.sin(self.rake*d2r) * tractions.ravel(order='C')
return t
default=.2, help='Width of the taper as a fraction of the plane.')
choices=('tukey', ), default='tukey', help='Type of the taper, default "tukey"')
if self.type == 'tukey': x = tukey_window(nx, self.width) y = tukey_window(ny, self.width) return (x[:, num.newaxis] * y).ravel()[:, num.newaxis]
raise AttributeError('unknown type %s' % self.type)
help='Depth where the taper begins [km]')
help='Depth where taper stops, and drops to 0. [km]')
choices=('linear', ), default='linear', help='Type of the taper, default "linear"')
assert self.depth_stop > self.depth_start depths = num.array([p.depth for p in patches])
if self.type == 'linear': slope = self.depth_stop - self.depth_start depths -= self.depth_stop depths /= -slope depths[depths > 1.] = 1. depths[depths < 0.] = 0. return depths[:, num.newaxis]
'''Plot choosen traction model for quick inspection
:param tractions: traction field or traction composition to be displayed :type tractions: :py:class:`pyrocko.gf.tractions.TractionField` :param nx: number of patches along strike :type nx: optional, int :param ny: number of patches down dip :type ny: optional, int :param depth: depth of the rupture plane center in [m] :type depth: optional, float :param component: choice, which component of the traction is displayed. Choose between: - "tx": along strike tractions - "ty": up dip tractions - "tz": normal tractions - "absolut": length of the traction vector :type component: optional, str ''' import matplotlib.pyplot as plt from pyrocko.modelling.okada import OkadaSource
comp2idx = dict( tx=0, ty=1, tz=2)
source = OkadaSource( lat=0., lon=0., depth=depth, al1=-20*km, al2=20*km, aw1=-15*km, aw2=15*km, strike=120., dip=90., rake=90., slip=5.)
patches, _ = source.discretize(nx, ny) tractions = tractions.get_tractions(nx, ny, patches)
if component in comp2idx: tractions = tractions[:, comp2idx[component]].reshape(nx, ny) elif component == 'absolut': tractions = num.linalg.norm(tractions, axis=1).reshape(nx, ny) else: raise ValueError('given component is not valid.')
fig = plt.figure() ax = fig.gca()
ax.imshow(tractions)
plt.show()
if __name__ == '__main__': tractions = TractionComposition( components=[ UniformTractions(traction=45e3), RectangularTaper(), DepthTaper(depth_start=10.*km, depth_stop=30.*km) ])
plot_tractions(tractions)
'AbstractTractionField', 'TractionField', 'TractionComposition', 'UniformTractions', 'HomogeneousTractions', 'DirectedTractions', 'FractalTractions', 'SelfSimilarTractions', 'RectangularTaper', 'DepthTaper', 'plot_tractions'] |