Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/fomosto/qseis.py: 83%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6import numpy as num 

7import logging 

8import os 

9import shutil 

10import math 

11import copy 

12import signal 

13 

14from tempfile import mkdtemp 

15from subprocess import Popen, PIPE 

16from os.path import join as pjoin 

17 

18from pyrocko import gf 

19from pyrocko import trace, util, cake 

20from pyrocko.moment_tensor import MomentTensor, symmat6 

21from pyrocko.guts import Float, Int, Tuple, List, Complex, Bool, Object, String 

22 

23km = 1e3 

24 

25guts_prefix = 'pf' 

26 

27Timing = gf.meta.Timing 

28 

29logger = logging.getLogger('pyrocko.fomosto.qseis') 

30 

31# how to call the programs 

32program_bins = { 

33 'qseis.2006': 'fomosto_qseis2006', 

34 'qseis.2006a': 'fomosto_qseis2006a', 

35 'qseis.2006b': 'fomosto_qseis2006b', 

36} 

37 

38 

39def have_backend(): 

40 have_any = False 

41 for cmd in [[exe] for exe in program_bins.values()]: 

42 try: 

43 p = Popen(cmd, stdout=PIPE, stderr=PIPE, stdin=PIPE) 

44 (stdout, stderr) = p.communicate() 

45 have_any = True 

46 

47 except OSError: 

48 pass 

49 

50 return have_any 

51 

52 

53qseis_components = 'r t z v'.split() 

54qseis_greenf_names = ('ex', 'ss', 'ds', 'cl', 'fz', 'fh') 

55 

56 

57def nextpow2(i): 

58 return 2**int(math.ceil(math.log(i)/math.log(2.))) 

59 

60 

61def str_float_vals(vals): 

62 return ' '.join('%e' % val for val in vals) 

63 

64 

65def str_int_vals(vals): 

66 return ' '.join('%i' % val for val in vals) 

67 

68 

69def str_str_vals(vals): 

70 return ' '.join("'%s'" % val for val in vals) 

71 

72 

73def scl(cs): 

74 if not cs: 

75 return '\n#' 

76 

77 return '\n'+' '.join('(%e,%e)' % (c.real, c.imag) for c in cs) 

78 

79 

80def cake_model_to_config(mod): 

81 k = 1000. 

82 srows = [] 

83 ref_depth = 0. 

84 for i, row in enumerate(mod.to_scanlines()): 

85 depth, vp, vs, rho, qp, qs = row 

86 if i == 0: 

87 ref_depth = depth 

88 

89 row = [(depth-ref_depth)/k, vp/k, vs/k, rho/k, qp, qs] 

90 srows.append('%i %s' % (i+1, str_float_vals(row))) 

91 

92 return '\n'.join(srows), len(srows), ref_depth 

93 

94 

95class QSeisSourceMech(Object): 

96 pass 

97 

98 

99class QSeisSourceMechMT(QSeisSourceMech): 

100 mnn = Float.T(default=1.0) 

101 mee = Float.T(default=1.0) 

102 mdd = Float.T(default=1.0) 

103 mne = Float.T(default=0.0) 

104 mnd = Float.T(default=0.0) 

105 med = Float.T(default=0.0) 

106 

107 def string_for_config(self): 

108 return '1 %(mnn)15e %(mee)15e %(mdd)15e ' \ 

109 '%(mne)15e %(med)15e %(mnd)15e ' % self.__dict__ 

110 

111 

112class QSeisSourceMechSDR(QSeisSourceMech): 

113 m_iso = Float.T(default=0.0) 

114 m_clvd = Float.T(default=0.0) 

115 m_dc = Float.T(default=1.0e9) 

116 strike = Float.T(default=0.0) 

117 dip = Float.T(default=90.0) 

118 rake = Float.T(default=0.0) 

119 

120 def string_for_config(self): 

121 return '2 %(m_iso)15e %(m_clvd)15e %(m_dc)15e ' \ 

122 '%(strike)15e %(dip)15e %(rake)15e ' % self.__dict__ 

123 

124 

125class QSeisPropagationFilter(Object): 

126 min_depth = Float.T(default=0.0) 

127 max_depth = Float.T(default=0.0) 

128 filtered_phase = Int.T(default=0) 

129 

130 def string_for_config(self): 

131 return '%(min_depth)15e %(max_depth)15e ' \ 

132 '%(filtered_phase)i' % self.__dict__ 

133 

134 

135class QSeisPoleZeroFilter(Object): 

136 constant = Complex.T(default=(1+0j)) 

137 poles = List.T(Complex.T()) 

138 zeros = List.T(Complex.T()) 

