1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5import sys 

6import re 

7import numpy as num 

8# import logger 

9from pyrocko import cake, util, orthodrome 

10from pyrocko.plot import cake_plot as plot 

11from optparse import OptionParser, OptionGroup 

12import matplotlib.pyplot as plt 

13from pyrocko.plot import mpl_init, mpl_papersize, mpl_margins 

14 

15r2d = cake.r2d 

16 

17 

18class Anon(dict): 

19 

20 def __getattr__(self, x): 

21 return self[x] 

22 

23 def getn(self, *keys): 

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

25 

26 def gett(self, *keys): 

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

28 

29 

30def process_color(s, phase_colors): 

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

32 if not m: 

33 return s 

34 

35 sphase = m.group(1) 

36 color = m.group(2) 

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

38 

39 return sphase 

40 

41 

42def optparse( 

43 required=(), 

44 optional=(), 

45 args=sys.argv, 

46 usage='%prog [options]', 

47 descr=None): 

48 

49 want = required + optional 

50 

51 parser = OptionParser( 

52 prog='cake', 

53 usage=usage, 

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

55 add_help_option=False, 

56 formatter=util.BetterHelpFormatter()) 

57 

58 parser.add_option( 

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

60 

61 if 'phases' in want: 

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

63 

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

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

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

67want to define phases in Cake's syntax. 

68 

69''') 

70 group.add_option( 

71 '--phase', '--phases', dest='phases', action="append", 

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

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

74 

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

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

77 

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

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

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

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

82upward direction. 

83 

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

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

86 

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

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

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

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

91 

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

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

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

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

96 

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

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

99 

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

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

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

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

104character, respectively. 

105 

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

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

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

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

110assumed when any of the values given exceeds 1.0. 

111''') 

112 

