Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/apps/cake.py: 77%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Cake - travel-times and ray paths for 1D layered media. 

8''' 

9 

10import sys 

11import re 

12import numpy as num 

13# import logger 

14from pyrocko import cake, util, orthodrome 

15from pyrocko.plot import cake_plot as plot 

16from optparse import OptionParser, OptionGroup 

17import matplotlib.pyplot as plt 

18from pyrocko.plot import mpl_init, mpl_papersize, mpl_margins 

19 

20r2d = cake.r2d 

21 

22 

23class Anon(dict): 

24 

25 def __getattr__(self, x): 

26 return self[x] 

27 

28 def getn(self, *keys): 

29 return Anon([(k, self[k]) for k in keys]) 

30 

31 def gett(self, *keys): 

32 return tuple([self[k] for k in keys]) 

33 

34 

35def process_color(s, phase_colors): 

36 m = re.match(r'^([^{]+){([^}]*)}$', s) 

37 if not m: 

38 return s 

39 

40 sphase = m.group(1) 

41 color = m.group(2) 

42 phase_colors[sphase] = plot.str_to_mpl_color(color) 

43 

44 return sphase 

45 

46 

47def optparse( 

48 required=(), 

49 optional=(), 

50 args=sys.argv, 

51 usage='%prog [options]', 

52 descr=None): 

53 

54 want = required + optional 

55 

56 parser = OptionParser( 

57 prog='cake', 

58 usage=usage, 

59 description=descr.capitalize()+'.', 

60 add_help_option=False, 

61 formatter=util.BetterHelpFormatter()) 

62 

63 parser.add_option( 

64 '-h', '--help', action='help', help='Show help message and exit.') 

65 

66 if 'phases' in want: 

67 group = OptionGroup(parser, 'Phases', ''' 

68 

69Seismic phase arrivals may be either specified as traditional phase names (e.g. 

70P, S, PP, PcP, ...) or in Cake's own syntax which is more powerful. Use the 

71--classic option, for traditional phase names. Use the --phase option if you 

72want to define phases in Cake's syntax. 

73 

74''') 

75 group.add_option( 

76 '--phase', '--phases', dest='phases', action='append', 

77 default=[], metavar='PHASE1,PHASE2,...', 

78 help='''Comma separated list of seismic phases in Cake\'s syntax. 

79 

80The definition of a seismic propagation path in Cake's phase syntax is a string 

81consisting of an alternating sequence of "legs" and "knees". 

82 

83A "leg" represents seismic wave propagation without any conversions, 

84encountering only super-critical reflections. Legs are denoted by "P", "p", 

85"S", or "s". The capital letters are used when the take-off of the "leg" is 

86in downward direction, while the lower case letters indicate a take-off in 

87upward direction. 

88 

89A "knee" is an interaction with an interface. It can be a mode conversion, a 

90reflection, or propagation as a headwave or diffracted wave. 

91 

92 * conversion is simply denoted as: "(INTERFACE)" or "DEPTH" 

93 * upperside reflection: "v(INTERFACE)" or "vDEPTH" 

94 * underside reflection: "^(INTERFACE)" or "^DEPTH" 

95 * normal kind headwave or diffracted wave: "v_(INTERFACE)" or "v_DEPTH" 

96 

97The interface may be given by name or by depth: INTERFACE is the name of an 

98interface defined in the model, DEPTH is the depth of an interface in [km] (the 

99interface closest to that depth is chosen). If two legs appear consecutively 

100without an explicit "knee", surface interaction is assumed. 

101 

102The preferred standard interface names in cake are "conrad", "moho", "cmb" 

103(core-mantle boundary), and "cb" (inner core boundary). 

104 

105The phase definition may end with a backslash "\\", to indicate that the ray 

106should arrive at the receiver from above instead of from below. It is possible 

107to restrict the maximum and minimum depth of a "leg" by appending 

108"<(INTERFACE)" or "<DEPTH" or ">(INTERFACE)" or ">DEPTH" after the leg 

109character, respectively. 

110 

111When plotting rays or travel-time curves, the color can be set by appending 

112"{COLOR}" to the phase definition, where COLOR is the name of a color or an RGB 