139 

140 def string_for_config(self, version=None): 

141 if version in ('2006a', '2006b'): 

142 return '(%e,%e)\n%i%s\n%i%s' % ( 

143 self.constant.real, self.constant.imag, 

144 len(self.zeros), scl(self.zeros), 

145 len(self.poles), scl(self.poles)) 

146 elif version == '2006': 

147 return '%e\n%i%s\n%i%s' % ( 

148 abs(self.constant), 

149 len(self.zeros), scl(self.zeros), 

150 len(self.poles), scl(self.poles)) 

151 

152 

153class QSeisConfig(Object): 

154 

155 qseis_version = String.T(default='2006') 

156 time_region = Tuple.T(2, Timing.T(), default=( 

157 Timing('-10'), Timing('+890'))) 

158 

159 cut = Tuple.T(2, Timing.T(), optional=True) 

160 fade = Tuple.T(4, Timing.T(), optional=True) 

161 relevel_with_fade_in = Bool.T(default=False) 

162 

163 sw_algorithm = Int.T(default=0) 

164 slowness_window = Tuple.T(4, Float.T(default=0.0)) 

165 wavenumber_sampling = Float.T(default=2.5) 

166 aliasing_suppression_factor = Float.T(default=0.1) 

167 source_disk_radius = Float.T(optional=True) 

168 

169 filter_surface_effects = Int.T(default=0) 

170 filter_shallow_paths = Int.T(default=0) 

171 filter_shallow_paths_depth = Float.T(default=0.0) 

172 propagation_filters = List.T(QSeisPropagationFilter.T()) 

173 receiver_filter = QSeisPoleZeroFilter.T(optional=True) 

174 

175 sw_flat_earth_transform = Int.T(default=0) 

176 

177 gradient_resolution_vp = Float.T(default=0.0) 

178 gradient_resolution_vs = Float.T(default=0.0) 

179 gradient_resolution_density = Float.T(default=0.0) 

180 

181 wavelet_duration_samples = Float.T(default=0.0) 

182 wavelet_type = Int.T(default=2) 

183 user_wavelet_samples = List.T(Float.T()) 

184 

185 def items(self): 

186 return dict(self.T.inamevals(self)) 

187 

188 

189class QSeisConfigFull(QSeisConfig): 

190 

191 time_start = Float.T(default=0.0) 

192 time_reduction_velocity = Float.T(default=0.0) 

193 time_window = Float.T(default=900.0) 

194 

195 source_depth = Float.T(default=10.0) 

196 source_mech = QSeisSourceMech.T( 

197 optional=True, 

198 default=QSeisSourceMechMT.D()) 

199 

200 receiver_depth = Float.T(default=0.0) 

201 receiver_distances = List.T(Float.T()) 

202 nsamples = Int.T(default=256) 

203 

204 gf_sw_source_types = Tuple.T(6, Int.T(), default=(1, 1, 1, 1, 0, 0)) 

205 

206 gf_filenames = Tuple.T(6, String.T(), default=qseis_greenf_names) 

207 

208 seismogram_filename = String.T(default='seis') 

209 

210 receiver_azimuths = List.T(Float.T()) 

211 

212 earthmodel_1d = gf.meta.Earthmodel1D.T(optional=True) 

213 earthmodel_receiver_1d = gf.meta.Earthmodel1D.T(optional=True) 

214 

215 @staticmethod 

216 def example(): 

217 conf = QSeisConfigFull() 

218 conf.receiver_distances = [2000.] 

219 conf.receiver_azimuths = [0.] 

220 conf.time_start = -10.0 

221 conf.time_reduction_velocity = 15.0 

222 conf.earthmodel_1d = cake.load_model().extract(depth_max='cmb') 

223 conf.earthmodel_receiver_1d = None 

224 conf.sw_flat_earth_transform = 1 

225 return conf 

226 

227 def get_output_filenames(self, rundir): 

228 return [pjoin(rundir, self.seismogram_filename+'.t'+c) 

229 for c in qseis_components] 

230 

231 def get_output_filenames_gf(self, rundir): 

232 return [pjoin(rundir, fn+'.t'+c) 

233 for fn in self.gf_filenames for c in qseis_components] 

234 

235 def string_for_config(self): 

236 

237 def aggregate(xx): 

238 return len(xx), '\n'.join( 

239 [''] + [x.string_for_config() for x in xx]) 

240 

241 assert len(self.receiver_distances) > 0 

242 assert len(self.receiver_distances) == len(self.receiver_azimuths) 

243 assert self.earthmodel_1d is not None 

244 

245 d = self.__dict__.copy() 

246 