113 group.add_option( 

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

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

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

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

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

119as in --phases.''') 

120 

121 parser.add_option_group(group) 

122 

123 if 'model' in want: 

124 group = OptionGroup(parser, 'Model') 

125 group.add_option( 

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

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

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

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

130 

131 group.add_option( 

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

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

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

135 'nd).') 

136 group.add_option( 

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

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

139 group.add_option( 

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

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

142 

143 parser.add_option_group(group) 

144 

145 if any(x in want for x in ( 

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

147 

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

149 if 'zstart' in want: 

150 group.add_option( 

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

152 metavar='FLOAT', 

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

154 if 'zstop' in want: 

155 group.add_option( 

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

157 metavar='FLOAT', 

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

159 if 'distances' in want: 

160 group.add_option( 

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

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

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

164 group.add_option( 

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

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

167 group.add_option( 

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

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

170 parser.add_option_group(group) 

171 

172 if 'material' in want: 

173 group = OptionGroup( 

174 parser, 'Material', 

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

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

177 group.add_option( 

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

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

180 group.add_option( 

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

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

183 group.add_option( 

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

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

186 group.add_option( 

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

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

189 group.add_option( 

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

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

192 group.add_option( 

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

194 metavar='FLOAT', 

195 help='Poisson ratio') 

196 group.add_option( 

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

198 metavar='FLOAT', 

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

200 group.add_option( 

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

202 metavar='FLOAT', 

203 help='Shear modulus [GPa]') 

204 group.add_option( 

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

206 help='Bulk attenuation Qk') 

207 group.add_option( 

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

209 help='Shear attenuation Qmu') 

210 parser.add_option_group(group) 

211 

212 if any(x in want for x in ( 

213 'vred', 'as_degrees', 'accuracy', 'slowness', 'interface', 

214 'aspect', 'shade_model')): 

215 

216 group = OptionGroup(parser, 'General') 

217 if 'vred' in want: 

218 group.add_option( 

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

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

221 

222 if 'as_degrees' in want: 

223 group.add_option( 

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

225 default=False, 

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

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

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

229 

230 if 'accuracy' in want: 

231 group.add_option( 

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

233 metavar='MAXIMUM_RELATIVE_RMS', default=0.002, 

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

235 

236 if 'slowness' in want: 

237 group.add_option( 

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

239 default=0.0, 

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

241 

242 if 'interface' in want: 

243 group.add_option( 

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

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

246 

247 if 'aspect' in want: 

248 group.add_option( 

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

250 help='Aspect ratio for plot') 

251 

252 if 'shade_model' in want: 

253 group.add_option( 

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

255 default=True, 

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

257 

258 parser.add_option_group(group) 

259 

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

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

262 if 'output_format' in want: 

263 group.add_option( 

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

265 default='textual', 

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

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

268 'default: textual)') 

269 if 'save' in want: 

270 group.add_option( 

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

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

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

274 default='') 

275 if 'size' in want: 

276 group.add_option( 

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

278 default='a4', 

279 help='gives size of returned plot, use \'a5\' or \'a4\'') 

280 if 'show' in want: 

281 group.add_option( 

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

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

284 

285 parser.add_option_group(group) 

286 

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

288 parser.print_help() 

289 

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

291 

292 if len(args) != 2: 

293 parser.error( 

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

295 

296 d = {} 

297 as_degrees = False 

298 if 'as_degrees' in want: 

299 as_degrees = options.as_degrees 

300 d['as_degrees'] = as_degrees 

301 

302 if 'accuracy' in want: 

303 d['accuracy'] = options.accuracy 

304 

305 if 'output_format' in want: 

306 d['output_format'] = options.output_format 

307 

308 if 'save' in want: 

309 d['save'] = options.save 

310 

311 if 'size' in want: 

312 d['size'] = options.size 

313 

314 if 'show' in want: 

315 d['show'] = options.show 

316 

317 if 'aspect' in want: 

318 d['aspect'] = options.aspect 

319 

320 if 'shade_model' in want: 

321 d['shade_model'] = options.shade_model 

322 

323 if 'phases' in want: 

324 phases = [] 

325 phase_colors = {} 

326 try: 

327 for ss in options.phases: 

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

329 s = process_color(s, phase_colors) 

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

331 

332 for pp in options.classic_phases: 

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

334 p = process_color(p, phase_colors) 

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

336 

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

338 parser.error(e) 

339 

340 if not phases and 'phases' in required: 

341 s = process_color('P', phase_colors) 

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

343 

344 if phases: 

345 d['phase_colors'] = phase_colors 

346 d['phases'] = phases 

347 

348 if 'model' in want: 

349 if options.model_filename: 

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

351 options.model_filename, options.model_format) 

352 

353 if options.crust2loc or options.crust2profile: 

354 if options.crust2loc: 

355 try: 

356 args = tuple( 

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

358 except Exception: 

359 parser.error( 

360 'format for --crust2loc option is ' 

361 '"LATITUDE,LONGITUDE"') 

362 elif options.crust2profile: 

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

364 else: 

365 assert False 

366 

367 if 'model' in d: 

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

369 else: 

370 from pyrocko import crust2x2 

371 profile = crust2x2.get_profile(*args) 

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

373 cake.from_crust2x2_profile(profile)) 

374 

375 if 'vred' in want: 

376 d['vred'] = options.vred 

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

378 if not as_degrees: 

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

380 

381 if 'distances' in want: 

382 distances = None 

383 if options.sdist: 

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

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

386 if len(ssn) != 3: 

387 parser.error( 

388 'format for distances is ' 

389 '"min_distance:max_distance:n_distances"') 

390 

391 distances = num.linspace( 

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

393 else: 

394 distances = num.array( 

395 list(map( 

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

397 

398 if not as_degrees: 

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

400 

401 if options.sloc and options.rloc: 

402 try: 

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

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

405 except Exception: 

406 parser.error( 

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

408 '"LATITUDE,LONGITUDE"') 

409 

410 distance_sr = orthodrome.distance_accurate50m_numpy( 

411 slat, slon, rlat, rlon) 

412 

413 distance_sr *= r2d / cake.earthradius 

414 if distances is not None: 

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

416 else: 

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

418 

419 if distances is not None: 

420 d['distances'] = distances 

421 else: 

422 if 'distances' not in required: 

423 d['distances'] = None 

424 

425 if 'slowness' in want: 

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

427 if not as_degrees: 

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

429 

430 if 'interface' in want: 

431 if options.interface: 

432 try: 

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

434 except ValueError: 

435 d['interface'] = options.interface 

436 

437 else: 

438 d['interface'] = None 

439 

440 if 'zstart' in want: 

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

442 

443 if 'zstop' in want: 

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

445 

446 if 'material' in want: 

447 md = {} 

448 userfactor = dict( 

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

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

451 

452 for k in userfactor.keys(): 

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

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

455 

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

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

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

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

460 

461 if md: 

462 try: 

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

464 except cake.InvalidArguments as e: 

465 parser.error(str(e)) 

466 

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

468 if k not in want: 

469 del d[k] 

470 

471 for k in required: 

472 if k not in d: 

473 if k == 'model': 

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

475 

476 elif k == 'distances': 

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

478 / cake.earthradius * r2d 

479 

480 elif k == 'phases': 

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

482 

483 else: 

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

485 

486 return Anon(d) 

487 

488 

489def my_simplify_model(mod, accuracy): 

490 mod_simple = mod.simplify(max_rel_error=accuracy) 

491 cake.write_nd_model_fh(mod_simple, sys.stdout) 

492 

493 

494def d2u(d): 

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

496 

497 

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

499 for fmt in try_fmts: 

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

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

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

503 return s 

504 

505 return s 

506 

507 

508def scatter_out_fmt(d, m, v): 

509 if v < 1e-8: 

510 return ' ', ' ' 

511 else: 

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

513 mini_fmt(v*100, 5) 

514 

515 

516def scatter_in_fmt(d, m, dwant): 

517 if d == dwant: 

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

519 else: 

520 return ' ' 

521 

522 

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

524 if interface is not None: 

525 discontinuities = [model.discontinuity(interface)] 

526 else: 

527 discontinuities = model.discontinuities() 

528 

529 for discontinuity in discontinuities: 

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

531 print() 

532 cols = [] 

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

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

535 p_critical = discontinuity.critical_ps(in_mode)[ 

536 in_direction == cake.UP] 

537 

538 if p_critical is None or p >= p_critical: 

539 continue 

540 

541 vals = [] 

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

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

544 vals.append( 

545 (out_direction, out_mode, discontinuity.efficiency( 

546 in_direction, out_direction, 

547 in_mode, out_mode, p))) 

548 

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

550 continue 

551 

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

553 

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

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

556 line1 = '%s %5s %5s' % ( 

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

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

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

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

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

562 line5 = '%s %5s %5s' % ( 

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

564 line3 = '-' * len(line1) 

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

566 

567 for cols in zip(*cols): 

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

569 

570 print() 

571 print() 

572 

573 

574def print_arrivals( 

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

576 as_degrees=False): 

577 

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

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

580 

581 if as_degrees: 

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

583 else: 

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

585 

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

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

588 

589 print(hline) 

590 print(uline) 

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

592 

593 for ray in model.arrivals( 

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

595 

596 if as_degrees: 

597 sd = ray.x 

598 slow = ray.p/cake.r2d 

599 else: 

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

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

602 

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

604 

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

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

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

608 space)) + tuple( 

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

610 

611 

612def plot_init(size, save, show): 

613 fontsize = 9 

614 mpl_init() 

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

616 

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

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

619 labelpos(axes, 2., 1.5) 

620 

621 showplt = bool(show or not save) 

622 

623 return fig, axes, showplt 

624 

625 

626class CakeError(Exception): 

627 pass 

628 

629 

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

631 if save: 

632 try: 

633 fig.savefig(save) 

634 if show: 

635 plt.show() 

636 except OSError as e: 

637 raise CakeError(str(e)) 

638 except ValueError as e: 

639 raise CakeError(str(e)) 

640 

641 

642def main(args=None): 

643 

644 if args is None: 

645 args = sys.argv[1:] 

646 

647 subcommand_descriptions = { 

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

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

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

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

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

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

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

655 'plot-model': 'plot velocity model', 

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

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

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

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

660 

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

662 

663Subcommands: 

664 

665 print %(print)s 

666 arrivals %(arrivals)s 

667 paths %(paths)s 

668 plot-xt %(plot_xt)s 

669 plot-xp %(plot_xp)s 

670 plot-rays %(plot_rays)s 

671 plot %(plot)s 

672 plot-model %(plot_model)s 

673 list-models %(list_models)s 

674 list-phase-map %(list_phase_map)s 

675 simplify-model %(simplify_model)s 

676 scatter %(scatter)s 

677 

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

679 

680 cake <subcommand> --help 

681 

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

683 

684 usage_sub = 'cake %s [options]' 

685 if len(args) < 1: 

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

687 

688 command = args[0] 

689 

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

691 descr = subcommand_descriptions.get(command, None) 

692 subusage = usage_sub % command 

693 

694 if command == 'print': 

695 c = optparse( 

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

697 usage=subusage, descr=descr, args=args) 

698 

699 if 'model' in c: 

700 if c.output_format == 'textual': 

701 print(c.model) 

702 print() 

703 elif c.output_format == 'nd': 

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

705 

706 if 'phases' in c: 

707 for phase in c.phases: 

708 print(phase) 

709 print() 

710 

711 if 'material' in c: 

712 print(c.material.describe()) 

713 print() 

714 

715 elif command == 'arrivals': 

716 c = optparse( 

717 ('model', 'phases', 'distances'), 

718 ('zstart', 'zstop', 'as_degrees'), 

719 usage=subusage, descr=descr, args=args) 

720 

721 print_arrivals( 

722 c.model, 

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

724 

725 elif command == 'paths': 

726 c = optparse( 

727 ('model', 'phases'), 

728 ('zstart', 'zstop', 'as_degrees'), 

729 usage=subusage, descr=descr, args=args) 

730 

731 mod = c.model 

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

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

734 

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

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

737 c = optparse( 

738 ('model', 'phases',), 

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

740 'phase_colors', 'save', 'size', 'show'), 

741 usage=subusage, descr=descr, args=args) 

742 else: 

743 c = optparse( 

744 ('model', 'phases'), 

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

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

747 usage=subusage, descr=descr, args=args) 

748 

749 mod = c.model 

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

751 

752 if c.distances is not None: 

753 arrivals = mod.arrivals( 

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

755 else: 

756 arrivals = None 

757 

758 fig, axes, showplt = plot_init(c.size, c.save, c.show) 

759 

760 if command == 'plot-xp': 

761 plot.my_xp_plot( 

762 paths, c.zstart, c.zstop, c.distances, 

763 c.as_degrees, show=showplt, phase_colors=c.phase_colors) 

764 

765 elif command == 'plot-xt': 

766 plot.my_xt_plot( 

767 paths, c.zstart, c.zstop, c.distances, c.as_degrees, 

768 vred=c.vred, show=showplt, 

769 phase_colors=c.phase_colors) 

770 

771 elif command == 'plot-rays': 

772 if c.as_degrees: 

773 plot.my_rays_plot_gcs( 

774 mod, paths, arrivals, c.zstart, c.zstop, c.distances, 

775 show=showplt, phase_colors=c.phase_colors) 

776 

777 else: 

778 plot.my_rays_plot( 

779 mod, paths, arrivals, c.zstart, c.zstop, c.distances, 

780 show=showplt, aspect=c.aspect, shade_model=c.shade_model, 

781 phase_colors=c.phase_colors) 

782 

783 elif command == 'plot': 

784 plot.my_combi_plot( 

785 mod, paths, arrivals, c.zstart, c.zstop, c.distances, 

786 c.as_degrees, show=showplt, vred=c.vred, 

787 phase_colors=c.phase_colors) 

788 

789 try: 

790 plot_end(save=c.save, fig=fig, show=c.show) 

791 except CakeError as e: 

792 exit('cake.py: %s' % str(e)) 

793 

794 elif command in ('plot-model',): 

795 c = optparse( 

796 ('model',), ('save', 'size', 'show'), 

797 usage=subusage, descr=descr, args=args) 

798 mod = c.model 

799 fig, axes, showplt = plot_init(c.size, c.save, c.show) 

800 plot.my_model_plot(mod, show=showplt) 

801 try: 

802 plot_end(save=c.save, fig=fig, show=c.show) 

803 except CakeError as e: 

804 exit('cake.py: %s' % str(e)) 

805 

806 elif command in ('simplify-model',): 

807 c = optparse(('model',), ('accuracy',), 

808 usage=subusage, descr=descr, args=args) 

809 my_simplify_model(c.model, c.accuracy) 

810 

811 elif command in ('list-models',): 

812 c = optparse((), (), usage=subusage, descr=descr, args=args) 

813 for x in cake.builtin_models(): 

814 print(x) 

815 

816 elif command in ('list-phase-map',): 

817 c = optparse((), (), usage=subusage, descr=descr, args=args) 

818 defs = cake.PhaseDef.classic_definitions() 

819 for k in sorted(defs.keys()): 

820 print('%-15s: %s' % (k, ', '.join(defs[k]))) 

821 

822 elif command in ('scatter',): 

823 c = optparse( 

824 ('model',), 

825 ('slowness', 'interface', 'as_degrees'), 

826 usage=subusage, descr=descr, args=args) 

827 

828 print_scatter(c.model, p=c.slowness, interface=c.interface) 

829 

830 elif command in ('help-options',): 

831 optparse( 

832 (), 

833 ('model', 'accuracy', 'slowness', 'interface', 'phases', 

834 'distances', 'zstart', 'zstop', 'distances', 'as_degrees', 

835 'material', 'vred', 'save'), 

836 usage='cake help-options', descr='list all available options', 

837 args=args) 

838 

839 elif command in ('--help', '-h', 'help'): 

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

841 

842 else: 

843 sys.exit('cake: no such subcommand: %s' % command) 

844 

845 

846if __name__ == '__main__': 

847 main()