1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

# 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.))