247 # fixing these switches here to reduce the amount of wrapper code 

248 d['sw_distance_unit'] = 1 # always give distances in [km] 

249 d['sw_t_reduce'] = 1 # time reduction always as velocity [km/s] 

250 d['sw_equidistant'] = 0 # always give all distances and azimuths 

251 d['sw_irregular_azimuths'] = 1 

252 

253 d['n_distances'] = len(self.receiver_distances) 

254 d['str_distances'] = str_float_vals(self.receiver_distances) 

255 d['str_azimuths'] = str_float_vals(self.receiver_azimuths) 

256 

257 model_str, nlines, ref_depth = cake_model_to_config(self.earthmodel_1d) 

258 d['n_model_lines'] = nlines 

259 d['model_lines'] = model_str 

260 

261 if self.earthmodel_receiver_1d: 

262 model_str, nlines, ref_depth2 = cake_model_to_config( 

263 self.earthmodel_receiver_1d) 

264 assert ref_depth == ref_depth2 

265 else: 

266 model_str = '# no receiver side model' 

267 nlines = 0 

268 

269 d['n_model_receiver_lines'] = nlines 

270 d['model_receiver_lines'] = model_str 

271 

272 d['str_slowness_window'] = str_float_vals(self.slowness_window) 

273 

274 if self.qseis_version == '2006b': 

275 sdr = self.source_disk_radius 

276 d['str_source_disk_radius'] \ 

277 = '\n %e |dble: source_radius;' % ( 

278 sdr if sdr is not None else -1.0) 

279 else: 

280 if self.source_disk_radius is not None: 

281 raise QSeisError( 

282 'This version of QSEIS does not support the ' 

283 '`source_disk_radius` parameter.') 

284 

285 d['str_source_disk_radius'] = '' 

286 

287 if self.propagation_filters and ref_depth != 0.0: 

288 raise QSeisError( 

289 'Earth model must start with zero depth if ' 

290 'propagation_filters are set.') 

291 

292 d['n_depth_ranges'], d['str_depth_ranges'] = \ 

293 aggregate(self.propagation_filters) 

294 

295 if self.wavelet_type == 0: # user wavelet 

296 d['str_w_samples'] = '\n' \ 

297 + '%i\n' % len(self.user_wavelet_samples) \ 

298 + str_float_vals(self.user_wavelet_samples) 

299 else: 

300 d['str_w_samples'] = '' 

301 

302 if self.receiver_filter: 

303 d['str_receiver_filter'] = self.receiver_filter.string_for_config( 

304 self.qseis_version) 

305 else: 

306 if self.qseis_version in ('2006a', '2006b'): 

307 d['str_receiver_filter'] = '(1.0,0.0)\n0\n#\n0' 

308 else: 

309 d['str_receiver_filter'] = '1.0\n0\n#\n0' 

310 

311 d['str_gf_sw_source_types'] = str_int_vals(self.gf_sw_source_types) 

312 d['str_gf_filenames'] = str_str_vals(self.gf_filenames) 

313 

314 if self.source_mech: 

315 d['str_source'] = "%s '%s'" % ( 

316 self.source_mech.string_for_config(), 

317 self.seismogram_filename) 

318 else: 

319 d['str_source'] = '0' 

320 

321 d['source_depth_rel'] = d['source_depth'] - ref_depth / km 

322 d['receiver_depth_rel'] = d['receiver_depth'] - ref_depth / km 

323 

