1#!/usr/bin/env python
2# http://pyrocko.org - GPLv3
3#
4# The Pyrocko Developers, 21st Century
5# ---|P------/S----------~Lg----------
6import sys
7import re
8import numpy as num
9# import logger
10from pyrocko import cake, util, orthodrome
11from pyrocko.plot import cake_plot as plot
12from optparse import OptionParser, OptionGroup
13import matplotlib.pyplot as plt
14from pyrocko.plot import mpl_init, mpl_papersize, mpl_margins
16r2d = cake.r2d
19class Anon(dict):
21 def __getattr__(self, x):
22 return self[x]
24 def getn(self, *keys):
25 return Anon([(k, self[k]) for k in keys])
27 def gett(self, *keys):
28 return tuple([self[k] for k in keys])
31def process_color(s, phase_colors):
32 m = re.match(r'^([^{]+){([^}]*)}$', s)
33 if not m:
34 return s
36 sphase = m.group(1)
37 color = m.group(2)
38 phase_colors[sphase] = plot.str_to_mpl_color(color)
40 return sphase
43def optparse(
44 required=(),
45 optional=(),
46 args=sys.argv,
47 usage='%prog [options]',
48 descr=None):
50 want = required + optional
52 parser = OptionParser(
53 prog='cake',
54 usage=usage,
55 description=descr.capitalize()+'.',
56 add_help_option=False,
57 formatter=util.BetterHelpFormatter())
59 parser.add_option(
60 '-h', '--help', action='help', help='Show help message and exit.')
62 if 'phases' in want:
63 group = OptionGroup(parser, 'Phases', '''
65Seismic phase arrivals may be either specified as traditional phase names (e.g.
66P, S, PP, PcP, ...) or in Cake's own syntax which is more powerful. Use the
67--classic option, for traditional phase names. Use the --phase option if you
68want to define phases in Cake's syntax.
70''')
71 group.add_option(
72 '--phase', '--phases', dest='phases', action='append',
73 default=[], metavar='PHASE1,PHASE2,...',
74 help='''Comma separated list of seismic phases in Cake\'s syntax.
76The definition of a seismic propagation path in Cake's phase syntax is a string
77consisting of an alternating sequence of "legs" and "knees".
79A "leg" represents seismic wave propagation without any conversions,
80encountering only super-critical reflections. Legs are denoted by "P", "p",
81"S", or "s". The capital letters are used when the take-off of the "leg" is
82in downward direction, while the lower case letters indicate a take-off in
83upward direction.
85A "knee" is an interaction with an interface. It can be a mode conversion, a
86reflection, or propagation as a headwave or diffracted wave.
88 * conversion is simply denoted as: "(INTERFACE)" or "DEPTH"
89 * upperside reflection: "v(INTERFACE)" or "vDEPTH"
90 * underside reflection: "^(INTERFACE)" or "^DEPTH"
91 * normal kind headwave or diffracted wave: "v_(INTERFACE)" or "v_DEPTH"
93The interface may be given by name or by depth: INTERFACE is the name of an
94interface defined in the model, DEPTH is the depth of an interface in [km] (the
95interface closest to that depth is chosen). If two legs appear consecutively
96without an explicit "knee", surface interaction is assumed.
98The preferred standard interface names in cake are "conrad", "moho", "cmb"
99(core-mantle boundary), and "cb" (inner core boundary).
101The phase definition may end with a backslash "\\", to indicate that the ray
102should arrive at the receiver from above instead of from below. It is possible
103to restrict the maximum and minimum depth of a "leg" by appending
104"<(INTERFACE)" or "<DEPTH" or ">(INTERFACE)" or ">DEPTH" after the leg
105character, respectively.
107When plotting rays or travel-time curves, the color can be set by appending
108"{COLOR}" to the phase definition, where COLOR is the name of a color or an RGB
109or RGBA color tuple in the format "R/G/B" or "R/G/B/A", respectively. The
110values can be normalized to the range [0, 1] or to [0, 255]. The latter is only
111assumed when any of the values given exceeds 1.0.
112''')
114 group.add_option(
115 '--classic', dest='classic_phases', action='append',
116 default=[], metavar='PHASE1,PHASE2,...',
117 help='''Comma separated list of seismic phases in classic
118nomenclature. Run "cake list-phase-map" for a list of available
119phase names. When plotting, color can be specified in the same way
120as in --phases.''')
122 parser.add_option_group(group)
124 if 'model' in want:
125 group = OptionGroup(parser, 'Model')
126 group.add_option(
127 '--model', dest='model_filename', metavar='(NAME or FILENAME)',
128 help='Use builtin model named NAME or user model from file '
129 'FILENAME. By default, the "ak135-f-continental.m" model is '
130 'used. Run "cake list-models" for a list of builtin models.')
132 group.add_option(
133 '--format', dest='model_format', metavar='FORMAT',
134 choices=['nd', 'hyposat'], default='nd',
135 help='Set model file format (available: nd, hyposat; default: '
136 'nd).')
137 group.add_option(
138 '--crust2loc', dest='crust2loc', metavar='LAT,LON',
139 help='Set model from CRUST2.0 profile at location (LAT,LON).')
140 group.add_option(
141 '--crust2profile', dest='crust2profile', metavar='KEY',
142 help='Set model from CRUST2.0 profile with given KEY.')
144 parser.add_option_group(group)
146 if any(x in want for x in (
147 'zstart', 'zstop', 'distances', 'sloc', 'rloc')):
149 group = OptionGroup(parser, 'Source-receiver geometry')
150 if 'zstart' in want:
151 group.add_option(
152 '--sdepth', dest='sdepth', type='float', default=0.0,
153 metavar='FLOAT',
154 help='Source depth [km] (default: 0)')
155 if 'zstop' in want:
156 group.add_option(
157 '--rdepth', dest='rdepth', type='float', default=0.0,
158 metavar='FLOAT',
159 help='Receiver depth [km] (default: 0)')
160 if 'distances' in want:
161 group.add_option(
162 '--distances', dest='sdist', metavar='DISTANCES',
163 help='Surface distances as "start:stop:n" or '
164 '"dist1,dist2,..." [km]')
165 group.add_option(
166 '--sloc', dest='sloc', metavar='LAT,LON',
167 help='Source location (LAT,LON).')
168 group.add_option(
169 '--rloc', dest='rloc', metavar='LAT,LON',
170 help='Receiver location (LAT,LON).')
171 parser.add_option_group(group)
173 if 'material' in want:
174 group = OptionGroup(
175 parser, 'Material',
176 'An isotropic elastic material may be specified by giving '
177 'a combination of some of the following options. ')
178 group.add_option(
179 '--vp', dest='vp', default=None, type='float', metavar='FLOAT',
180 help='P-wave velocity [km/s]')
181 group.add_option(
182 '--vs', dest='vs', default=None, type='float', metavar='FLOAT',
183 help='S-wave velocity [km/s]')
184 group.add_option(
185 '--rho', dest='rho', default=None, type='float', metavar='FLOAT',
186 help='density [g/cm**3]')
187 group.add_option(
188 '--qp', dest='qp', default=None, type='float', metavar='FLOAT',
189 help='P-wave attenuation Qp (default: 1456)')
190 group.add_option(
191 '--qs', dest='qs', default=None, type='float', metavar='FLOAT',
192 help='S-wave attenuation Qs (default: 600)')
193 group.add_option(
194 '--poisson', dest='poisson', default=None, type='float',
195 metavar='FLOAT',
196 help='Poisson ratio')
197 group.add_option(
198 '--lambda', dest='lame_lambda', default=None, type='float',
199 metavar='FLOAT',
200 help='Lame parameter lambda [GPa]')
201 group.add_option(
202 '--mu', dest='lame_mu', default=None, type='float',
203 metavar='FLOAT',
204 help='Shear modulus [GPa]')
205 group.add_option(
206 '--qk', dest='qk', default=None, type='float', metavar='FLOAT',
207 help='Bulk attenuation Qk')
208 group.add_option(
209 '--qmu', dest='qmu', default=None, type='float', metavar='FLOAT',
210 help='Shear attenuation Qmu')
211 parser.add_option_group(group)
213 if any(x in want for x in (
214 'vred', 'as_degrees', 'accuracy', 'slowness', 'interface',
215 'aspect', 'shade_model')):
217 group = OptionGroup(parser, 'General')
218 if 'vred' in want:
219 group.add_option(
220 '--vred', dest='vred', type='float', metavar='FLOAT',
221 help='Velocity for time reduction in plot [km/s]')
223 if 'as_degrees' in want:
224 group.add_option(
225 '--degrees', dest='as_degrees', action='store_true',
226 default=False,
227 help='Distances are in [deg] instead of [km], velocities in '
228 '[deg/s] instead of [km/s], slownesses in [s/deg] '
229 'instead of [s/km].')
231 if 'accuracy' in want:
232 group.add_option(
233 '--accuracy', dest='accuracy', type='float',
234 metavar='MAXIMUM_RELATIVE_RMS', default=0.002,
235 help='Set accuracy for model simplification.')
237 if 'slowness' in want:
238 group.add_option(
239 '--slowness', dest='slowness', type='float', metavar='FLOAT',
240 default=0.0,
241 help='Select surface slowness [s/km] (default: 0)')
243 if 'interface' in want:
244 group.add_option(
245 '--interface', dest='interface', metavar='(NAME or DEPTH)',
246 help='Name or depth [km] of interface to select')
248 if 'aspect' in want:
249 group.add_option(
250 '--aspect', dest='aspect', type='float', metavar='FLOAT',
251 help='Aspect ratio for plot')
253 if 'shade_model' in want:
254 group.add_option(
255 '--no-shade-model', dest='shade_model', action='store_false',
256 default=True,
257 help='Suppress shading of earth model layers')
259 parser.add_option_group(group)
261 if any(x in want for x in ('output_format', 'save', 'size', 'show')):
262 group = OptionGroup(parser, 'Output', 'Output specifications')
263 if 'output_format' in want:
264 group.add_option(
265 '--output-format', dest='output_format', metavar='FORMAT',
266 default='textual',
267 choices=('textual', 'nd'),
268 help='Set model output format (available: textual, nd, '
269 'default: textual)')
270 if 'save' in want:
271 group.add_option(
272 '-s', '--save', dest='save', metavar='PATH',
273 help='saves plot to .png (default) or other py-supported\
274 endings without showing, use --show or -u for showing',
275 default='')
276 if 'size' in want:
277 group.add_option(
278 '--size', dest='size', type='string',
279 default='a4',
280 help="gives size of returned plot, use 'a5' or 'a4'")
281 if 'show' in want:
282 group.add_option(
283 '-u', '--show', dest='show', action='store_true',
284 help='shows plot when saving (-u for unhide)')
286 parser.add_option_group(group)
288 if usage == 'cake help-options':
289 parser.print_help()
291 (options, args) = parser.parse_args(args)
293 if len(args) != 2:
294 parser.error(
295 'Cake arguments should look like "--option" or "--option=...".')
297 d = {}
298 as_degrees = False
299 if 'as_degrees' in want:
300 as_degrees = options.as_degrees
301 d['as_degrees'] = as_degrees
303 if 'accuracy' in want:
304 d['accuracy'] = options.accuracy
306 if 'output_format' in want:
307 d['output_format'] = options.output_format
309 if 'save' in want:
310 d['save'] = options.save
312 if 'size' in want:
313 d['size'] = options.size
315 if 'show' in want:
316 d['show'] = options.show
318 if 'aspect' in want:
319 d['aspect'] = options.aspect
321 if 'shade_model' in want:
322 d['shade_model'] = options.shade_model
324 if 'phases' in want:
325 phases = []
326 phase_colors = {}
327 try:
328 for ss in options.phases:
329 for s in ss.split(','):
330 s = process_color(s, phase_colors)
331 phases.append(cake.PhaseDef(s))
333 for pp in options.classic_phases:
334 for p in pp.split(','):
335 p = process_color(p, phase_colors)
336 phases.extend(cake.PhaseDef.classic(p))
338 except (cake.PhaseDefParseError, cake.UnknownClassicPhase) as e:
339 parser.error(e)
341 if not phases and 'phases' in required:
342 s = process_color('P', phase_colors)
343 phases.append(cake.PhaseDef(s))
345 if phases:
346 d['phase_colors'] = phase_colors
347 d['phases'] = phases
349 if 'model' in want:
350 if options.model_filename:
351 d['model'] = cake.load_model(
352 options.model_filename, options.model_format)
354 if options.crust2loc or options.crust2profile:
355 if options.crust2loc:
356 try:
357 args = tuple(
358 [float(x) for x in options.crust2loc.split(',')])
359 except Exception:
360 parser.error(
361 'format for --crust2loc option is '
362 '"LATITUDE,LONGITUDE"')
363 elif options.crust2profile:
364 args = (options.crust2profile.upper(),)
365 else:
366 assert False
368 if 'model' in d:
369 d['model'] = d['model'].replaced_crust(args)
370 else:
371 from pyrocko import crust2x2
372 profile = crust2x2.get_profile(*args)
373 d['model'] = cake.LayeredModel.from_scanlines(
374 cake.from_crust2x2_profile(profile))
376 if 'vred' in want:
377 d['vred'] = options.vred
378 if d['vred'] is not None:
379 if not as_degrees:
380 d['vred'] *= r2d * cake.km / cake.earthradius
382 if 'distances' in want:
383 distances = None
384 if options.sdist:
385 if options.sdist.find(':') != -1:
386 ssn = options.sdist.split(':')
387 if len(ssn) != 3:
388 parser.error(
389 'format for distances is '
390 '"min_distance:max_distance:n_distances"')
392 distances = num.linspace(
393 float(ssn[0]), float(ssn[1]), int(ssn[2]))
394 else:
395 distances = num.array(
396 list(map(
397 float, options.sdist.split(','))), dtype=float)
399 if not as_degrees:
400 distances *= r2d * cake.km / cake.earthradius
402 if options.sloc and options.rloc:
403 try:
404 slat, slon = tuple([float(x) for x in options.sloc.split(',')])
405 rlat, rlon = tuple([float(x) for x in options.rloc.split(',')])
406 except Exception:
407 parser.error(
408 'format for --sloc and --rloc options is '
409 '"LATITUDE,LONGITUDE"')
411 distance_sr = orthodrome.distance_accurate50m_numpy(
412 slat, slon, rlat, rlon)
414 distance_sr *= r2d / cake.earthradius
415 if distances is not None:
416 distances = num.concatenate((distances, [distance_sr]))
417 else:
418 distances = num.array([distance_sr], dtype=float)
420 if distances is not None:
421 d['distances'] = distances
422 else:
423 if 'distances' not in required:
424 d['distances'] = None
426 if 'slowness' in want:
427 d['slowness'] = options.slowness/cake.d2r
428 if not as_degrees:
429 d['slowness'] /= cake.km*cake.m2d
431 if 'interface' in want:
432 if options.interface:
433 try:
434 d['interface'] = float(options.interface)*cake.km
435 except ValueError:
436 d['interface'] = options.interface
438 else:
439 d['interface'] = None
441 if 'zstart' in want:
442 d['zstart'] = options.sdepth*cake.km
444 if 'zstop' in want:
445 d['zstop'] = options.rdepth*cake.km
447 if 'material' in want:
448 md = {}
449 userfactor = dict(
450 vp=1000., vs=1000., rho=1000., qp=1., qs=1., qmu=1., qk=1.,
451 lame_lambda=1.0e9, lame_mu=1.0e9, poisson=1.)
453 for k in userfactor.keys():
454 if getattr(options, k) is not None:
455 md[k] = getattr(options, k) * userfactor[k]
457 if not (bool('lame_lambda' in md) == bool('lame_mu' in md)):
458 parser.error('lambda and mu must be specified both.')
459 if 'lame_lambda' in md and 'lame_mu' in md:
460 md['lame'] = md.pop('lame_lambda'), md.pop('lame_mu')
462 if md:
463 try:
464 d['material'] = cake.Material(**md)
465 except cake.InvalidArguments as e:
466 parser.error(str(e))
468 for k in list(d.keys()):
469 if k not in want:
470 del d[k]
472 for k in required:
473 if k not in d:
474 if k == 'model':
475 d['model'] = cake.load_model('ak135-f-continental.m')
477 elif k == 'distances':
478 d['distances'] = num.linspace(10*cake.km, 100*cake.km, 10) \
479 / cake.earthradius * r2d
481 elif k == 'phases':
482 d['phases'] = list(map(cake.PhaseDef, 'Pp'))
484 else:
485 parser.error('missing %s' % k)
487 return Anon(d)
490def my_simplify_model(mod, accuracy):
491 mod_simple = mod.simplify(max_rel_error=accuracy)
492 cake.write_nd_model_fh(mod_simple, sys.stdout)
495def d2u(d):
496 return dict((k.replace('-', '_'), v) for (k, v) in d.items())
499def mini_fmt(v, d=5, try_fmts='fe'):
500 for fmt in try_fmts:
501 for i in range(d, -1, -1):
502 s = ('%%%i.%i%s' % (d, i, fmt)) % v
503 if len(s) <= d and (v == 0.0 or float(s) != 0.0):
504 return s
506 return s
509def scatter_out_fmt(d, m, v):
510 if v < 1e-8:
511 return ' ', ' '
512 else:
513 return '%s%s' % ('\\/'[d == cake.UP], 'SP'[m == cake.P]), \
514 mini_fmt(v*100, 5)
517def scatter_in_fmt(d, m, dwant):
518 if d == dwant:
519 return '%s%s' % ('\\/'[d == cake.UP], 'SP'[m == cake.P])
520 else:
521 return ' '
524def print_scatter(model, p=0.0, interface=None):
525 if interface is not None:
526 discontinuities = [model.discontinuity(interface)]
527 else:
528 discontinuities = model.discontinuities()
530 for discontinuity in discontinuities:
531 print('%s (%g km)' % (discontinuity, discontinuity.z/cake.km))
532 print()
533 cols = []
534 for in_direction in (cake.DOWN, cake.UP):
535 for in_mode in (cake.P, cake.S):
536 p_critical = discontinuity.critical_ps(in_mode)[
537 in_direction == cake.UP]
539 if p_critical is None or p >= p_critical:
540 continue
542 vals = []
543 for out_direction in (cake.UP, cake.DOWN):
544 for out_mode in (cake.S, cake.P):
545 vals.append(
546 (out_direction, out_mode, discontinuity.efficiency(
547 in_direction, out_direction,
548 in_mode, out_mode, p)))
550 if all(x[-1] == 0.0 for x in vals):
551 continue
553 sout = [scatter_out_fmt(d, m, v) for (d, m, v) in vals]
555 sin1 = scatter_in_fmt(in_direction, in_mode, cake.DOWN)
556 sin2 = scatter_in_fmt(in_direction, in_mode, cake.UP)
557 line1 = '%s %5s %5s' % (
558 ' '*len(sin1), sout[0][1], sout[1][1])
559 line2 = '%s %-5s %-5s' % (
560 sin1, sout[0][0], sout[1][0])
561 line4 = '%s %-5s %-5s' % (
562 sin2, sout[2][0], sout[3][0])
563 line5 = '%s %5s %5s' % (
564 ' '*len(sin2), sout[2][1], sout[3][1])
565 line3 = '-' * len(line1)
566 cols.append((line1, line2, line3, line4, line5))
568 for cols in zip(*cols):
569 print(' ' + ' '.join(cols))
571 print()
572 print()
575def print_arrivals(
576 model, distances=[], phases=cake.PhaseDef('P'), zstart=0.0, zstop=0.0,
577 as_degrees=False):
579 headers = 'slow dist time take inci effi spre phase used'.split()
580 space = (7, 5, 6, 4, 4, 4, 4, 17, 17)
582 if as_degrees:
583 units = 's/deg deg s deg deg % %'.split()
584 else:
585 units = 's/km km s deg deg % %'.split()
587 hline = ' '.join(x.ljust(s) for (x, s) in zip(headers, space))
588 uline = ' '.join(('%s' % x).ljust(s) for (x, s) in zip(units, space))
590 print(hline)
591 print(uline)
592 print('-' * len(hline))
594 for ray in model.arrivals(
595 distances=distances, phases=phases, zstart=zstart, zstop=zstop):
597 if as_degrees:
598 sd = ray.x
599 slow = ray.p/cake.r2d
600 else:
601 sd = ray.x*(cake.d2r*cake.earthradius/cake.km)
602 slow = ray.p/(r2d*cake.d2m/cake.km)
604 su = '(%s)' % ray.path.used_phase(p=ray.p, eps=1.0).used_repr()
606 print(' '.join(tuple(mini_fmt(x, s).rjust(s) for (x, s) in zip(
607 (slow, sd, ray.t, ray.takeoff_angle(), ray.incidence_angle(),
608 100*ray.efficiency(), 100*ray.spreading()*ray.surface_sphere()),
609 space)) + tuple(
610 x.ljust(17) for x in (ray.path.phase.definition(), su))))
613def plot_init(size, save, show):
614 fontsize = 9
615 mpl_init()
616 fig = plt.figure(figsize=mpl_papersize(size, 'landscape'))
618 labelpos = mpl_margins(fig, w=7., h=5., units=fontsize)
619 axes = fig.add_subplot(1, 1, 1)
620 labelpos(axes, 2., 1.5)
622 showplt = bool(show or not save)
624 return fig, axes, showplt
627class CakeError(Exception):
628 pass
631def plot_end(save, fig, show=True):
632 if save:
633 try:
634 fig.savefig(save)
635 if show:
636 plt.show()
637 except OSError as e:
638 raise CakeError(str(e))
639 except ValueError as e:
640 raise CakeError(str(e))
643def main(args=None):
645 if args is None:
646 args = sys.argv[1:]
648 subcommand_descriptions = {
649 'print': 'get information on model/phase/material properties',
650 'arrivals': 'print list of phase arrivals',
651 'paths': 'print ray path details',
652 'plot-xt': 'plot traveltime vs distance curves',
653 'plot-xp': 'plot ray parameter vs distance curves',
654 'plot-rays': 'plot ray propagation paths',
655 'plot': 'plot combination of ray and traveltime curves',
656 'plot-model': 'plot velocity model',
657 'list-models': 'list builtin velocity models',
658 'list-phase-map': 'show translation table for classic phase names',
659 'simplify-model': 'create a simplified version of a layered model',
660 'scatter': 'show details about scattering at model interfaces'}
662 usage = '''cake <subcommand> [options]
664Subcommands:
666 print %(print)s
667 arrivals %(arrivals)s
668 paths %(paths)s
669 plot-xt %(plot_xt)s
670 plot-xp %(plot_xp)s
671 plot-rays %(plot_rays)s
672 plot %(plot)s
673 plot-model %(plot_model)s
674 list-models %(list_models)s
675 list-phase-map %(list_phase_map)s
676 simplify-model %(simplify_model)s
677 scatter %(scatter)s
679To get further help and a list of available options for any subcommand run:
681 cake <subcommand> --help
683'''.strip() % d2u(subcommand_descriptions)
685 usage_sub = 'cake %s [options]'
686 if len(args) < 1:
687 sys.exit('Usage: %s' % usage)
689 command = args[0]
691 args[0:0] = ['cake']
692 descr = subcommand_descriptions.get(command, None)
693 subusage = usage_sub % command
695 if command == 'print':
696 c = optparse(
697 (), ('model', 'phases', 'material', 'output_format'),
698 usage=subusage, descr=descr, args=args)
700 if 'model' in c:
701 if c.output_format == 'textual':
702 print(c.model)
703 print()
704 elif c.output_format == 'nd':
705 cake.write_nd_model_fh(c.model, sys.stdout)
707 if 'phases' in c:
708 for phase in c.phases:
709 print(phase)
710 print()
712 if 'material' in c:
713 print(c.material.describe())
714 print()
716 elif command == 'arrivals':
717 c = optparse(
718 ('model', 'phases', 'distances'),
719 ('zstart', 'zstop', 'as_degrees'),
720 usage=subusage, descr=descr, args=args)
722 print_arrivals(
723 c.model,
724 **c.getn('zstart', 'zstop', 'phases', 'distances', 'as_degrees'))
726 elif command == 'paths':
727 c = optparse(
728 ('model', 'phases'),
729 ('zstart', 'zstop', 'as_degrees'),
730 usage=subusage, descr=descr, args=args)
732 mod = c.model
733 for path in mod.gather_paths(**c.getn('phases', 'zstart', 'zstop')):
734 print(path.describe(path.endgaps(c.zstart, c.zstop), c.as_degrees))
736 elif command in ('plot-xt', 'plot-xp', 'plot-rays', 'plot'):
737 if command in ('plot-xt', 'plot'):
738 c = optparse(
739 ('model', 'phases',),
740 ('zstart', 'zstop', 'distances', 'as_degrees', 'vred',
741 'phase_colors', 'save', 'size', 'show'),
742 usage=subusage, descr=descr, args=args)
743 else:
744 c = optparse(
745 ('model', 'phases'),
746 ('zstart', 'zstop', 'distances', 'as_degrees', 'aspect',
747 'shade_model', 'phase_colors', 'save', 'size', 'show'),
748 usage=subusage, descr=descr, args=args)
750 mod = c.model
751 paths = mod.gather_paths(**c.getn('phases', 'zstart', 'zstop'))
753 if c.distances is not None:
754 arrivals = mod.arrivals(
755 **c.getn('phases', 'zstart', 'zstop', 'distances'))
756 else:
757 arrivals = None
759 fig, axes, showplt = plot_init(c.size, c.save, c.show)
761 if command == 'plot-xp':
762 plot.my_xp_plot(
763 paths, c.zstart, c.zstop, c.distances,
764 c.as_degrees, show=showplt, phase_colors=c.phase_colors)
766 elif command == 'plot-xt':
767 plot.my_xt_plot(
768 paths, c.zstart, c.zstop, c.distances, c.as_degrees,
769 vred=c.vred, show=showplt,
770 phase_colors=c.phase_colors)
772 elif command == 'plot-rays':
773 if c.as_degrees:
774 plot.my_rays_plot_gcs(
775 mod, paths, arrivals, c.zstart, c.zstop, c.distances,
776 show=showplt, phase_colors=c.phase_colors)
778 else:
779 plot.my_rays_plot(
780 mod, paths, arrivals, c.zstart, c.zstop, c.distances,
781 show=showplt, aspect=c.aspect, shade_model=c.shade_model,
782 phase_colors=c.phase_colors)
784 elif command == 'plot':
785 plot.my_combi_plot(
786 mod, paths, arrivals, c.zstart, c.zstop, c.distances,
787 c.as_degrees, show=showplt, vred=c.vred,
788 phase_colors=c.phase_colors)
790 try:
791 plot_end(save=c.save, fig=fig, show=c.show)
792 except CakeError as e:
793 exit('cake.py: %s' % str(e))
795 elif command in ('plot-model',):
796 c = optparse(
797 ('model',), ('save', 'size', 'show'),
798 usage=subusage, descr=descr, args=args)
799 mod = c.model
800 fig, axes, showplt = plot_init(c.size, c.save, c.show)
801 plot.my_model_plot(mod, show=showplt)
802 try:
803 plot_end(save=c.save, fig=fig, show=c.show)
804 except CakeError as e:
805 exit('cake.py: %s' % str(e))
807 elif command in ('simplify-model',):
808 c = optparse(('model',), ('accuracy',),
809 usage=subusage, descr=descr, args=args)
810 my_simplify_model(c.model, c.accuracy)
812 elif command in ('list-models',):
813 c = optparse((), (), usage=subusage, descr=descr, args=args)
814 for x in cake.builtin_models():
815 print(x)
817 elif command in ('list-phase-map',):
818 c = optparse((), (), usage=subusage, descr=descr, args=args)
819 defs = cake.PhaseDef.classic_definitions()
820 for k in sorted(defs.keys()):
821 print('%-15s: %s' % (k, ', '.join(defs[k])))
823 elif command in ('scatter',):
824 c = optparse(
825 ('model',),
826 ('slowness', 'interface', 'as_degrees'),
827 usage=subusage, descr=descr, args=args)
829 print_scatter(c.model, p=c.slowness, interface=c.interface)
831 elif command in ('help-options',):
832 optparse(
833 (),
834 ('model', 'accuracy', 'slowness', 'interface', 'phases',
835 'distances', 'zstart', 'zstop', 'distances', 'as_degrees',
836 'material', 'vred', 'save'),
837 usage='cake help-options', descr='list all available options',
838 args=args)
840 elif command in ('--help', '-h', 'help'):
841 sys.exit('Usage: %s' % usage)
843 else:
844 sys.exit('cake: no such subcommand: %s' % command)
847if __name__ == '__main__':
848 main()