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
raise NotImplementedError
raise NotImplementedError
AbstractTractionField.T(), default=[], help='Ordered list of tractions')
else: raise AttributeError( 'Component %s has an invalid operation %s.' % (comp, comp.operation))
logger.debug('adding traction component') self.components.append(field)
default=1., help='Uniform traction in strike, dip and normal direction [Pa]')
npatches = nx * ny return num.full((npatches, 3), self.traction)
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)
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].')
return (num.random.random(self.rank) * 2. - 1.) * num.pi
default=None, optional=True, help='index of the numpy random state used. If None, an arbitrary ' 'random state is 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]')
raise AttributeError( 'patches needs to be given for this traction field')
# Create random data and get spectrum and power spectrum
# Get 0-centered wavenumbers (k_rad == 0.) is in the centre
# Define wavenumber bins
# Set amplitudes within wavenumber bins to power spec * 1 / k_max
# Multiply spectrum by amplitudes and inverse fft into demeaned noise
default=.2, help='Width of the taper as a fraction of the plane.') choices=('tukey', ), default='tukey', help='Type of the taper, default "tukey"')
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]
import matplotlib.pyplot as plt from pyrocko.modelling.okada import OkadaSource
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) tractions = tractions[:, 0].reshape(nx, ny)
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) |