324 template = '''# autogenerated QSEIS input by qseis.py 

325# 

326# This is the input file of FORTRAN77 program "qseis06" for calculation of 

327# synthetic seismograms based on a layered halfspace earth model. 

328# 

329# by 

330# Rongjiang Wang <wang@gfz-potsdam.de> 

331# GeoForschungsZentrum Potsdam 

332# Telegrafenberg, D-14473 Potsdam, Germany 

333# 

334# Last modified: Potsdam, Nov., 2006 

335# 

336# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 

337# If not specified, SI Unit System is used overall! 

338# 

339# Coordinate systems: 

340# cylindrical (z,r,t) with z = downward, 

341# r = from source outward, 

342# t = azmuth angle from north to east; 

343# cartesian (x,y,z) with x = north, 

344# y = east, 

345# z = downward; 

346# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 

347# 

348# SOURCE PARAMETERS 

349# ================= 

350# 1. source depth [km] 

351#------------------------------------------------------------------------------ 

352 %(source_depth_rel)e |dble: source_depth; 

353#------------------------------------------------------------------------------ 

354# 

355# RECEIVER PARAMETERS 

356# =================== 

357# 1. receiver depth [km] 

358# 2. switch for distance sampling role (1/0 = equidistant/irregular); switch 

359# for unit used (1/0 = km/deg) 

360# 3. number of distance samples 

361# 4. if equidistant, then start and end trace distance (> 0); else distance 

362# list (please order the receiver distances from small to large) 

363# 5. (reduced) time begin [sec] & length of time window [sec], number of time 

364# samples (<= 2*nfmax in qsglobal.h) 

365# 6. switch for unit of the following time reduction parameter: 1 = velocity 

366# [km/sec], 0 = slowness [sec/deg]; time reduction parameter 

367#------------------------------------------------------------------------------ 

368 %(receiver_depth_rel)e |dble: receiver_depth; 

369 %(sw_equidistant)i %(sw_distance_unit)i |int: sw_equidistant, sw_d_unit; 

370 %(n_distances)i |int: no_distances; 

371 %(str_distances)s |dble: d_1,d_n; or d_1,d_2, ...(no comments in between!); 

372 %(time_start)e %(time_window)e %(nsamples)i |dble: t_start,t_window; int: no_t_samples; 

373 %(sw_t_reduce)i %(time_reduction_velocity)e |int: sw_t_reduce; dble: t_reduce; 

374#------------------------------------------------------------------------------ 

375# 

376# WAVENUMBER INTEGRATION PARAMETERS 

377# ================================= 

378# 1. select slowness integration algorithm (0 = suggested for full wave-field 

379# modelling; 1 or 2 = suggested when using a slowness window with narrow 

380# taper range - a technique for suppressing space-domain aliasing); 

381# 2. 4 parameters for low and high slowness (Note 1) cut-offs [s/km] with 

382# tapering: 0 < slw1 < slw2 defining cosine taper at the lower end, and 0 < 

383# slw3 < slw4 defining the cosine taper at the higher end. default values 

384# will be used in case of inconsistent input of the cut-offs (possibly with 

385# much more computational effort); 

386# 3. parameter for sampling rate of the wavenumber integration (1 = sampled 

387# with the spatial Nyquist frequency, 2 = sampled with twice higher than 

388# the Nyquist, and so on: the larger this parameter, the smaller the space- 

389# domain aliasing effect, but also the more computation effort); 

390# 4. the factor for suppressing time domain aliasing (> 0 and <= 1) (Note 2). 

391#------------------------------------------------------------------------------ 

392 %(sw_algorithm)i |int: sw_algorithm; 

393 %(str_slowness_window)s |dble: slw(1-4); 

394 %(wavenumber_sampling)e |dble: sample_rate; 

395 %(aliasing_suppression_factor)e |dble: supp_factor;%(str_source_disk_radius)s 

396#------------------------------------------------------------------------------ 

397# 

398# OPTIONS FOR PARTIAL SOLUTIONS 

399# (only applied to the source-site structure) 

400# =========================================== 

401# 

402# 1. switch for filtering free surface effects (0 = with free surface, i.e., 

403# do not select this filter; 1 = without free surface; 2 = without free 

404# surface but with correction on amplitude and wave form. Note switch 2 

405# can only be used for receivers at the surface) 

406# 2. switch for filtering waves with a shallow penetration depth (concerning 

407# their whole trace from source to receiver), penetration depth limit [km] 

408# 

409# if this option is selected, waves whose travel path never exceeds the 

410# given depth limit will be filtered ("seismic nuting"). the condition for 

411# selecting this filter is that the given shallow path depth limit should 

412# be larger than both source and receiver depth. 

413# 

414# 3. number of depth ranges where the following selected up/down-sp2oing P or 

415# SV waves should be filtered 

416# 4. the 1. depth range: upper and lower depth [km], switch for filtering P 

417# or SV wave in this depth range: 

418# 

419# switch no: 1 2 3 4 other 

420# filtered phase: P(up) P(down) SV(up) SV(down) Error 

421# 

422# 5. the 2. ... 

423# 

424# The partial solution options are useful tools to increase the numerical 

425# significance of desired wave phases. Especially when the desired phases 

426# are smaller than the undesired phases, these options should be selected 

427# and carefully combined. 

428#------------------------------------------------------------------------------ 

429 %(filter_surface_effects)i |int: isurf; 

430 %(filter_shallow_paths)i %(filter_shallow_paths_depth)e |int: sw_path_filter; dble:shallow_depth_limit; 

431 %(n_depth_ranges)i %(str_depth_ranges)s 

432#------------------------------------------------------------------------------ 

433# 

434# SOURCE TIME FUNCTION (WAVELET) PARAMETERS (Note 3) 

435# ================================================== 

436# 1. wavelet duration [unit = time sample rather than sec!], that is about 

437# equal to the half-amplitude cut-off period of the wavelet (> 0. if <= 0, 

438# then default value = 2 time samples will be used), and switch for the 

439# wavelet form (0 = user's own wavelet; 1 = default wavelet: normalized 

440# square half-sinusoid for simulating a physical delta impulse; 2 = tapered 

441# Heaviside wavelet, i.e. integral of wavelet 1) 

442# 2. IF user's own wavelet is selected, then number of the wavelet time samples 

443# (<= 1024), and followed by 

444# 3. equidistant wavelet time samples 

445# 4 ...(continue) (! no comment lines allowed between the time sample list!) 

446# IF default, delete line 2, 3, 4 ... or comment them out! 

447#------------------------------------------------------------------------------ 

448 %(wavelet_duration_samples)e %(wavelet_type)i%(str_w_samples)s 

449#------------------------------------------------------------------------------ 

450# 

451# FILTER PARAMETERS OF RECEIVERS (SEISMOMETERS OR HYDROPHONES) 

452# ============================================================ 

453# 1. constant coefficient (normalization factor) 

454# 2. number of roots (<= nrootmax in qsglobal.h) 

455# 3. list of the root positions in the complex format (Re,Im). If no roots, 

456# comment out this line 

457# 4. number of poles (<= npolemax in qsglobal.h) 

458# 5. list of the pole positions in the complex format (Re,Im). If no poles, 

459# comment out this line 

460#------------------------------------------------------------------------------ 

461 %(str_receiver_filter)s 

462#------------------------------------------------------------------------------ 

463# 

464# OUTPUT FILES FOR GREEN'S FUNCTIONS (Note 4) 

465# =========================================== 

466# 1. selections of source types (yes/no = 1/0) 

467# 2. file names of Green's functions (please give the names without extensions, 

468# which will be appended by the program automatically: *.tz, *.tr, *.tt 

469# and *.tv are for the vertical, radial, tangential, and volume change (for 

470# hydrophones) components, respectively) 

471#------------------------------------------------------------------------------ 

472# explosion strike-slip dip-slip clvd single_f_v single_f_h 

473#------------------------------------------------------------------------------ 

474 %(str_gf_sw_source_types)s 

475 %(str_gf_filenames)s 

476#------------------------------------------------------------------------------ 

477# OUTPUT FILES FOR AN ARBITRARY POINT DISLOCATION SOURCE 

478# (for applications to earthquakes) 

479# ====================================================== 

480# 1. selection (0 = not selected; 1 or 2 = selected), if (selection = 1), then 

481# the 6 moment tensor elements [N*m]: Mxx, Myy, Mzz, Mxy, Myz, Mzx (x is 

482# northward, y is eastward and z is downard); else if (selection = 2), then 

483# Mis [N*m] = isotropic moment part = (MT+MN+MP)/3, Mcl = CLVD moment part 

484# = (2/3)(MT+MP-2*MN), Mdc = double-couple moment part = MT-MN, Strike [deg], 

485# Dip [deg] and Rake [deg]. 

486# 

487# Note: to use this option, the Green's functions above should be computed 

488# (selection = 1) if they do not exist already. 

489# 

490# north(x) 

491# / 

492# /\ strike 

493# *-----------------------> east(y) 

494# |\ \ 

495# |-\ \ 

496# | \ fault plane \ 

497# |90 \ \ 

498# |-dip\ \ 

499# | \ \ 

500# | \ \ 

501# downward(z) \-----------------------\\ 

502# 

503# 2. switch for azimuth distribution of the stations (0 = uniform azimuth, 

504# else = irregular azimuth angles) 

505# 3. list of the azimuth angles [deg] for all stations given above (if the 

506# uniform azimuth is selected, then only one azimuth angle is required) 

507# 

508#------------------------------------------------------------------------------ 

509# Mis Mcl Mdc Strike Dip Rake File 

510#------------------------------------------------------------------------------ 

511# 2 0.00 1.00 6.0E+19 120.0 30.0 25.0 'seis' 

512#------------------------------------------------------------------------------ 

513# Mxx Myy Mzz Mxy Myz Mzx File 

514#------------------------------------------------------------------------------ 

515%(str_source)s 

516%(sw_irregular_azimuths)i 

517%(str_azimuths)s 

518#------------------------------------------------------------------------------ 

519# 

520# GLOBAL MODEL PARAMETERS (Note 5) 

521# ================================ 

522# 1. switch for flat-earth-transform 

523# 2. gradient resolution [%%] of vp, vs, and ro (density), if <= 0, then default 

524# values (depending on wave length at cut-off frequency) will be used 

525#------------------------------------------------------------------------------ 

526 %(sw_flat_earth_transform)i |int: sw_flat_earth_transform; 

527 %(gradient_resolution_vp)e %(gradient_resolution_vs)e %(gradient_resolution_density)e |dble: vp_res, vs_res, ro_res; 

528#------------------------------------------------------------------------------ 

529# 

530# LAYERED EARTH MODEL 

531# (SHALLOW SOURCE + UNIFORM DEEP SOURCE/RECEIVER STRUCTURE) 

532# ========================================================= 

533# 1. number of data lines of the layered model (source site) 

534#------------------------------------------------------------------------------ 

535 %(n_model_lines)i |int: no_model_lines; 

536#------------------------------------------------------------------------------ 

537# 

538# MULTILAYERED MODEL PARAMETERS (source site) 

539# =========================================== 

540# no depth[km] vp[km/s] vs[km/s] ro[g/cm^3] qp qs 

541#------------------------------------------------------------------------------ 

542%(model_lines)s 

543#------------------------------------------------------------------------------ 

544# 

545# LAYERED EARTH MODEL 

546# (ONLY THE SHALLOW RECEIVER STRUCTURE) 

547# ===================================== 

548# 1. number of data lines of the layered model 

549# 

550# Note: if the number = 0, then the receiver site is the same as the 

551# source site, else different receiver-site structure is considered. 

552# please be sure that the lowest interface of the receiver-site 

553# structure given given below can be found within the source-site 

554# structure, too. 

555# 

556#------------------------------------------------------------------------------ 

557 %(n_model_receiver_lines)i |int: no_model_lines; 

558#------------------------------------------------------------------------------ 

559# 

560# MULTILAYERED MODEL PARAMETERS (shallow receiver-site structure) 

561# =============================================================== 

562# no depth[km] vp[km/s] vs[km/s] ro[g/cm^3] qp qs 

563#------------------------------------------------------------------------------ 

564%(model_receiver_lines)s 

565#---------------------------------end of all inputs---------------------------- 

566 

567 

568Note 1: 

569 

570The slowness is defined by inverse value of apparent wave velocity = sin(i)/v 

571with i = incident angle and v = true wave velocity. 

572 

573Note 2: 

574 

575The suppression of the time domain aliasing is achieved by using the complex 

576frequency technique. The suppression factor should be a value between 0 and 1. 

577If this factor is set to 0.1, for example, the aliasing phase at the reduced 

578time begin is suppressed to 10%%. 

579 

580Note 3: 

581 

582The default basic wavelet function (option 1) is (2/tau)*sin^2(pi*t/tau), 

583for 0 < t < tau, simulating physical delta impuls. Its half-amplitude cut-off 

584frequency is 1/tau. To avoid high-frequency noise, tau should not be smaller 

585than 4-5 time samples. 

586 

587Note 4: 

588 

589 Double-Couple m11/ m22/ m33/ m12/ m23/ m31 Azimuth_Factor_(tz,tr,tv)/(tt) 

590 ============================================================================ 

591 explosion 1.0/ 1.0/ 1.0/ -- / -- / -- 1.0 / 0.0 

592 strike-slip -- / -- / -- / 1.0/ -- / -- sin(2*azi) / cos(2*azi) 

593 1.0/-1.0/ -- / -- / -- / -- cos(2*azi) / -sin(2*azi) 

594 dip-slip -- / -- / -- / -- / -- / 1.0 cos(azi) / sin(azi) 

595 -- / -- / -- / -- / 1.0/ -- sin(azi) / -cos(azi) 

596 clvd -0.5/-0.5/ 1.0/ -- / -- / -- 1.0 / 0.0 

597 ============================================================================ 

598 Single-Force fx / fy / fz Azimuth_Factor_(tz,tr,tv)/(tt) 

599 ============================================================================ 

600 fz -- / -- / 1.0 1.0 / 0.0 

601 fx 1.0/ -- / -- cos(azi) / sin(azi) 

602 fy -- / 1.0/ -- sin(azi) / -cos(azi) 

603 ============================================================================ 

604 

605Note 5: 

606 

607Layers with a constant gradient will be discretized with a number of homogeneous 

608sublayers. The gradient resolutions are then used to determine the maximum 

609allowed thickness of the sublayers. If the resolutions of Vp, Vs and Rho 

610(density) require different thicknesses, the smallest is first chosen. If this 

611is even smaller than 1%% of the characteristic wavelength, then the latter is 

612taken finally for the sublayer thickness. 

613''' # noqa 

