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
15r2d = cake.r2d
18class Anon(dict):
20 def __getattr__(self, x):
21 return self[x]
23 def getn(self, *keys):
24 return Anon([(k, self[k]) for k in keys])
26 def gett(self, *keys):
27 return tuple([self[k] for k in keys])
30def process_color(s, phase_colors):
31 m = re.match(r'^([^{]+){([^}]*)}$', s)
32 if not m:
33 return s
35 sphase = m.group(1)
36 color = m.group(2)
37 phase_colors[sphase] = plot.str_to_mpl_color(color)
39 return sphase
42def optparse(
43 required=(),
44 optional=(),
45 args=sys.argv,
46 usage='%prog [options]',
47 descr=None):
49 want = required + optional
51 parser = OptionParser(
52 prog='cake',
53 usage=usage,
54 description=descr.capitalize()+'.',
55 add_help_option=False,
56 formatter=util.BetterHelpFormatter())
58 parser.add_option(
59 '-h', '--help', action='help', help='Show help message and exit.')
61 if 'phases' in want:
62 group = OptionGroup(parser, 'Phases', '''
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.
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.
75The definition of a seismic propagation path in Cake's phase syntax is a string
76consisting of an alternating sequence of "legs" and "knees".
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.
84A "knee" is an interaction with an interface. It can be a mode conversion, a
85reflection, or propagation as a headwave or diffracted wave.
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"
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.
97The preferred standard interface names in cake are "conrad", "moho", "cmb"
98(core-mantle boundary), and "cb" (inner core boundary).
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.
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''')
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.''')
121 parser.add_option_group(group)
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.')
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.')
143 parser.add_option_group(group)
145 if any(x in want for x in (
146 'zstart', 'zstop', 'distances', 'sloc', 'rloc')):
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)
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)
212 if any(x in want for x in (
213 'vred', 'as_degrees', 'accuracy', 'slowness', 'interface',
214 'aspect', 'shade_model')):
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]')
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].')
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.')
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)')
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')
247 if 'aspect' in want:
248 group.add_option(
249 '--aspect', dest='aspect', type='float', metavar='FLOAT',
250 help='Aspect ratio for plot')
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')
258 parser.add_option_group(group)
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)')
285 parser.add_option_group(group)
287 if usage == 'cake help-options':
288 parser.print_help()
290 (options, args) = parser.parse_args(args)
292 if len(args) != 2:
293 parser.error(
294 'Cake arguments should look like "--option" or "--option=...".')
296 d = {}
297 as_degrees = False
298 if 'as_degrees' in want:
299 as_degrees = options.as_degrees
300 d['as_degrees'] = as_degrees
302 if 'accuracy' in want:
303 d['accuracy'] = options.accuracy
305 if 'output_format' in want:
306 d['output_format'] = options.output_format
308 if 'save' in want:
309 d['save'] = options.save
311 if 'size' in want:
312 d['size'] = options.size
314 if 'show' in want:
315 d['show'] = options.show
317 if 'aspect' in want:
318 d['aspect'] = options.aspect
320 if 'shade_model' in want:
321 d['shade_model'] = options.shade_model
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))
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))
337 except (cake.PhaseDefParseError, cake.UnknownClassicPhase) as e:
338 parser.error(e)
340 if not phases and 'phases' in required:
341 s = process_color('P', phase_colors)
342 phases.append(cake.PhaseDef(s))
344 if phases:
345 d['phase_colors'] = phase_colors
346 d['phases'] = phases
348 if 'model' in want:
349 if options.model_filename:
350 d['model'] = cake.load_model(
351 options.model_filename, options.model_format)
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
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))
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
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"')
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)
398 if not as_degrees:
399 distances *= r2d * cake.km / cake.earthradius
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"')
410 distance_sr = orthodrome.distance_accurate50m_numpy(
411 slat, slon, rlat, rlon)
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)
419 if distances is not None:
420 d['distances'] = distances
421 else:
422 if 'distances' not in required:
423 d['distances'] = None
425 if 'slowness' in want:
426 d['slowness'] = options.slowness/cake.d2r
427 if not as_degrees:
428 d['slowness'] /= cake.km*cake.m2d
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
437 else:
438 d['interface'] = None
440 if 'zstart' in want:
441 d['zstart'] = options.sdepth*cake.km
443 if 'zstop' in want:
444 d['zstop'] = options.rdepth*cake.km
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.)
452 for k in userfactor.keys():
453 if getattr(options, k) is not None:
454 md[k] = getattr(options, k) * userfactor[k]
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')
461 if md:
462 try:
463 d['material'] = cake.Material(**md)
464 except cake.InvalidArguments as e:
465 parser.error(str(e))
467 for k in list(d.keys()):
468 if k not in want:
469 del d[k]
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')
476 elif k == 'distances':
477 d['distances'] = num.linspace(10*cake.km, 100*cake.km, 10) \
478 / cake.earthradius * r2d
480 elif k == 'phases':
481 d['phases'] = list(map(cake.PhaseDef, 'Pp'))
483 else:
484 parser.error('missing %s' % k)
486 return Anon(d)
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)
494def d2u(d):
495 return dict((k.replace('-', '_'), v) for (k, v) in d.items())
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
505 return s
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)
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 ' '
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()
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]
538 if p_critical is None or p >= p_critical:
539 continue
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)))
549 if all(x[-1] == 0.0 for x in vals):
550 continue
552 sout = [scatter_out_fmt(d, m, v) for (d, m, v) in vals]
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))
567 for cols in zip(*cols):
568 print(' ' + ' '.join(cols))
570 print()
571 print()
574def print_arrivals(
575 model, distances=[], phases=cake.PhaseDef('P'), zstart=0.0, zstop=0.0,
576 as_degrees=False):
578 headers = 'slow dist time take inci effi spre phase used'.split()
579 space = (7, 5, 6, 4, 4, 4, 4, 17, 17)
581 if as_degrees:
582 units = 's/deg deg s deg deg % %'.split()
583 else:
584 units = 's/km km s deg deg % %'.split()
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))
589 print(hline)
590 print(uline)
591 print('-' * len(hline))
593 for ray in model.arrivals(
594 distances=distances, phases=phases, zstart=zstart, zstop=zstop):
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)
603 su = '(%s)' % ray.path.used_phase(p=ray.p, eps=1.0).used_repr()
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))))
612def plot_init(size, save, show):
613 fontsize = 9
614 mpl_init()
615 fig = plt.figure(figsize=mpl_papersize(size, 'landscape'))
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)
621 showplt = bool(show or not save)
623 return fig, axes, showplt
626class CakeError(Exception):
627 pass
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))
642def main(args=None):
644 if args is None:
645 args = sys.argv[1:]
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'}
661 usage = '''cake <subcommand> [options]
663Subcommands:
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
678To get further help and a list of available options for any subcommand run:
680 cake <subcommand> --help
682'''.strip() % d2u(subcommand_descriptions)
684 usage_sub = 'cake %s [options]'
685 if len(args) < 1:
686 sys.exit('Usage: %s' % usage)
688 command = args[0]
690 args[0:0] = ['cake']
691 descr = subcommand_descriptions.get(command, None)
692 subusage = usage_sub % command
694 if command == 'print':
695 c = optparse(
696 (), ('model', 'phases', 'material', 'output_format'),
697 usage=subusage, descr=descr, args=args)
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)
706 if 'phases' in c:
707 for phase in c.phases:
708 print(phase)
709 print()
711 if 'material' in c:
712 print(c.material.describe())
713 print()
715 elif command == 'arrivals':
716 c = optparse(
717 ('model', 'phases', 'distances'),
718 ('zstart', 'zstop', 'as_degrees'),
719 usage=subusage, descr=descr, args=args)
721 print_arrivals(
722 c.model,
723 **c.getn('zstart', 'zstop', 'phases', 'distances', 'as_degrees'))
725 elif command == 'paths':
726 c = optparse(
727 ('model', 'phases'),
728 ('zstart', 'zstop', 'as_degrees'),
729 usage=subusage, descr=descr, args=args)
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))
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)
749 mod = c.model
750 paths = mod.gather_paths(**c.getn('phases', 'zstart', 'zstop'))
752 if c.distances is not None:
753 arrivals = mod.arrivals(
754 **c.getn('phases', 'zstart', 'zstop', 'distances'))
755 else:
756 arrivals = None
758 fig, axes, showplt = plot_init(c.size, c.save, c.show)
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)
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)
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)
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)
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)
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))
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))
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)
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)
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])))
822 elif command in ('scatter',):
823 c = optparse(
824 ('model',),
825 ('slowness', 'interface', 'as_degrees'),
826 usage=subusage, descr=descr, args=args)
828 print_scatter(c.model, p=c.slowness, interface=c.interface)
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)
839 elif command in ('--help', '-h', 'help'):
840 sys.exit('Usage: %s' % usage)
842 else:
843 sys.exit('cake: no such subcommand: %s' % command)
846if __name__ == '__main__':
847 main()