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