614 

615 return (template % d).encode('ascii') 

616 

617 

618class QSeisError(gf.store.StoreError): 

619 pass 

620 

621 

622class Interrupted(gf.store.StoreError): 

623 def __str__(self): 

624 return 'Interrupted.' 

625 

626 

627class QSeisRunner(object): 

628 

629 def __init__(self, tmp=None, keep_tmp=False): 

630 self.tempdir = mkdtemp(prefix='qseisrun-', dir=tmp) 

631 self.keep_tmp = keep_tmp 

632 self.config = None 

633 

634 def run(self, config): 

635 self.config = config 

636 

637 input_fn = pjoin(self.tempdir, 'input') 

638 

639 with open(input_fn, 'wb') as f: 

640 input_str = config.string_for_config() 

641 logger.debug('===== begin qseis input =====\n' 

642 '%s===== end qseis input =====' % input_str.decode()) 

643 f.write(input_str) 

644 

645 program = program_bins['qseis.%s' % config.qseis_version] 

646 

647 old_wd = os.getcwd() 

648 

649 os.chdir(self.tempdir) 

650 

651 interrupted = [] 

652 

653 def signal_handler(signum, frame): 

654 os.kill(proc.pid, signal.SIGTERM) 

655 interrupted.append(True) 

656 

657 original = signal.signal(signal.SIGINT, signal_handler) 