113or RGBA color tuple in the format "R/G/B" or "R/G/B/A", respectively. The 

114values can be normalized to the range [0, 1] or to [0, 255]. The latter is only 

115assumed when any of the values given exceeds 1.0. 

116''') 

117 

118 group.add_option( 

119 '--classic', dest='classic_phases', action='append', 

120 default=[], metavar='PHASE1,PHASE2,...', 

121 help='''Comma separated list of seismic phases in classic 

122nomenclature. Run "cake list-phase-map" for a list of available 

123phase names. When plotting, color can be specified in the same way 

124as in --phases.''') 

125 

126 parser.add_option_group(group) 

127 

128 if 'model' in want: 

129 group = OptionGroup(parser, 'Model') 

130 group.add_option( 

131 '--model', dest='model_filename', metavar='(NAME or FILENAME)', 

132 help='Use builtin model named NAME or user model from file ' 

133 'FILENAME. By default, the "ak135-f-continental.m" model is ' 

134 'used. Run "cake list-models" for a list of builtin models.') 

135 

136 group.add_option( 

137 '--format', dest='model_format', metavar='FORMAT', 

138 choices=['nd', 'hyposat'], default='nd', 

139 help='Set model file format (available: nd, hyposat; default: ' 

140 'nd).') 

141 group.add_option( 

142 '--crust2loc', dest='crust2loc', metavar='LAT,LON', 

143 help='Set model from CRUST2.0 profile at location (LAT,LON).') 

144 group.add_option( 

145 '--crust2profile', dest='crust2profile', metavar='KEY', 

146 help='Set model from CRUST2.0 profile with given KEY.') 

147 

148 parser.add_option_group(group) 

149 

150 if any(x in want for x in ( 

151 'zstart', 'zstop', 'distances', 'sloc', 'rloc')): 

152 

153 group = OptionGroup(parser, 'Source-receiver geometry') 

154 if 'zstart' in want: 

155 group.add_option( 

156 '--sdepth', dest='sdepth', type='float', default=0.0, 

157 metavar='FLOAT', 

158 help='Source depth [km] (default: 0)') 

159 if 'zstop' in want: 

160 group.add_option( 

161 '--rdepth', dest='rdepth', type='float', default=0.0, 

162 metavar='FLOAT', 

163 help='Receiver depth [km] (default: 0)') 

164 if 'distances' in want: 

165 group.add_option( 

166 '--distances', dest='sdist', metavar='DISTANCES', 

167 help='Surface distances as "start:stop:n" or ' 

168 '"dist1,dist2,..." [km]') 

169 group.add_option( 

170 '--sloc', dest='sloc', metavar='LAT,LON', 

171 help='Source location (LAT,LON).') 

172 group.add_option( 

173 '--rloc', dest='rloc', metavar='LAT,LON', 

174 help='Receiver location (LAT,LON).') 

175 parser.add_option_group(group) 

176 

177 if 'material' in want: 

178 group = OptionGroup( 

179 parser, 'Material', 

180 'An isotropic elastic material may be specified by giving ' 

181 'a combination of some of the following options. ') 

182 group.add_option( 

183 '--vp', dest='vp', default=None, type='float', metavar='FLOAT', 

184 help='P-wave velocity [km/s]') 

185 group.add_option( 

186 '--vs', dest='vs', default=None, type='float', metavar='FLOAT', 

187 help='S-wave velocity [km/s]') 

188 group.add_option( 

189 '--rho', dest='rho', default=None, type='float', metavar='FLOAT', 

190 help='density [g/cm**3]') 

191 group.add_option( 

192 '--qp', dest='qp', default=None, type='float', metavar='FLOAT', 

193 help='P-wave attenuation Qp (default: 1456)') 

194 group.add_option( 

195 '--qs', dest='qs', default=None, type='float', metavar='FLOAT', 

196 help='S-wave attenuation Qs (default: 600)') 

197 group.add_option( 

198 '--poisson', dest='poisson', default=None, type='float', 

199 metavar='FLOAT', 

200 help='Poisson ratio') 

201 group.add_option( 

202 '--lambda', dest='lame_lambda', default=None, type='float', 

203 metavar='FLOAT', 

204 help='Lame parameter lambda [GPa]') 

205 group.add_option( 

206 '--mu', dest='lame_mu', default=None, type='float', 

207 metavar='FLOAT', 

208 help='Shear modulus [GPa]') 

209 group.add_option( 

210 '--qk', dest='qk', default=None, type='float', metavar='FLOAT', 

211 help='Bulk attenuation Qk') 

212 group.add_option( 

213 '--qmu', dest='qmu', default=None, type='float', metavar='FLOAT', 

214 help='Shear attenuation Qmu') 

215 parser.add_option_group(group) 

216 

217 if any(x in want for x in ( 

218 'vred', 'as_degrees', 'accuracy', 'slowness', 'interface', 

219 'aspect', 'shade_model')): 

220 

221 group = OptionGroup(parser, 'General') 

222 if 'vred' in want: 

223 group.add_option( 

224 '--vred', dest='vred', type='float', metavar='FLOAT', 

225 help='Velocity for time reduction in plot [km/s]') 

226 

227 if 'as_degrees' in want: 

228 group.add_option( 

229 '--degrees', dest='as_degrees', action='store_true', 

230 default=False, 

231 help='Distances are in [deg] instead of [km], velocities in ' 

232 '[deg/s] instead of [km/s], slownesses in [s/deg] ' 

233 'instead of [s/km].') 

234 

235 if 'accuracy' in want: 

236 group.add_option( 

237 '--accuracy', dest='accuracy', type='float', 

238 metavar='MAXIMUM_RELATIVE_RMS', default=0.002, 

239 help='Set accuracy for model simplification.') 

240 

241 if 'slowness' in want: 

242 group.add_option( 

243 '--slowness', dest='slowness', type='float', metavar='FLOAT', 

244 default=0.0, 

245 help='Select surface slowness [s/km] (default: 0)') 

246 

247 if 'interface' in want: 

248 group.add_option( 

249 '--interface', dest='interface', metavar='(NAME or DEPTH)', 

250 help='Name or depth [km] of interface to select') 

251 

252 if 'aspect' in want: 

253 group.add_option( 

254 '--aspect', dest='aspect', type='float', metavar='FLOAT', 

255 help='Aspect ratio for plot') 

256 

257 if 'shade_model' in want: 

258 group.add_option( 

259 '--no-shade-model', dest='shade_model', action='store_false', 

260 default=True, 

261 help='Suppress shading of earth model layers') 

262 

263 parser.add_option_group(group) 

264 

265 if any(x in want for x in ('output_format', 'save', 'size', 'show')): 

266 group = OptionGroup(parser, 'Output', 'Output specifications') 

267 if 'output_format' in want: 

268 group.add_option( 

269 '--output-format', dest='output_format', metavar='FORMAT', 

270 default='textual', 

271 choices=('textual', 'nd'), 

272 help='Set model output format (available: textual, nd, ' 

273 'default: textual)') 

274 if 'save' in want: 

275 group.add_option( 

276 '-s', '--save', dest='save', metavar='PATH', 

277 help='saves plot to .png (default) or other py-supported\ 

278 endings without showing, use --show or -u for showing', 

279 default='') 

280 if 'size' in want: 

281 group.add_option( 

282 '--size', dest='size', type='string', 

283 default='a4', 

284 help="gives size of returned plot, use 'a5' or 'a4'") 

285 if 'show' in want: 

286 group.add_option( 

287 '-u', '--show', dest='show', action='store_true', 

288 help='shows plot when saving (-u for unhide)') 

289 

290 parser.add_option_group(group) 

291 

292 if usage == 'cake help-options': 

293 parser.print_help() 

294 

295 (options, args) = parser.parse_args(args) 

296 

297 if len(args) != 2: 

298 parser.error( 

299 'Cake arguments should look like "--option" or "--option=...".') 

300 

301 d = {} 

302 as_degrees = False 

303 if 'as_degrees' in want: 

304 as_degrees = options.as_degrees 

305 d['as_degrees'] = as_degrees 

306 

307 if 'accuracy' in want: 

308 d['accuracy'] = options.accuracy 

309 

310 if 'output_format' in want: 

311 d['output_format'] = options.output_format 

312 

313 if 'save' in want: 

314 d['save'] = options.save 

315 

316 if 'size' in want: 

317 d['size'] = options.size 

318 

319 if 'show' in want: 

320 d['show'] = options.show 

321 

322 if 'aspect' in want: 

323 d['aspect'] = options.aspect 

324 

325 if 'shade_model' in want: 

326 d['shade_model'] = options.shade_model 

327 

328 if 'phases' in want: 

329 phases = [] 

330 phase_colors = {} 

331 try: 

332 for ss in options.phases: 

333 for s in ss.split(','): 

334 s = process_color(s, phase_colors) 

335 phases.append(cake.PhaseDef(s)) 

336 

337 for pp in options.classic_phases: 

338 for p in pp.split(','): 

339 p = process_color(p, phase_colors) 

340 phases.extend(cake.PhaseDef.classic(p)) 

341 

342 except (cake.PhaseDefParseError, cake.UnknownClassicPhase) as e: 

343 parser.error(e) 

344 

345 if not phases and 'phases' in required: 

346 s = process_color('P', phase_colors) 

347 phases.append(cake.PhaseDef(s)) 

348 

349 if phases: 

350 d['phase_colors'] = phase_colors 

351 d['phases'] = phases 

352 

353 if 'model' in want: 

354 if options.model_filename: 

355 d['model'] = cake.load_model( 

356 options.model_filename, options.model_format) 

357 

358 if options.crust2loc or options.crust2profile: 

359 if options.crust2loc: 

360 try: 

361 args = tuple( 

362 [float(x) for x in options.crust2loc.split(',')]) 

363 except Exception: 

364 parser.error( 

365 'format for --crust2loc option is ' 

366 '"LATITUDE,LONGITUDE"') 

367 elif options.crust2profile: 

368 args = (options.crust2profile.upper(),) 

369 else: 

370 assert False 

371 

372 if 'model' in d: 

373 d['model'] = d['model'].replaced_crust(args) 

374 else: 

375 from pyrocko import crust2x2 

376 profile = crust2x2.get_profile(*args) 

377 d['model'] = cake.LayeredModel.from_scanlines( 

378 cake.from_crust2x2_profile(profile)) 

379 

380 if 'vred' in want: 

381 d['vred'] = options.vred 

382 if d['vred'] is not None: 

383 if not as_degrees: 

384 d['vred'] *= r2d * cake.km / cake.earthradius 

385 

386 if 'distances' in want: 

387 distances = None 

388 if options.sdist: 

389 if options.sdist.find(':') != -1: 

390 ssn = options.sdist.split(':') 

391 if len(ssn) != 3: 

392 parser.error( 

393 'format for distances is ' 

394 '"min_distance:max_distance:n_distances"') 

395 

396 distances = num.linspace( 

397 float(ssn[0]), float(ssn[1]), int(ssn[2])) 

398 else: 

399 distances = num.array( 

400 list(map( 

401 float, options.sdist.split(','))), dtype=float) 

402 

403 if not as_degrees: 

404 distances *= r2d * cake.km / cake.earthradius 

405 

406 if options.sloc and options.rloc: 

407 try: 

408 slat, slon = tuple([float(x) for x in options.sloc.split(',')]) 

409 rlat, rlon = tuple([float(x) for x in options.rloc.split(',')]) 

410 except Exception: 

411 parser.error( 

412 'format for --sloc and --rloc options is ' 

413 '"LATITUDE,LONGITUDE"') 

414 

415 distance_sr = orthodrome.distance_accurate50m_numpy( 

416 slat, slon, rlat, rlon) 

417 

418 distance_sr *= r2d / cake.earthradius 

419 if distances is not None: 

420 distances = num.concatenate((distances, [distance_sr])) 

421 else: 

422 distances = num.array([distance_sr], dtype=float) 

423 

424 if distances is not None: 

425 d['distances'] = distances 

426 else: 

427 if 'distances' not in required: 

428 d['distances'] = None 

429 

430 if 'slowness' in want: 

431 d['slowness'] = options.slowness/cake.d2r 

432 if not as_degrees: 

433 d['slowness'] /= cake.km*cake.m2d 

434 

435 if 'interface' in want: 

436 if options.interface: 

437 try: 

438 d['interface'] = float(options.interface)*cake.km 

439 except ValueError: 

440 d['interface'] = options.interface 

441 

442 else: 

443 d['interface'] = None 

444 

445 if 'zstart' in want: 

446 d['zstart'] = options.sdepth*cake.km 

447 

448 if 'zstop' in want: 

449 d['zstop'] = options.rdepth*cake.km 

450 

451 if 'material' in want: 

452 md = {} 

453 userfactor = dict( 

454 vp=1000., vs=1000., rho=1000., qp=1., qs=1., qmu=1., qk=1., 

455 lame_lambda=1.0e9, lame_mu=1.0e9, poisson=1.) 

456 

457 for k in userfactor.keys(): 

458 if getattr(options, k) is not None: 

459 md[k] = getattr(options, k) * userfactor[k] 

460 

461 if not (bool('lame_lambda' in md) == bool('lame_mu' in md)): 

462 parser.error('lambda and mu must be specified both.') 

463 if 'lame_lambda' in md and 'lame_mu' in md: 

464 md['lame'] = md.pop('lame_lambda'), md.pop('lame_mu') 

465 

466 if md: 

467 try: 

468 d['material'] = cake.Material(**md) 

469 except cake.InvalidArguments as e: 

470 parser.error(str(e)) 

471 

472 for k in list(d.keys()): 

473 if k not in want: 

474 del d[k] 

475 

476 for k in required: 

477 if k not in d: 

478 if k == 'model': 

479 d['model'] = cake.load_model('ak135-f-continental.m') 

480 

481 elif k == 'distances': 

482 d['distances'] = num.linspace(10*cake.km, 100*cake.km, 10) \ 

483 / cake.earthradius * r2d 

484 

485 elif k == 'phases': 

486 d['phases'] = list(map(cake.PhaseDef, 'Pp')) 

487 

488 else: 

489 parser.error('missing %s' % k) 

490 

491 return Anon(d) 

492 

493 

494def my_simplify_model(mod, accuracy): 

495 mod_simple = mod.simplify(max_rel_error=accuracy) 

496 cake.write_nd_model_fh(mod_simple, sys.stdout) 

497 

498 

499def d2u(d): 

500 return dict((k.replace('-', '_'), v) for (k, v) in d.items()) 

501 

502 

503def mini_fmt(v, d=5, try_fmts='fe'): 

504 for fmt in try_fmts: 

505 for i in range(d, -1, -1): 

506 s = ('%%%i.%i%s' % (d, i, fmt)) % v 

507 if len(s) <= d and (v == 0.0 or float(s) != 0.0): 

508 return s 

509 

510 return s 

511 

512 

513def scatter_out_fmt(d, m, v): 

514 if v < 1e-8: 

515 return ' ', ' ' 

516 else: 

517 return '%s%s' % ('\\/'[d == cake.UP], 'SP'[m == cake.P]), \ 

518 mini_fmt(v*100, 5) 

519 

520 

521def scatter_in_fmt(d, m, dwant): 

522 if d == dwant: 

523 return '%s%s' % ('\\/'[d == cake.UP], 'SP'[m == cake.P]) 

524 else: 

525 return ' ' 

526 

527 

528def print_scatter(model, p=0.0, interface=None): 

529 if interface is not None: 

530 discontinuities = [model.discontinuity(interface)] 

531 else: 

532 discontinuities = model.discontinuities() 

533 

534 for discontinuity in discontinuities: 

535 print('%s (%g km)' % (discontinuity, discontinuity.z/cake.km)) 

536 print() 

537 cols = [] 

538 for in_direction in (cake.DOWN, cake.UP): 

539 for in_mode in (cake.P, cake.S): 

540 p_critical = discontinuity.critical_ps(in_mode)[ 

541 in_direction == cake.UP] 

542 

543 if p_critical is None or p >= p_critical: 

544 continue 

545 

546 vals = [] 

547 for out_direction in (cake.UP, cake.DOWN): 

548 for out_mode in (cake.S, cake.P): 

549 vals.append( 

550 (out_direction, out_mode, discontinuity.efficiency( 

551 in_direction, out_direction, 

552 in_mode, out_mode, p))) 

553 

554 if all(x[-1] == 0.0 for x in vals): 

555 continue 

556 

557 sout = [scatter_out_fmt(d, m, v) for (d, m, v) in vals] 

558 

559 sin1 = scatter_in_fmt(in_direction, in_mode, cake.DOWN) 

560 sin2 = scatter_in_fmt(in_direction, in_mode, cake.UP) 

561 line1 = '%s %5s %5s' % ( 

562 ' '*len(sin1), sout[0][1], sout[1][1]) 

563 line2 = '%s %-5s %-5s' % ( 

564 sin1, sout[0][0], sout[1][0]) 

565 line4 = '%s %-5s %-5s' % ( 

566 sin2, sout[2][0], sout[3][0]) 

567 line5 = '%s %5s %5s' % ( 

568 ' '*len(sin2), sout[2][1], sout[3][1]) 

569 line3 = '-' * len(line1) 

570 cols.append((line1, line2, line3, line4, line5)) 

571 

572 for cols in zip(*cols): 

573 print(' ' + ' '.join(cols)) 

574 

575 print() 

576 print() 

577 

578 

579def print_arrivals( 

580 model, distances=[], phases=cake.PhaseDef('P'), zstart=0.0, zstop=0.0, 

581 as_degrees=False): 

582 

583 headers = 'slow dist time take inci effi spre phase used'.split() 

584 space = (7, 5, 6, 4, 4, 4, 4, 17, 17) 

585 

586 if as_degrees: 

587 units = 's/deg deg s deg deg % %'.split() 

588 else: 

589 units = 's/km km s deg deg % %'.split() 

590 

591 hline = ' '.join(x.ljust(s) for (x, s) in zip(headers, space)) 

592 uline = ' '.join(('%s' % x).ljust(s) for (x, s) in zip(units, space)) 

593 

594 print(hline) 

595 print(uline) 

596 print('-' * len(hline)) 

597 

598 for ray in model.arrivals( 

599 distances=distances, phases=phases, zstart=zstart, zstop=zstop): 

600 

601 if as_degrees: 

602 sd = ray.x 

603 slow = ray.p/cake.r2d 

604 else: 

605 sd = ray.x*(cake.d2r*cake.earthradius/cake.km) 

606 slow = ray.p/(r2d*cake.d2m/cake.km) 

607 

608 su = '(%s)' % ray.path.used_phase(p=ray.p, eps=1.0).used_repr() 

609 

610 print(' '.join(tuple(mini_fmt(x, s).rjust(s) for (x, s) in zip( 

611 (slow, sd, ray.t, ray.takeoff_angle(), ray.incidence_angle(), 

612 100*ray.efficiency(), 100*ray.spreading()*ray.surface_sphere()), 

613 space)) + tuple( 

614 x.ljust(17) for x in (ray.path.phase.definition(), su)))) 

615 

616 

617def plot_init(size, save, show): 

618 fontsize = 9 

619 mpl_init() 

620 fig = plt.figure(figsize=mpl_papersize(size, 'landscape')) 

621 

622 labelpos = mpl_margins(fig, w=7., h=5., units=fontsize) 

623 axes = fig.add_subplot(1, 1, 1) 

624 labelpos(axes, 2., 1.5) 

625 

626 showplt = bool(show or not save) 

627 

628 return fig, axes, showplt 

629 

630 

631class CakeError(Exception): 

632 pass 

633 

634 

635def plot_end(save, fig, show=True): 

636 if save: 

637 try: 

638 fig.savefig(save) 

639 if show: 

640 plt.show() 

641 except OSError as e: 

642 raise CakeError(str(e)) 

643 except ValueError as e: 

644 raise CakeError(str(e)) 

645 

646 

647def main(args=None): 

648 ''' 

