UPDATE:: Migrate demo notebook to current API

Author: Malcolm White

jupyter/demo.ipynb

@@ -40,16 +40,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def tt_head(xs, ys, zs, xr, yr, zr, y1, v1, v2):\n",
-    "    theta = np.arcsin(v1/v2)\n",
-    "    l1    = (y1 - ys) / np.cos(theta)\n",
-    "    l2    = np.sqrt((xr-xs)**2 + (zr-zs)**2) - np.tan(theta) * (2 * y1 - ys - yr)\n",
-    "    l3    = (y1 - yr) / np.cos(theta)\n",
-    "    return ((l1 + l3) / v1 + l2 / v2)\n",
-    "\n",
-    "def tt_direct(xs, ys, zs, xr, yr, zr, v1):\n",
-    "    return (np.sqrt((xr-xs)**2 + (yr-ys)**2 + (zr-zs)**2) / v1)\n",
-    "\n",
     "def plot(solver, ix=None, iy=None, iz=None, attr='uu', rays=None, cbar_label='Travel-time [s]'):\n",
     "    grid                = pykonal.GridND(ndim=3)\n",
     "    grid.min_coords     = solver.pgrid.min_coords\n",
@@ -140,6 +130,36 @@
     "    fig.tight_layout()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Define function to instantiate EikonalSolver with uniform velocity model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def init_solver():\n",
+    "    # Initialize the solver\n",
+    "    solver = pykonal.EikonalSolver()\n",
+    "\n",
+    "    # Initialize the velocity grid with a uniform velocity model.\n",
+    "    # EikonalSolver.vgrid.min_coords specifies the minimum coordinates of the velocity grid\n",
+    "    solver.vgrid.min_coords     = 0, 0, 0    # xmin, ymin, zmin\n",
+    "    # EikonalSolver.vgrid.node_intervals specifies the spacing between velocity grid nodes\n",
+    "    solver.vgrid.node_intervals = 1, 1, 1    # dx, dy, dz\n",
+    "    # EikonalSolver.vgrid.npts specifies the number of grid nodes along each axis\n",
+    "    solver.vgrid.npts           = 11, 11, 11 # nx, ny, nz\n",
+    "    # EikonalSolver.vv holds the velocity at each grid node and should be a numpy.ndarray\n",
+    "    # with shape == EikonalSolver.vgrid.npts.\n",
+    "    solver.vv                   = np.ones(solver.vgrid.npts)\n",
+    "    return (solver)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -154,29 +174,13 @@
    "outputs": [],
    "source": [
     "# Initialize the solver\n",
-    "solver = pykonal.EikonalSolver()\n",
-    "\n",
-    "# Initialize the velocity grid with a uniform velocity model.\n",
-    "# EikonalSolver.vgrid.min_coords specifies the minimum coordinates of the velocity grid\n",
-    "solver.vgrid.min_coords     = 0, 0, 0    # xmin, ymin, zmin\n",
-    "# EikonalSolver.vgrid.node_intervals specifies the spacing between velocity grid nodes\n",
-    "solver.vgrid.node_intervals = 1, 1, 1    # dx, dy, dz\n",
-    "# EikonalSolver.vgrid.npts specifies the number of grid nodes along each axis\n",
-    "solver.vgrid.npts           = 11, 11, 11 # nx, ny, nz\n",
-    "# EikonalSolver.vv holds the velocity at each grid node and should be a numpy.ndarray\n",
-    "# with shape == EikonalSolver.vgrid.npts.\n",
-    "solver.vv                   = np.ones(solver.vgrid.npts)\n",
-    "\n",
-    "# Initialize the propagation grid to coincide with velocity grid exactly.\n",
-    "# Attributes of the propagation grid are analogous to those of the velocity grid.\n",
-    "solver.pgrid.min_coords     = solver.vgrid.min_coords\n",
-    "solver.pgrid.node_intervals = solver.vgrid.node_intervals\n",
-    "solver.pgrid.npts           = solver.vgrid.npts\n",
+    "solver = init_solver()\n",
     "\n",
     "# Add a source in the center of the computational domain.\n",
-    "# Source locations need to be specified in physical coordinates.\n",
-    "src = (5, 5, 5)\n",
-    "solver.add_source(src)\n",
+    "src_idx = (5, 5, 5)            # The source location as an array index\n",
+    "solver.uu[src_idx]     = 0     # Set the traveltime at the source location to zero\n",
+    "solver.is_far[src_idx] = False # Set the is_far flag to False for the source node\n",
+    "solver.close.push(*src_idx)    # Push the source index onto the close heap\n",
     "\n",
     "# Solve the Eikonal equation.\n",
     "solver.solve()\n",
@@ -199,12 +203,21 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# Initialize the solver\n",
+    "solver = init_solver()\n",
     "# Decrease the node interval by a factor of 2\n",
     "solver.pgrid.node_intervals = solver.vgrid.node_intervals / 2\n",
     "# And increase the number of points by a factor of 2, making sure to not go beyond the\n",
     "# boundaries of the velocity grid\n",
     "solver.pgrid.npts           = solver.vgrid.npts * 2 - 1\n",
     "\n",
+    "\n",
+    "# Add a source in the center of the computational domain.\n",
+    "src_idx = (10, 10, 10)            # The source location as an array index\n",
+    "solver.uu[src_idx]     = 0     # Set the traveltime at the source location to zero\n",
+    "solver.is_far[src_idx] = False # Set the is_far flag to False for the source node\n",
+    "solver.close.push(*src_idx)    # Push the source index onto the close heap\n",
+    "\n",
     "# Solve the Eikonal equation again.\n",
     "solver.solve()\n",
     "\n",
@@ -226,8 +239,16 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "solver = init_solver()\n",
     "solver.pgrid.node_intervals = solver.vgrid.node_intervals / 10\n",
     "solver.pgrid.npts           = solver.vgrid.npts * 10 - 9\n",
+    "\n",
+    "# Add a source in the center of the computational domain.\n",
+    "src_idx = (50, 50, 50)            # The source location as an array index\n",
+    "solver.uu[src_idx]     = 0     # Set the traveltime at the source location to zero\n",
+    "solver.is_far[src_idx] = False # Set the is_far flag to False for the source node\n",
+    "solver.close.push(*src_idx)    # Push the source index onto the close heap\n",
+    "\n",
     "solver.solve()\n",
     "plot(solver)"
    ]