658 try: 

659 try: 

660 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE) 

661 

662 except OSError: 

663 os.chdir(old_wd) 

664 raise QSeisError( 

665 '''could not start qseis executable: "%s" 

666Available fomosto backends and download links to the modelling codes are listed 

667on 

668 

669 https://pyrocko.org/docs/current/apps/fomosto/backends.html 

670 

671''' % program) 

672 

673 (output_str, error_str) = proc.communicate(b'input\n') 

674 

675 finally: 

676 signal.signal(signal.SIGINT, original) 

677 

678 if interrupted: 

679 raise KeyboardInterrupt() 

680 

681 logger.debug('===== begin qseis output =====\n' 

682 '%s===== end qseis output =====' % output_str.decode()) 

683 

684 errmess = [] 

685 if proc.returncode != 0: 

686 errmess.append( 

687 'qseis had a non-zero exit state: %i' % proc.returncode) 

688 

689 if error_str: 

690 logger.warning( 

691 'qseis emitted something via stderr:\n\n%s' 

692 % error_str.decode()) 

693 

694 # errmess.append('qseis emitted something via stderr') 

695 

696 if output_str.lower().find(b'error') != -1: 

697 errmess.append("the string 'error' appeared in qseis output") 

698 

699 if errmess: 

700 self.keep_tmp = True 

