# Traveltime calculation and raytracing with `cake`¶

## Calculate synthetic traveltimes¶

Here we will excercise two example how to calculate traveltimes for the phases `P` and `Pg` for different earth velocity models.

Modules covered in this example:

The first example is minimalistic and will give you a simple travel time table.

Download `cake_ray_tracing.py`

```'''
Calculate P-phase arrivals.
'''
from __future__ import print_function
from pyrocko import cake
import numpy as num

km = 1000.

# Load builtin 'prem-no-ocean' model (medium resolution)

# Source depth [m].
source_depth = 300. * km

# Distances as a numpy array [deg].
distances = num.linspace(1500, 3000, 16)*km * cake.m2d

# Define the phase to use.
Phase = cake.PhaseDef('P')

# calculate distances and arrivals and print them:
print('distance [km]      time [s]')
for arrival in model.arrivals(distances, phases=Phase, zstart=source_depth):
print('%13g %13g' % (arrival.x*cake.d2m/km, arrival.t))
```

The second code snippet includes some lines to plot a simple traveltime figure.

Download `cake_ray_tracing.py`

```import numpy as num
from matplotlib import pyplot as plt

from pyrocko.plot import mpl_init, mpl_margins, mpl_papersize
from pyrocko import cake

fontsize = 9.
km = 1000.

phases = cake.PhaseDef.classic('Pg')
source_depth = 20.*km

distances = num.linspace(100.*km, 1000.*km, 100)

data = []
for distance in distances:
rays = mod.arrivals(
phases=phases, distances=[distance*cake.m2d], zstart=source_depth)
for ray in rays[:1]:
data.append((distance, ray.t))

phase_distances, phase_time = num.array(data, dtype=float).T

# Plot the arrival times
mpl_init(fontsize=fontsize)
fig = plt.figure(figsize=mpl_papersize('a5', 'landscape'))
labelpos = mpl_margins(fig, w=7., h=5., units=fontsize)
labelpos(axes, 2., 1.5)

axes.set_xlabel('Distance [km]')
axes.set_ylabel('Time [s]')

axes.plot(phase_distances/km, phase_time, 'o', ms=3., color='black')

fig.savefig('cake_first_arrivals.pdf')
```

## Travel time table interpolation¶

This example demonstrates how to interpolate and query travel time tables.

Classes covered in this example:

Download `cake_raytracing.py`

```import numpy as num
import matplotlib.pyplot as plt
from pyrocko import spit, cake
from pyrocko.gf import meta

# Define a list of phases.
phase_defs = [meta.TPDef(id='stored:depth_p', definition='p'),
meta.TPDef(id='stored:P', definition='P')]

# Load a velocity model. In this example use the default AK135.

# Time and space tolerance thresholds defining the accuracy of the
# :py:class:`pyrocko.spit.SPTree`.
t_tolerance = 0.1                           # in seconds
x_tolerance = num.array((500., 500.))       # in meters

# Boundaries of the grid.
xmin = 0.
xmax = 20000
zmin = 0.
zmax = 11000
x_bounds = num.array(((xmin, xmax), (zmin, zmax)))

# In this example the receiver is located at the surface.

interpolated_tts = {}

for phase_def in phase_defs:
v_horizontal = phase_def.horizontal_velocities

def evaluate(args):
'''Calculate arrival using source and receiver location
defined by *args*. To be evaluated by the SPTree instance.'''
source_depth, x = args

t = []

# Calculate arrivals
rays = mod.arrivals(
phases=phase_def.phases,
distances=[x*cake.m2d],
zstart=source_depth,

for ray in rays:
t.append(ray.t)

for v in v_horizontal:
t.append(x/(v*1000.))
if t:
return min(t)
else:
return None

# Creat a :py:class:`pyrocko.spit.SPTree` interpolator.
sptree = spit.SPTree(
f=evaluate,
ftol=t_tolerance,
xbounds=x_bounds,
xtols=x_tolerance)

# Store the result in a dictionary which is later used to retrieve an
# SPTree (value) for each phase_id (key).
interpolated_tts[phase_def.id] = sptree

# Dump the sptree for later reuse:
sptree.dump(filename='sptree_%s.yaml' % phase_def.id.split(':'))

# Define a :py:class:`pyrocko.gf.meta.Timing` instance.
timing = meta.Timing('first(depth_p|P)')

# If only one interpolated onset is need at a time you can retrieve
# that value as follows:
# First argument has to be a function which takes a requested *phase_id*
# and returns the associated :py:class:`pyrocko.spit.SPTree` instance.
# Second argument is a tuple of distance and source depth.
z_want = 5000.
x_want = 2000.
one_onset = timing.evaluate(lambda x: interpolated_tts[x],
(z_want, x_want))
print('a single arrival: ', one_onset)

# But if you have many locations for which you would like to calculate the
# onset time the following is the preferred way as it is much faster
# on large coordinate arrays.
# x_want is now an array of 1000 distances
x_want = num.linspace(0, xmax, 1000)

# Coords is set up of distances-depth-pairs
coords = num.array((x_want, num.tile(z_want, x_want.shape))).T

# *interpolate_many* then interpolates onset times for each of these
# pairs.
tts = interpolated_tts["stored:depth_p"].interpolate_many(coords)

# Plot distance vs. onset time
plt.plot(x_want, tts, '.')
plt.xlabel('Distance [m]')
plt.ylabel('Travel Time [s]')
plt.show()
```