1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import math
7import logging
8import numpy as num
9from . import trace
10from .guts import Float, Object
11from . import ahfullgreen_ext as ext
13logger = logging.getLogger('pyrocko.fomosto.ahfullgreen')
16guts_prefix = 'pf'
19class AhfullgreenError(Exception):
20 pass
23def make_seismogram(
24 vp, vs, density, qp, qs, x, f, m6,
25 quantity, deltat, stf=None, wanted_components='ned',
26 want_far=True, want_intermediate=True, want_near=True,
27 npad_levelling=40, out_alignment=0.):
29 if stf is None:
30 stf = AhfullgreenSTFImpulse()
32 x = num.asarray(x, float)
33 f = num.asarray(f, float)
34 m6 = num.asarray(m6, float)
36 r = math.sqrt(num.sum(x**2))
38 tp = r / vp
39 ts = r / vs
41 if ts < tp:
42 raise AhfullgreenError('unsupported material properties: ts < tp')
44 tpad = stf.t_cutoff() or deltat * 10.
46 tstart = tp - tpad - npad_levelling * deltat
47 tstart = out_alignment + round((tstart - out_alignment) / deltat) * deltat
49 nt = trace.nextpow2(int(math.ceil(
50 (ts - tp + 2 * tpad + 2*npad_levelling * deltat) / deltat)))
52 nspec = nt // 2 + 1
54 specs = []
55 for component in 'ned':
56 if component in wanted_components:
57 specs.append(num.zeros(nspec, dtype=complex))
58 else:
59 specs.append(None)
61 oc_c = {
62 'displacement': 1, # treated in post processing
63 'velocity': 1,
64 'acceleration': 2}[quantity]
66 out_spec_delta = float(2.0 * math.pi / (nt*deltat))
67 out_spec_offset = 0.0
69 omega = out_spec_offset + out_spec_delta * num.arange(nspec)
70 coeffs_stf = stf(omega/(2.*math.pi)).astype(complex)
71 coeffs_stf *= num.exp(1.0j * omega * tstart)
73 omega_max = 2.0 * math.pi * 0.5 / deltat
74 omega_cut = omega_max * 0.75
75 icut = int(num.ceil((omega_cut - out_spec_offset) / out_spec_delta))
77 coeffs_stf[icut:] *= 0.5 + 0.5 * num.cos(
78 math.pi * num.minimum(
79 1.0, (omega[icut:] - omega_cut) / (omega_max - omega_cut)))
81 if num.all(x == 0.0):
82 logger.warn(
83 'Source and receiver are at the same position -> setting GF for '
84 'this combination to zero.')
85 else:
86 ext.add_seismogram(
87 float(vp), float(vs), float(density), float(qp), float(qs),
88 x, f, m6, oc_c, out_spec_delta, out_spec_offset,
89 specs[0], specs[1], specs[2],
90 want_far, want_intermediate, want_near)
92 outs = []
93 for i, component in enumerate('ned'):
94 if component not in wanted_components:
95 outs.append(None)
97 out = num.fft.irfft(coeffs_stf * specs[i], nt)
98 out /= deltat
99 assert out.size // 2 + 1 == specs[i].size
101 m1 = num.mean(
102 out[:npad_levelling] * num.linspace(1., 0., npad_levelling))
104 out -= m1 * 2.
106 if quantity == 'displacement':
107 out = num.cumsum(out) * deltat
109 outs.append(out)
111 outs_wanted = []
112 for component in wanted_components:
113 i = 'ned'.find(component)
114 if i != -1:
115 outs_wanted.append(outs[i])
116 else:
117 outs_wanted.append(None)
119 return tstart, outs_wanted
122def add_seismogram(
123 vp, vs, density, qp, qs, x, f, m6,
124 quantity, deltat, out_offset,
125 out_n, out_e, out_d, stf=None,
126 want_far=True, want_intermediate=True, want_near=True,
127 npad_levelling=40):
129 ns = [out.size for out in (out_n, out_e, out_d) if out is not None]
131 if not all(n == ns[0] for n in ns):
132 raise AhfullgreenError('Length of component arrays are not identical.')
134 n = ns[0]
136 wanted_components = ''.join(
137 (c if out is not None else '-')
138 for (out, c) in zip((out_n, out_e, out_d), 'ned'))
140 tstart, temps = make_seismogram(
141 vp, vs, density, qp, qs, x, f, m6,
142 quantity, deltat, stf=stf,
143 wanted_components=wanted_components,
144 want_far=want_far,
145 want_intermediate=want_intermediate,
146 want_near=want_near,
147 npad_levelling=npad_levelling, out_alignment=out_offset)
149 for i, out in enumerate((out_n, out_e, out_d)):
150 if out is None:
151 continue
153 temp = temps[i]
155 ntemp = temp.size
157 tmin = max(out_offset, tstart)
158 tmax = min(
159 out_offset + (n-1) * deltat,
160 tstart + (ntemp-1) * deltat)
162 def ind(t, t0):
163 return int(round((t-t0)/deltat))
165 out[ind(tmin, out_offset):ind(tmax, out_offset)+1] \
166 += temp[ind(tmin, tstart):ind(tmax, tstart)+1]
168 out[:ind(tmin, out_offset)] += 0.
169 out[ind(tmax, out_offset)+1:] += temp[ind(tmax, tstart)]
172class AhfullgreenSTF(Object):
173 pass
176class AhfullgreenSTFImpulse(AhfullgreenSTF):
178 def t_cutoff(self):
179 return None
181 def __call__(self, f):
182 return num.ones(f.size, dtype=complex)
185class AhfullgreenSTFGauss(AhfullgreenSTF):
187 tau = Float.T(default=1.0)
189 def t_cutoff(self):
190 return self.tau * 2.
192 def __call__(self, f):
193 omega = f * 2. * math.pi
195 return num.exp(-(omega**2 * self.tau**2 / 8.))