@@ -245,8 +266,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "solver.clear_sources()\n",
-    "solver.add_source((0, 5, 10))\n",
+    "solver = init_solver()\n",
+    "solver.pgrid.node_intervals = solver.vgrid.node_intervals / 10\n",
+    "solver.pgrid.npts           = solver.vgrid.npts * 10 - 9\n",
+    "\n",
+    "src_idx = (25, 50, 75)            # The source location as an array index\n",
+    "solver.uu[src_idx]     = 0     # Set the traveltime at the source location to zero\n",
+    "solver.is_far[src_idx] = False # Set the is_far flag to False for the source node\n",
+    "solver.close.push(*src_idx)    # Push the source index onto the close heap\n",
+    "\n",
     "solver.solve()\n",
     "plot(solver)"
    ]
@@ -264,7 +292,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "solver.add_source((10, 5, 0))\n",
+    "solver = init_solver()\n",
+    "solver.pgrid.node_intervals = solver.vgrid.node_intervals / 10\n",
+    "solver.pgrid.npts           = solver.vgrid.npts * 10 - 9\n",
+    "\n",
+    "for src_idx in ((25, 50, 75), (75, 50, 25)):\n",
+    "    solver.uu[src_idx]     = 0     # Set the traveltime at the source location to zero\n",
+    "    solver.is_far[src_idx] = False # Set the is_far flag to False for the source node\n",
+    "    solver.close.push(*src_idx)    # Push the source index onto the close heap\n",
+    "\n",
     "solver.solve()\n",
     "plot(solver)"
    ]
@@ -282,7 +318,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "solver.add_source((5, 5, 5), t0=2.5)\n",
+    "solver = init_solver()\n",
+    "solver.pgrid.node_intervals = solver.vgrid.node_intervals / 10\n",
+    "solver.pgrid.npts           = solver.vgrid.npts * 10 - 9\n",
+    "\n",
+    "for src_idx, t0 in (((25, 50, 75), 0), ((75, 50, 25), 2.5)):\n",
+    "    solver.uu[src_idx]     = t0     # Set the traveltime at the source location to t0\n",
+    "    solver.is_far[src_idx] = False # Set the is_far flag to False for the source node\n",
+    "    solver.close.push(*src_idx)    # Push the source index onto the close heap\n",
+    "\n",
+    "solver.solve()\n",
+    "plot(solver)\n",
     "solver.solve()\n",
     "plot(solver)"
    ]
@@ -300,32 +346,24 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "solver = init_solver()\n",
+    "\n",
     "vy = np.linspace(1, 5, solver.vgrid.npts[1])\n",
     "for iy in range(len(vy)):\n",
     "    solver.vv[:,iy] = vy[iy]\n",
     "    \n",
-    "# This will plot the velocity model\n",
-    "plot(solver, attr='vv_p', cbar_label='Velocity [km/s]')\n",
-    "\n",
-    "solver.clear_sources()\n",
-    "solver.add_source((5, 5, 5))\n",
-    "solver.solve()\n",
-    "plot(solver)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "solver.vv[:, int(solver.vgrid.npts[1]//2):] = 5\n",
+    "solver.pgrid.node_intervals = solver.vgrid.node_intervals / 10\n",
+    "solver.pgrid.npts           = solver.vgrid.npts * 10 - 9\n",
     "\n",
-    "plot(solver, attr='vv_p', cbar_label='Velocity [km/s]')\n",
+    "src_idx = (25, 50, 75)\n",
+    "solver.uu[src_idx]     = 0     # Set the traveltime at the source location to zero\n",
+    "solver.is_far[src_idx] = False # Set the is_far flag to False for the source node\n",
+    "solver.close.push(*src_idx)    # Push the source index onto the close heap\n",
     "\n",
-    "solver.clear_sources()\n",
-    "solver.add_source((1, 3, 1))\n",
     "solver.solve()\n",
+    "\n",
+    "# This will plot the velocity model\n",
+    "plot(solver, attr='vvp', cbar_label='Velocity [km/s]')\n",
     "plot(solver)"
    ]
   },
@@ -343,7 +381,7 @@
    "outputs": [],
    "source": [
     "ray = solver.trace_ray((9.5, 0, 9.5))\n",
-    "plot(solver, rays=[ray])"
+    "plot(solver,rays=[ray])"
    ]
   },
   {
@@ -391,5 +429,5 @@
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }