|
# http://pyrocko.org - GPLv3 # # The Pyrocko Developers, 21st Century # ---|P------/S----------~Lg---------- from __future__ import absolute_import, division from builtins import zip
import math import numpy as num from . import trace from . import ahfullgreen_ext as ext
class AhfullgreenError(Exception): pass
def make_seismogram( vp, vs, density, qp, qs, x, f, m6, quantity, deltat, stf=None, wanted_components='ned', want_far=True, want_intermediate=True, want_near=True, npad_levelling=40, out_alignment=0.):
if stf is None: stf = Impulse()
x = num.asarray(x, num.float) f = num.asarray(f, num.float) m6 = num.asarray(m6, num.float)
r = math.sqrt(num.sum(x**2))
tp = r / vp ts = r / vs
if ts <= tp: raise AhfullgreenError('unsupported material properties')
tpad = stf.t_cutoff() or deltat * 10.
tstart = tp - tpad - npad_levelling * deltat tstart = out_alignment + round((tstart - out_alignment) / deltat) * deltat
nt = trace.nextpow2(int(math.ceil( (ts - tp + 2 * tpad + 2*npad_levelling * deltat) / deltat)))
nspec = nt // 2 + 1
specs = [] for component in 'ned': if component in wanted_components: specs.append(num.zeros(nspec, dtype=num.complex)) else: specs.append(None)
oc_c = { 'displacement': 1, # treated in post processing 'velocity': 1, 'acceleration': 2}[quantity]
out_spec_delta = float(2.0 * math.pi / (nt*deltat)) out_spec_offset = 0.0
omega = out_spec_offset + out_spec_delta * num.arange(nspec) coeffs_stf = stf(omega/(2.*math.pi)).astype(num.complex) coeffs_stf *= num.exp(1.0j * omega * tstart)
omega_max = 2.0 * math.pi * 0.5 / deltat omega_cut = omega_max * 0.75 icut = int(num.ceil((omega_cut - out_spec_offset) / out_spec_delta))
coeffs_stf[icut:] *= 0.5 + 0.5 * num.cos( math.pi * num.minimum( 1.0, (omega[icut:] - omega_cut) / (omega_max - omega_cut)))
ext.add_seismogram( float(vp), float(vs), float(density), float(qp), float(qs), x, f, m6, oc_c, out_spec_delta, out_spec_offset, specs[0], specs[1], specs[2], want_far, want_intermediate, want_near)
outs = [] for i, component in enumerate('ned'): if component not in wanted_components: outs.append(None)
out = num.fft.irfft(coeffs_stf * specs[i], nt) out /= deltat assert out.size // 2 + 1 == specs[i].size
m1 = num.mean( out[:npad_levelling] * num.linspace(1., 0., npad_levelling))
out -= m1 * 2.
if quantity == 'displacement': out = num.cumsum(out) * deltat
outs.append(out)
outs_wanted = [] for component in wanted_components: i = 'ned'.find(component) if i != -1: outs_wanted.append(outs[i]) else: outs_wanted.append(None)
return tstart, outs_wanted
def add_seismogram( vp, vs, density, qp, qs, x, f, m6, quantity, deltat, out_offset, out_n, out_e, out_d, stf=None, want_far=True, want_intermediate=True, want_near=True, npad_levelling=40):
ns = [out.size for out in (out_n, out_e, out_d) if out is not None]
if not all(n == ns[0] for n in ns): raise AhfullgreenError('length of component arrays must be identical')
n = ns[0]
wanted_components = ''.join( (c if out is not None else '-') for (out, c) in zip((out_n, out_e, out_d), 'ned'))
tstart, temps = make_seismogram( vp, vs, density, qp, qs, x, f, m6, quantity, deltat, stf=stf, wanted_components=wanted_components, want_far=want_far, want_intermediate=want_intermediate, want_near=want_near, npad_levelling=npad_levelling, out_alignment=out_offset)
for i, out in enumerate((out_n, out_e, out_d)): if out is None: continue
temp = temps[i]
ntemp = temp.size
tmin = max(out_offset, tstart) tmax = min( out_offset + (n-1) * deltat, tstart + (ntemp-1) * deltat)
def ind(t, t0): return int(round((t-t0)/deltat))
out[ind(tmin, out_offset):ind(tmax, out_offset)+1] \ += temp[ind(tmin, tstart):ind(tmax, tstart)+1]
out[:ind(tmin, out_offset)] += 0. out[ind(tmax, out_offset)+1:] += temp[ind(tmax, tstart)]
class Impulse(object):
def __init__(self): pass
def t_cutoff(self): return None
def __call__(self, f): return num.ones(f.size, dtype=num.complex)
class Gauss(object): def __init__(self, tau): self._tau = tau
def t_cutoff(self): return self._tau * 2.
def __call__(self, f): omega = f * 2. * math.pi
return num.exp(-(omega**2 * self._tau**2 / 8.)) |