701 

702 os.chdir(old_wd) 

703 raise QSeisError(''' 

704===== begin qseis input ===== 

705%s===== end qseis input ===== 

706===== begin qseis output ===== 

707%s===== end qseis output ===== 

708===== begin qseis error ===== 

709%s===== end qseis error ===== 

710%s 

711qseis has been invoked as "%s" 

712in the directory %s'''.lstrip() % ( 

713 input_str.decode(), output_str.decode(), error_str.decode(), 

714 '\n'.join(errmess), program, self.tempdir)) 

715 

716 self.qseis_output = output_str 

717 self.qseis_error = error_str 

718 

719 os.chdir(old_wd) 

720 

721 def get_traces(self, which='seis'): 

722 

723 if which == 'seis': 

724 fns = self.config.get_output_filenames(self.tempdir) 

725 components = qseis_components 

726 

727 elif which == 'gf': 

728 fns = self.config.get_output_filenames_gf(self.tempdir) 

729 components = [ 

730 fn+'.t'+c 

731 for fn in self.config.gf_filenames for c in qseis_components] 

732 else: 

733 raise Exception( 

734 'get_traces: which argument should be "seis" or "gf"') 

735 

736 traces = [] 

737 distances = self.config.receiver_distances 

738 azimuths = self.config.receiver_azimuths 

739 for comp, fn in zip(components, fns): 

740 if not os.path.exists(fn): 

741 continue 

742 

743 data = num.loadtxt(fn, skiprows=1, dtype=float) 

744 nsamples, ntraces = data.shape 

745 ntraces -= 1 

746 vred = self.config.time_reduction_velocity 

747 deltat = (data[-1, 0] - data[0, 0])/(nsamples-1) 

748 

749 for ireceiver, distance, azimuth in zip( 

750 range(ntraces), distances, azimuths): 

751 

752 tmin = self.config.time_start 

753 if vred != 0.0: 

754 tmin += distance / vred 

755 

756 tmin += deltat 

