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 numpy as num
9from . import trace
10from . import ahfullgreen_ext as ext
13class AhfullgreenError(Exception):
14 pass
17def make_seismogram(
18 vp, vs, density, qp, qs, x, f, m6,
19 quantity, deltat, stf=None, wanted_components='ned',
20 want_far=True, want_intermediate=True, want_near=True,
21 npad_levelling=40, out_alignment=0.):
23 if stf is None:
24 stf = Impulse()
26 x = num.asarray(x, float)
27 f = num.asarray(f, float)
28 m6 = num.asarray(m6, float)
30 r = math.sqrt(num.sum(x**2))
32 tp = r / vp
33 ts = r / vs
35 if ts <= tp:
36 raise AhfullgreenError('Unsupported material properties: ts <= tp')
38 tpad = stf.t_cutoff() or deltat * 10.
40 tstart = tp - tpad - npad_levelling * deltat
41 tstart = out_alignment + round((tstart - out_alignment) / deltat) * deltat
43 nt = trace.nextpow2(int(math.ceil(
44 (ts - tp + 2 * tpad + 2*npad_levelling * deltat) / deltat)))
46 nspec = nt // 2 + 1
48 specs = []
49 for component in 'ned':
50 if component in wanted_components:
51 specs.append(num.zeros(nspec, dtype=complex))
52 else:
53 specs.append(None)
55 oc_c = {
56 'displacement': 1, # treated in post processing
57 'velocity': 1,
58 'acceleration': 2}[quantity]
60 out_spec_delta = float(2.0 * math.pi / (nt*deltat))
61 out_spec_offset = 0.0
63 omega = out_spec_offset + out_spec_delta * num.arange(nspec)
64 coeffs_stf = stf(omega/(2.*math.pi)).astype(complex)
65 coeffs_stf *= num.exp(1.0j * omega * tstart)
67 omega_max = 2.0 * math.pi * 0.5 / deltat
68 omega_cut = omega_max * 0.75
69 icut = int(num.ceil((omega_cut - out_spec_offset) / out_spec_delta))
71 coeffs_stf[icut:] *= 0.5 + 0.5 * num.cos(
72 math.pi * num.minimum(
73 1.0, (omega[icut:] - omega_cut) / (omega_max - omega_cut)))
75 ext.add_seismogram(
76 float(vp), float(vs), float(density), float(qp), float(qs),
77 x, f, m6, oc_c, out_spec_delta, out_spec_offset,
78 specs[0], specs[1], specs[2], want_far, want_intermediate, want_near)
80 outs = []
81 for i, component in enumerate('ned'):
82 if component not in wanted_components:
83 outs.append(None)
85 out = num.fft.irfft(coeffs_stf * specs[i], nt)
86 out /= deltat
87 assert out.size // 2 + 1 == specs[i].size
89 m1 = num.mean(
90 out[:npad_levelling] * num.linspace(1., 0., npad_levelling))
92 out -= m1 * 2.
94 if quantity == 'displacement':
95 out = num.cumsum(out) * deltat
97 outs.append(out)
99 outs_wanted = []
100 for component in wanted_components:
101 i = 'ned'.find(component)
102 if i != -1:
103 outs_wanted.append(outs[i])
104 else:
105 outs_wanted.append(None)
107 return tstart, outs_wanted
110def add_seismogram(
111 vp, vs, density, qp, qs, x, f, m6,
112 quantity, deltat, out_offset,
113 out_n, out_e, out_d, stf=None,
114 want_far=True, want_intermediate=True, want_near=True,
115 npad_levelling=40):
117 ns = [out.size for out in (out_n, out_e, out_d) if out is not None]
119 if not all(n == ns[0] for n in ns):
120 raise AhfullgreenError('Length of component arrays are not identical.')
122 n = ns[0]
124 wanted_components = ''.join(
125 (c if out is not None else '-')
126 for (out, c) in zip((out_n, out_e, out_d), 'ned'))
128 tstart, temps = make_seismogram(
129 vp, vs, density, qp, qs, x, f, m6,
130 quantity, deltat, stf=stf,
131 wanted_components=wanted_components,
132 want_far=want_far,
133 want_intermediate=want_intermediate,
134 want_near=want_near,
135 npad_levelling=npad_levelling, out_alignment=out_offset)
137 for i, out in enumerate((out_n, out_e, out_d)):
138 if out is None:
139 continue
141 temp = temps[i]
143 ntemp = temp.size
145 tmin = max(out_offset, tstart)
146 tmax = min(
147 out_offset + (n-1) * deltat,
148 tstart + (ntemp-1) * deltat)
150 def ind(t, t0):
151 return int(round((t-t0)/deltat))
153 out[ind(tmin, out_offset):ind(tmax, out_offset)+1] \
154 += temp[ind(tmin, tstart):ind(tmax, tstart)+1]
156 out[:ind(tmin, out_offset)] += 0.
157 out[ind(tmax, out_offset)+1:] += temp[ind(tmax, tstart)]
160class Impulse(object):
162 def __init__(self):
163 pass
165 def t_cutoff(self):
166 return None
168 def __call__(self, f):
169 return num.ones(f.size, dtype=complex)
172class Gauss(object):
173 def __init__(self, tau):
174 self._tau = tau
176 def t_cutoff(self):
177 return self._tau * 2.
179 def __call__(self, f):
180 omega = f * 2. * math.pi
182 return num.exp(-(omega**2 * self._tau**2 / 8.))