649 CLI entry point for Pyrocko's ``cake`` app. 

650 ''' 

651 

652 if args is None: 

653 args = sys.argv[1:] 

654 

655 subcommand_descriptions = { 

656 'print': 'get information on model/phase/material properties', 

657 'arrivals': 'print list of phase arrivals', 

658 'paths': 'print ray path details', 

659 'plot-xt': 'plot traveltime vs distance curves', 

660 'plot-xp': 'plot ray parameter vs distance curves', 

661 'plot-rays': 'plot ray propagation paths', 

662 'plot': 'plot combination of ray and traveltime curves', 

663 'plot-model': 'plot velocity model', 

664 'list-models': 'list builtin velocity models', 

665 'list-phase-map': 'show translation table for classic phase names', 

666 'simplify-model': 'create a simplified version of a layered model', 

667 'scatter': 'show details about scattering at model interfaces'} 

668 

669 usage = '''cake <subcommand> [options] 

670 

671Subcommands: 

672 

673 print %(print)s 

674 arrivals %(arrivals)s 

675 paths %(paths)s 

676 plot-xt %(plot_xt)s 

677 plot-xp %(plot_xp)s 

678 plot-rays %(plot_rays)s 

679 plot %(plot)s 

680 plot-model %(plot_model)s 

681 list-models %(list_models)s 

