1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division
7import math
8import logging
9import numpy as num
10from . import trace
11from .guts import Float, Object
12from . import ahfullgreen_ext as ext
14logger = logging.getLogger('pyrocko.fomosto.ahfullgreen')
17guts_prefix = 'pf'
20class AhfullgreenError(Exception):
21 pass
24def make_seismogram(
25 vp, vs, density, qp, qs, x, f, m6,
26 quantity, deltat, stf=None, wanted_components='ned',
27 want_far=True, want_intermediate=True, want_near=True,
28 npad_levelling=40, out_alignment=0.):
30 if stf is None:
31 stf = AhfullgreenSTFImpulse()
33 x = num.asarray(x, float)
34 f = num.asarray(f, float)
35 m6 = num.asarray(m6, float)
37 r = math.sqrt(num.sum(x**2))
39 tp = r / vp
40 ts = r / vs
42 if ts < tp:
43 raise AhfullgreenError('unsupported material properties: ts < tp')
45 tpad = stf.t_cutoff() or deltat * 10.
47 tstart = tp - tpad - npad_levelling * deltat
48 tstart = out_alignment + round((tstart - out_alignment) / deltat) * deltat
50 nt = trace.nextpow2(int(math.ceil(
51 (ts - tp + 2 * tpad + 2*npad_levelling * deltat) / deltat)))
53 nspec = nt // 2 + 1
55 specs = []
56 for component in 'ned':
57 if component in wanted_components:
58 specs.append(num.zeros(nspec, dtype=complex))
59 else:
60 specs.append(None)
62 oc_c = {
63 'displacement': 1, # treated in post processing
64 'velocity': 1,
65 'acceleration': 2}[quantity]
67 out_spec_delta = float(2.0 * math.pi / (nt*deltat))
68 out_spec_offset = 0.0
70 omega = out_spec_offset + out_spec_delta * num.arange(nspec)
71 coeffs_stf = stf(omega/(2.*math.pi)).astype(complex)
72 coeffs_stf *= num.exp(1.0j * omega * tstart)
74 omega_max = 2.0 * math.pi * 0.5 / deltat
75 omega_cut = omega_max * 0.75
76 icut = int(num.ceil((omega_cut - out_spec_offset) / out_spec_delta))
78 coeffs_stf[icut:] *= 0.5 + 0.5 * num.cos(
79 math.pi * num.minimum(
80 1.0, (omega[icut:] - omega_cut) / (omega_max - omega_cut)))
82 if num.all(x == 0.0):
83 logger.warn(
84 'Source and receiver are at the same position -> setting GF for '
85 'this combination to zero.')
86 else:
87 ext.add_seismogram(
88 float(vp), float(vs), float(density), float(qp), float(qs),
89 x, f, m6, oc_c, out_spec_delta, out_spec_offset,
90 specs[0], specs[1], specs[2],
91 want_far, want_intermediate, want_near)
93 outs = []
94 for i, component in enumerate('ned'):
95 if component not in wanted_components:
96 outs.append(None)
98 out = num.fft.irfft(coeffs_stf * specs[i], nt)
99 out /= deltat
100 assert out.size // 2 + 1 == specs[i].size
102 m1 = num.mean(
103 out[:npad_levelling] * num.linspace(1., 0., npad_levelling))
105 out -= m1 * 2.
107 if quantity == 'displacement':
108 out = num.cumsum(out) * deltat
110 outs.append(out)
112 outs_wanted = []
113 for component in wanted_components:
114 i = 'ned'.find(component)
115 if i != -1:
116 outs_wanted.append(outs[i])
117 else:
118 outs_wanted.append(None)
120 return tstart, outs_wanted
123def add_seismogram(
124 vp, vs, density, qp, qs, x, f, m6,
125 quantity, deltat, out_offset,
126 out_n, out_e, out_d, stf=None,
127 want_far=True, want_intermediate=True, want_near=True,
128 npad_levelling=40):
130 ns = [out.size for out in (out_n, out_e, out_d) if out is not None]
132 if not all(n == ns[0] for n in ns):
133 raise AhfullgreenError('Length of component arrays are not identical.')
135 n = ns[0]
137 wanted_components = ''.join(
138 (c if out is not None else '-')
139 for (out, c) in zip((out_n, out_e, out_d), 'ned'))
141 tstart, temps = make_seismogram(
142 vp, vs, density, qp, qs, x, f, m6,
143 quantity, deltat, stf=stf,
144 wanted_components=wanted_components,
145 want_far=want_far,
146 want_intermediate=want_intermediate,
147 want_near=want_near,
148 npad_levelling=npad_levelling, out_alignment=out_offset)
150 for i, out in enumerate((out_n, out_e, out_d)):
151 if out is None:
152 continue
154 temp = temps[i]
156 ntemp = temp.size
158 tmin = max(out_offset, tstart)
159 tmax = min(
160 out_offset + (n-1) * deltat,
161 tstart + (ntemp-1) * deltat)
163 def ind(t, t0):
164 return int(round((t-t0)/deltat))
166 out[ind(tmin, out_offset):ind(tmax, out_offset)+1] \
167 += temp[ind(tmin, tstart):ind(tmax, tstart)+1]
169 out[:ind(tmin, out_offset)] += 0.
170 out[ind(tmax, out_offset)+1:] += temp[ind(tmax, tstart)]
173class AhfullgreenSTF(Object):
174 pass
177class AhfullgreenSTFImpulse(AhfullgreenSTF):
179 def t_cutoff(self):
180 return None
182 def __call__(self, f):
183 return num.ones(f.size, dtype=complex)
186class AhfullgreenSTFGauss(AhfullgreenSTF):
188 tau = Float.T(default=1.0)
190 def t_cutoff(self):
191 return self.tau * 2.
193 def __call__(self, f):
194 omega = f * 2. * math.pi
196 return num.exp(-(omega**2 * self.tau**2 / 8.))