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 .guts import Float, Object
11from . import ahfullgreen_ext as ext
14guts_prefix = 'pf'
17class AhfullgreenError(Exception):
18 pass
21def make_seismogram(
22 vp, vs, density, qp, qs, x, f, m6,
23 quantity, deltat, stf=None, wanted_components='ned',
24 want_far=True, want_intermediate=True, want_near=True,
25 npad_levelling=40, out_alignment=0.):
27 if stf is None:
28 stf = AhfullgreenSTFImpulse()
30 x = num.asarray(x, float)
31 f = num.asarray(f, float)
32 m6 = num.asarray(m6, float)
34 r = math.sqrt(num.sum(x**2))
36 tp = r / vp
37 ts = r / vs
39 if ts <= tp:
40 raise AhfullgreenError('Unsupported material properties: ts <= tp')
42 tpad = stf.t_cutoff() or deltat * 10.
44 tstart = tp - tpad - npad_levelling * deltat
45 tstart = out_alignment + round((tstart - out_alignment) / deltat) * deltat
47 nt = trace.nextpow2(int(math.ceil(
48 (ts - tp + 2 * tpad + 2*npad_levelling * deltat) / deltat)))
50 nspec = nt // 2 + 1
52 specs = []
53 for component in 'ned':
54 if component in wanted_components:
55 specs.append(num.zeros(nspec, dtype=complex))
56 else:
57 specs.append(None)
59 oc_c = {
60 'displacement': 1, # treated in post processing
61 'velocity': 1,
62 'acceleration': 2}[quantity]
64 out_spec_delta = float(2.0 * math.pi / (nt*deltat))
65 out_spec_offset = 0.0
67 omega = out_spec_offset + out_spec_delta * num.arange(nspec)
68 coeffs_stf = stf(omega/(2.*math.pi)).astype(complex)
69 coeffs_stf *= num.exp(1.0j * omega * tstart)
71 omega_max = 2.0 * math.pi * 0.5 / deltat
72 omega_cut = omega_max * 0.75
73 icut = int(num.ceil((omega_cut - out_spec_offset) / out_spec_delta))
75 coeffs_stf[icut:] *= 0.5 + 0.5 * num.cos(
76 math.pi * num.minimum(
77 1.0, (omega[icut:] - omega_cut) / (omega_max - omega_cut)))
79 ext.add_seismogram(
80 float(vp), float(vs), float(density), float(qp), float(qs),
81 x, f, m6, oc_c, out_spec_delta, out_spec_offset,
82 specs[0], specs[1], specs[2], want_far, want_intermediate, want_near)
84 outs = []
85 for i, component in enumerate('ned'):
86 if component not in wanted_components:
87 outs.append(None)
89 out = num.fft.irfft(coeffs_stf * specs[i], nt)
90 out /= deltat
91 assert out.size // 2 + 1 == specs[i].size
93 m1 = num.mean(
94 out[:npad_levelling] * num.linspace(1., 0., npad_levelling))
96 out -= m1 * 2.
98 if quantity == 'displacement':
99 out = num.cumsum(out) * deltat
101 outs.append(out)
103 outs_wanted = []
104 for component in wanted_components:
105 i = 'ned'.find(component)
106 if i != -1:
107 outs_wanted.append(outs[i])
108 else:
109 outs_wanted.append(None)
111 return tstart, outs_wanted
114def add_seismogram(
115 vp, vs, density, qp, qs, x, f, m6,
116 quantity, deltat, out_offset,
117 out_n, out_e, out_d, stf=None,
118 want_far=True, want_intermediate=True, want_near=True,
119 npad_levelling=40):
121 ns = [out.size for out in (out_n, out_e, out_d) if out is not None]
123 if not all(n == ns[0] for n in ns):
124 raise AhfullgreenError('Length of component arrays are not identical.')
126 n = ns[0]
128 wanted_components = ''.join(
129 (c if out is not None else '-')
130 for (out, c) in zip((out_n, out_e, out_d), 'ned'))
132 tstart, temps = make_seismogram(
133 vp, vs, density, qp, qs, x, f, m6,
134 quantity, deltat, stf=stf,
135 wanted_components=wanted_components,
136 want_far=want_far,
137 want_intermediate=want_intermediate,
138 want_near=want_near,
139 npad_levelling=npad_levelling, out_alignment=out_offset)
141 for i, out in enumerate((out_n, out_e, out_d)):
142 if out is None:
143 continue
145 temp = temps[i]
147 ntemp = temp.size
149 tmin = max(out_offset, tstart)
150 tmax = min(
151 out_offset + (n-1) * deltat,
152 tstart + (ntemp-1) * deltat)
154 def ind(t, t0):
155 return int(round((t-t0)/deltat))
157 out[ind(tmin, out_offset):ind(tmax, out_offset)+1] \
158 += temp[ind(tmin, tstart):ind(tmax, tstart)+1]
160 out[:ind(tmin, out_offset)] += 0.
161 out[ind(tmax, out_offset)+1:] += temp[ind(tmax, tstart)]
164class AhfullgreenSTF(Object):
165 pass
168class AhfullgreenSTFImpulse(AhfullgreenSTF):
170 def t_cutoff(self):
171 return None
173 def __call__(self, f):
174 return num.ones(f.size, dtype=complex)
177class AhfullgreenSTFGauss(AhfullgreenSTF):
179 tau = Float.T(default=1.0)
181 def t_cutoff(self):
182 return self.tau * 2.
184 def __call__(self, f):
185 omega = f * 2. * math.pi
187 return num.exp(-(omega**2 * self.tau**2 / 8.))