757 tr = trace.Trace( 

758 '', '%04i' % ireceiver, '', comp, 

759 tmin=tmin, deltat=deltat, ydata=data[:, ireceiver+1], 

760 meta=dict( 

761 ireceiver=ireceiver, 

762 distance=distance*km, 

763 azimuth=azimuth)) 

764 

765 traces.append(tr) 

766 

767 return traces 

768 

769 def __del__(self): 

770 if self.tempdir: 

771 if not self.keep_tmp: 

772 shutil.rmtree(self.tempdir) 

773 self.tempdir = None 

774 else: 

775 logger.warning( 

776 'not removing temporary directory: %s' % self.tempdir) 

777 

778 

779def in_array_almost_equal(x, arr, eps): 

780 return num.any(num.abs(x - arr) < eps) 

781 

782 

783class QSeisGFBuilder(gf.builder.Builder): 

784 def __init__(self, store_dir, step, shared, block_size=None, tmp=None, 

785 force=False): 

786 

787 self.store = gf.store.Store(store_dir, 'w') 

788 

789 ndistances = 100 

790 if self.store.config.component_scheme == 'elastic10_fd': 

791 ndistances = 33 

792 

793 if block_size is None: 

794 block_size = (1, 1, ndistances) 

795 

796 if len(self.store.config.ns) == 2: 

797 block_size = block_size[1:] 

798 

799 gf.builder.Builder.__init__( 

800 self, self.store.config, step, block_size=block_size, force=force) 

801 

802 baseconf = self.store.get_extra('qseis') 

803 

804 conf = QSeisConfigFull(**baseconf.items()) 

805 conf.earthmodel_1d = self.store.config.earthmodel_1d 

806 conf.earthmodel_receiver_1d = self.store.config.earthmodel_receiver_1d 

807 

808 deltat = 1.0/self.gf_config.sample_rate 

809 

810 if 'time_window_min' not in shared: 

811 d = self.store.make_timing_params( 

812 conf.time_region[0], conf.time_region[1], 

813 force=force) 

814 

815 if self.store.config.component_scheme == 'elastic10_fd': 

816 # using reduction velocity is not possible here 

817 tmax = math.ceil(d['tmax'] / deltat) * deltat 

818 tmin = math.floor(d['tmin'] / deltat) * deltat 

819 shared['time_window_min'] = tmax - tmin 

820 shared['time_start'] = tmin 

821 shared['time_reduction_velocity'] = 0.0 

822 else: 

823 shared['time_window_min'] = d['tlenmax_vred'] 

824 shared['time_start'] = d['tmin_vred'] 

825 shared['time_reduction_velocity'] = d['vred'] / km 

826 

827 time_window_min = shared['time_window_min'] 

828 conf.time_start = shared['time_start'] 

829 

830 conf.time_reduction_velocity = shared['time_reduction_velocity'] 

831 

832 conf.nsamples = nextpow2(int(round(time_window_min / deltat)) + 1) 

833 conf.time_window = (conf.nsamples-1)*deltat 

834 

835 self.qseis_config = conf 

836 

837 self.tmp = tmp 

838 if self.tmp is not None: 

839 util.ensuredir(self.tmp) 

840 

841 def cleanup(self): 

842 self.store.close() 

843 

844 def work_block(self, index): 

845 if len(self.store.config.ns) == 2: 

846 (sz, firstx), (sz, lastx), (ns, nx) = \ 

847 self.get_block_extents(index) 

848 

849 rz = self.store.config.receiver_depth 

850 else: 

851 (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \ 

852 self.get_block_extents(index) 

853 

854 conf = copy.deepcopy(self.qseis_config) 

855 

856 logger.info('Starting block %i / %i' % 

857 (index+1, self.nblocks)) 

858 

859 conf.source_depth = float(sz/km) 

860 conf.receiver_depth = float(rz/km) 

861 

862 component_scheme = self.store.config.component_scheme 

863 

864 runner = QSeisRunner(tmp=self.tmp) 

865 

866 dx = self.gf_config.distance_delta 

867 

868 distances = num.linspace(firstx, firstx + (nx-1)*dx, nx) 

869 if component_scheme == 'elastic10_fd': 

870 

871 fd_distances = [ 

872 distances - self.gf_config.fd_distance_delta, 

873 distances, 

874 distances + self.gf_config.fd_distance_delta] 

875 

876 n_fd_distances = distances.size 

877 

878 distances = num.concatenate(fd_distances) 

879 

880 distances = distances.tolist() 

881 

882 # add global max distance, because qseis does some adjustments with 

883 # this value 

884 distances.append(self.gf_config.distance_max) 

885 

886 nreceivers = len(distances)