1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import numpy as num
8from pyrocko import util
9from pyrocko.guts import Float
10from pyrocko import moment_tensor, gf
12from .base import SourceGenerator
13from ..error import ScenarioError
15km = 1e3
16guts_prefix = 'pf.scenario'
19class DCSourceGenerator(SourceGenerator):
21 depth_min = Float.T(default=0.0)
22 depth_max = Float.T(default=30*km)
24 strike = Float.T(optional=True)
25 dip = Float.T(optional=True)
26 rake = Float.T(optional=True)
27 perturbation_angle_std = Float.T(optional=True)
29 def get_source(self, ievent):
30 rstate = self.get_rstate(ievent)
31 time = self.time_min + rstate.uniform(
32 0., float(self.time_max - self.time_min)) # hptime aware
34 lat, lon, north_shift, east_shift, depth = self.get_coordinates(ievent)
35 depth = rstate.uniform(self.depth_min, self.depth_max)
36 magnitude = self.draw_magnitude(rstate)
38 if self.strike is None and self.dip is None and self.rake is None:
39 mt = moment_tensor.MomentTensor.random_dc(x=rstate.uniform(size=3))
40 else:
41 if None in (self.strike, self.dip, self.rake):
42 raise ScenarioError(
43 'DCSourceGenerator: '
44 'strike, dip, and rake must be used in combination.')
46 mt = moment_tensor.MomentTensor(
47 strike=self.strike, dip=self.dip, rake=self.rake)
49 if self.perturbation_angle_std is not None:
50 mt = mt.random_rotated(
51 self.perturbation_angle_std,
52 rstate=rstate)
54 (s, d, r), (_, _, _) = mt.both_strike_dip_rake()
56 source = gf.DCSource(
57 name='ev%04i' % ievent,
58 time=util.to_time_float(time),
59 lat=float(lat),
60 lon=float(lon),
61 north_shift=float(north_shift),
62 east_shift=float(east_shift),
63 depth=float(depth),
64 magnitude=float(magnitude),
65 strike=float(s),
66 dip=float(d),
67 rake=float(r))
69 return source
71 def add_map_artists(self, automap):
72 from pyrocko.plot import gmtpy
74 for source in self.get_sources():
75 event = source.pyrocko_event()
76 mt = event.moment_tensor.m_up_south_east()
78 xx = num.trace(mt) / 3.
79 mc = num.array([[xx, 0., 0.], [0., xx, 0.], [0., 0., xx]])
80 mc = mt - mc
81 mc = mc / event.moment_tensor.scalar_moment() * \
82 moment_tensor.magnitude_to_moment(5.0)
83 m6 = tuple(moment_tensor.to6(mc))
84 symbol_size = 20.
85 automap.gmt.psmeca(
86 S='%s%g' % ('d', symbol_size / gmtpy.cm),
87 in_rows=[
88 (source.effective_lon, source.effective_lat, 10)
89 + m6
90 + (1, 0, 0)],
91 M=True,
92 *automap.jxyr)