add eikonal examples and improve docs

doc/source/library/examples/cake_raytracing.rst

@@ -1,15 +1,15 @@
-Traveltime calculation and raytracing with ``cake``
-===================================================
+Traveltime calculation and raytracing
+=====================================
 
-Calculate synthetic traveltimes
--------------------------------
+Calculate traveltimes in layered media
+--------------------------------------
 
 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:
  * :py:class:`pyrocko.cake`
 
- The first example is minimalistic and will give you a simple travel time table. 
+ The first example is minimalistic and will give you a simple traveltime table.
 
 Download :download:`cake_ray_tracing.py </../../examples/cake_arrivals.py>`
 
@@ -25,13 +25,48 @@ Download :download:`cake_ray_tracing.py </../../examples/cake_first_arrivals.py>
     :language: python
 
 
-Travel time table interpolation
+
+Calculate traveltimes in heterogeneous media
+--------------------------------------------
+
+These examples demonstrate how to use the :py:mod:`pyrocko.modelling.eikonal`
+module to calculate first arrivals in heterogenous media.
+
+
+.. literalinclude :: /../../examples/eikonal_example1.py
+    :caption: :download:`eikonal_example1.py </../../examples/eikonal_example1.py>`
+    :language: python
+
+.. figure :: /static/eikonal_example1.png
+    :align: center
+    :width: 90%
+    :alt: output of eikonal_example1.py
+
+    First arrivals (contours) from a seismic source (star) at 15 km depth in a
+    5-layer crustal model where velocities increase with depth.
+
+.. literalinclude :: /../../examples/eikonal_example2.py
+    :caption: :download:`eikonal_example2.py </../../examples/eikonal_example2.py>`
+    :language: python
+
+.. figure :: /static/eikonal_example2.png
+    :align: center
+    :width: 90%
+    :alt: output of eikonal_example2.py
+
+    First arrivals (contours) from a distant seismic source in a 2-layer
+    crustal model with intrusions. The planar wave front entering from below is
+    simulated by a source at 10 km depth moving quickly from left to right with
+    a given constant speed.
+
+
+Traveltime table interpolation
 -------------------------------
 
-This example demonstrates how to interpolate and query travel time tables.
+This example demonstrates how to interpolate and query traveltime tables.
 
 Classes covered in this example:
- * :py:class:`pyrocko.spit.SPTree` (interpolation of travel time tables)
+ * :py:class:`pyrocko.spit.SPTree` (interpolation of traveltime tables)
  * :py:class:`pyrocko.gf.meta.TPDef` (phase definitions)
  * :py:class:`pyrocko.gf.meta.Timing` (onset definition to query the travel
    time tables)
@@ -40,4 +75,3 @@ Download :download:`cake_raytracing.py </../../examples/cake_raytracing.py>`
 
 .. literalinclude :: /../../examples/cake_raytracing.py
     :language: python
-

examples/eikonal_example1.py

@@ -0,0 +1,45 @@
+import numpy as num
+from matplotlib import pyplot as plt
+from pyrocko.modelling import eikonal
+
+km = 1000.
+nx, ny = 1500, 500                # grid size
+delta = 90*km / float(nx)         # grid spacing
+source_x, source_y = 0.0, 15*km   # source position
+
+# Indexing arrays
+x = num.arange(nx) * delta - 2*km
+y = num.arange(ny) * delta
+x2 = x[num.newaxis, :]
+y2 = y[:, num.newaxis]
+
+# Define layers with different speeds, roughly representing a crustal model.
+speeds = num.ones((ny, nx))
+nlayer = ny // 5
+speeds[0*nlayer:1*nlayer, :] = 2500.
+speeds[1*nlayer:2*nlayer, :] = 3500.
+speeds[2*nlayer:3*nlayer, :] = 5000.
+speeds[3*nlayer:4*nlayer, :] = 6000.
+speeds[4*nlayer:, :] = 8000.
+
+# Seeding points have non-negative times. Here we simply set one grid node to
+# zero. The solution to the eikonal equation is computed at all nodes where
+# times < 0.
+times = num.zeros((ny, nx)) - 1.0
+times[int(round(source_y/delta)), int(round((source_x-x[0])//delta))] = 0.0
+
+# Solve eikonal equation.
+eikonal.eikonal_solver_fmm_cartesian(speeds, times, delta)
+
+# Plot
+fig = plt.figure(figsize=(9.0, 4.0))
+
+axes = fig.add_subplot(1, 1, 1, aspect=1.0)
+axes.contourf(x/km, y/km, times)
+axes.invert_yaxis()
+axes.contourf(x/km, y/km, speeds, alpha=0.1, cmap='gray')
+axes.plot(source_x/km, source_y/km, '*', ms=20, color='white')
+axes.set_xlabel('Distance [km]')
+axes.set_ylabel('Depth [km]')
+# fig.savefig('eikonal_example1.png')
+plt.show()

examples/eikonal_example2.py

@@ -0,0 +1,42 @@
+import numpy as num
+from matplotlib import pyplot as plt
+from pyrocko.modelling import eikonal
+
+km = 1000.
+nx, ny = 1000, 500          # grid size
+delta = 20*km / float(nx)   # drid spacing
+
+# Indexing arrays
+x = num.arange(nx) * delta - 10.0
+y = num.arange(ny) * delta
+x2 = x[num.newaxis, :]
+y2 = y[:, num.newaxis]
+
+# Define layers and circles with different speeds, roughly representing a case
+# with two layers and intrusions.
+speeds = num.ones((ny, nx))
+r1 = num.sqrt((x2-0*km)**2 + (y2-2*km)**2)
+r2 = num.sqrt((x2-12*km)**2 + (y2-0*km)**2)
+nlayer = ny // 5
+speeds[r1 < 4*km] = 2.0
+speeds[r2 < 4*km] = 0.7
+speeds[:3*nlayer, :] *= 0.5
+
+# Seeding points have non-negative times. Here we
+times = num.zeros((ny, nx)) - 1.0
+times[-1, :] = (x-num.min(x)) * 0.1
+
+# Solve eikonal equation.
+eikonal.eikonal_solver_fmm_cartesian(speeds, times, delta)
+
+# Plot
+fig = plt.figure(figsize=(9.0, 4.0))
+
+axes = fig.add_subplot(1, 1, 1, aspect=1.0)
+axes.contourf(x/km, y/km, times)
+axes.invert_yaxis()
+axes.contourf(x/km, y/km, speeds, alpha=0.1, cmap='gray')
+axes.set_xlabel('Distance [km]')
+axes.set_ylabel('Depth [km]')
+# fig.savefig('eikonal_example2.png')
+plt.show()

src/modelling/eikonal.py

@@ -0,0 +1,37 @@
+# https://pyrocko.org - GPLv3
+#
+# The Pyrocko Developers, 21st Century
+# ---|P------/S----------~Lg----------
+
+from .. import eikonal_ext
+
+
+def eikonal_solver_fmm_cartesian(speeds, times, delta):
+    '''
+    Solve eikonal equation in 2D or 3D using the fast marching method.
+
+    This function implements the fast marching method (FMM) by [sethian1996]_.
+
+    :param speeds:
+        Velocities at the grid nodes.
+    :type speeds:
+        2D or 3D :py:class:`numpy.ndarray`
+
+    :param times:
+        Arrival times (input and output). The solution is obtained at nodes
+        where times is set to a negative value. Values of zero, or positive
+        values are used as seeding points.
+    :type times:
+        2D or 3D :py:class:`numpy.ndarray`, same shape as `speeds`
+
+    :param delta:
+        Grid spacing.
+    :type delta:
+        float
+
+    .. [sethian1996] Sethian, James A. "A fast marching level set method for
+        monotonically advancing fronts." Proceedings of the National Academy of
+        Sciences 93.4 (1996): 1591-1595. https://doi.org/10.1073/pnas.93.4.1591
+    '''
+
+    return eikonal_ext.eikonal_solver_fmm_cartesian(speeds, times, delta)