682 list-phase-map %(list_phase_map)s 

683 simplify-model %(simplify_model)s 

684 scatter %(scatter)s 

685 

686To get further help and a list of available options for any subcommand run: 

687 

688 cake <subcommand> --help 

689 

690'''.strip() % d2u(subcommand_descriptions) 

691 

692 usage_sub = 'cake %s [options]' 

693 if len(args) < 1: 

694 sys.exit('Usage: %s' % usage) 

695 

696 command = args[0] 

697 

698 args[0:0] = ['cake'] 

699 descr = subcommand_descriptions.get(command, None) 

700 subusage = usage_sub % command 

701 

702 if command == 'print': 

703 c = optparse( 

704 (), ('model', 'phases', 'material', 'output_format'), 

705 usage=subusage, descr=descr, args=args) 

706 

707 if 'model' in c: 

708 if c.output_format == 'textual': 

709 print(c.model) 

710 print() 

711 elif c.output_format == 'nd': 

712 cake.write_nd_model_fh(c.model, sys.stdout) 

713 

714 if 'phases' in c: 

715 for phase in c.phases: 

716 print(phase) 

717 print() 

718 

719 if 'material' in c: 

720 print(c.material.describe()) 

721 print() 

722 

723 elif command == 'arrivals': 

724 c = optparse( 

725 ('model', 'phases', 'distances'), 

726 ('zstart', 'zstop', 'as_degrees'), 

727 usage=subusage, descr=descr, args=args) 

728 

729 print_arrivals( 

730 c.model, 

731 **c.getn('zstart', 'zstop', 'phases', 'distances', 'as_degrees')) 

732 

733 elif command == 'paths': 

734 c = optparse( 

735 ('model', 'phases'), 

736 ('zstart', 'zstop', 'as_degrees'), 

737 usage=subusage, descr=descr, args=args) 

738 

739 mod = c.model 

740 for path in mod.gather_paths(**c.getn('phases', 'zstart', 'zstop')): 

741 print(path.describe(path.endgaps(c.zstart, c.zstop), c.as_degrees)) 

742 

743 elif command in ('plot-xt', 'plot-xp', 'plot-rays', 'plot'): 

744 if command in ('plot-xt', 'plot'): 

745 c = optparse( 

746 ('model', 'phases',), 

747 ('zstart', 'zstop', 'distances', 'as_degrees', 'vred', 

748 'phase_colors', 'save', 'size', 'show'), 

749 usage=subusage, descr=descr, args=args) 

750 else: 

751 c = optparse( 

752 ('model', 'phases'), 

753 ('zstart', 'zstop', 'distances', 'as_degrees', 'aspect', 

754 'shade_model', 'phase_colors', 'save', 'size', 'show'), 

755 usage=subusage, descr=descr, args=args) 

756 

757 mod = c.model 

758 paths = mod.gather_paths(**c.getn('phases', 'zstart', 'zstop')) 

759 

760 if c.distances is not None: 

761 arrivals = mod.arrivals( 

762 **c.getn('phases', 'zstart', 'zstop', 'distances')) 

763 else: 

764 arrivals = None