diff --git a/book/3_diffusion_fields/em_diffusion.md b/book/3_diffusion_fields/em_diffusion.md
index 00d54f7..26e3532 100644
--- a/book/3_diffusion_fields/em_diffusion.md
+++ b/book/3_diffusion_fields/em_diffusion.md
@@ -208,7 +208,7 @@ $$
\begin{aligned}
-\partial_y H_z + \sigma E_x &= 0, \\
\partial_x H_z + \sigma E_y &= 0, \\
-\partial_x H_y - \partial_y H_x + \sigma E_z &= -J_z^e, \\
+-\partial_x H_y + \partial_y H_x + \sigma E_z &= -J_z^e, \\
\partial_y E_z + \mu\,\partial_t H_x &= 0, \\
-\partial_x E_z + \mu\,\partial_t H_y &= 0, \\
\partial_x E_y - \partial_y E_x + \mu\,\partial_t H_z &= 0 .
@@ -218,7 +218,7 @@ $$
The fields that connect to the source are grouped as
$$
-\partial_x H_y - \partial_y H_x + \sigma E_z = -J_z^e,
+-\partial_x H_y + \partial_y H_x + \sigma E_z = -J_z^e,
$$ (eq:TEEy)
$$
@@ -248,7 +248,7 @@ We are not going to derive solutions but give the electric and magnetic fields f
$$
\begin{aligned}
-E_z(x,y,t) &= \frac{\mu I}{4\pi t}\exp\left(-\frac{\sigma\mu\varrho^2}{4t}\right), \\
+E_z(x,y,t) &= -\frac{\mu I}{4\pi t}\exp\left(-\frac{\sigma\mu\varrho^2}{4t}\right), \\
H_x(x,y,t) &= -\frac{\mu y I}{2\pi\varrho^2}\exp\left(-\frac{\sigma\mu\varrho^2}{4t}\right), \\
H_y(x,y,t) &= \frac{\mu x I}{2\pi\varrho^2}\exp\left(-\frac{\sigma\mu\varrho^2}{4t}\right).
\end{aligned}
diff --git a/book/3_diffusion_fields/labs/figures/ELoopwsxz.mp4 b/book/3_diffusion_fields/labs/figures/ELoopwsxz.mp4
new file mode 100644
index 0000000..75c8522
Binary files /dev/null and b/book/3_diffusion_fields/labs/figures/ELoopwsxz.mp4 differ
diff --git a/book/3_diffusion_fields/labs/figures/make_week06_figures.py b/book/3_diffusion_fields/labs/figures/make_week06_figures.py
new file mode 100644
index 0000000..45a98b6
--- /dev/null
+++ b/book/3_diffusion_fields/labs/figures/make_week06_figures.py
@@ -0,0 +1,240 @@
+"""Regenerate the four model-setup schematics of the Week 6 lab.
+
+Run: python make_week06_figures.py
+
+Same palette and conventions as make_figures.py: transparent SVG, grey ink so the
+theme can invert it in dark mode, TU Delft cyan for the field under study and
+orange for the source. Depth runs downward throughout, as it does on every figure
+in the chapter. Long captions go in figure coordinates, so they never fight the
+equal aspect ratio the circular panels need.
+"""
+import matplotlib
+matplotlib.use("Agg")
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.patches import Rectangle, FancyArrowPatch, Circle, Polygon, Arc
+
+INK, GND, SRC, FLD = "#6b7280", "#8a9199", "#EC6842", "#00A6D6"
+
+
+def arrow(ax, p0, p1, color=SRC, lw=2.0, scale=13, z=5):
+ ax.add_patch(FancyArrowPatch(p0, p1, arrowstyle="-|>", mutation_scale=scale,
+ color=color, lw=lw, zorder=z, shrinkA=0, shrinkB=0))
+
+
+def bare(ax, xlim, ylim, equal=False):
+ ax.set_xlim(*xlim); ax.set_ylim(*ylim)
+ ax.set_xticks([]); ax.set_yticks([])
+ for s in ax.spines.values():
+ s.set_visible(False)
+ if equal:
+ ax.set_aspect("equal")
+
+
+def caption(fig, text, y=0.035, size=8.5):
+ fig.text(0.5, y, text, ha="center", va="bottom", color=INK, fontsize=size)
+
+
+# ---------------------------------------------------------------------------
+# Part 1: an infinite current sheet, seen edge-on
+# ---------------------------------------------------------------------------
+def sheet():
+ fig, ax = plt.subplots(figsize=(6.8, 3.1))
+ ax.add_patch(Rectangle((-1.0, -1.0), 2.0, 1.0, fc=GND, ec="none", alpha=0.16))
+ ax.text(0.97, -0.96, "conductor $\\sigma,\\ \\mu$", ha="right", va="bottom",
+ color=INK, fontsize=9)
+
+ # the sheet itself
+ ax.plot([-1.0, 1.0], [0, 0], color=SRC, lw=3, solid_capstyle="butt", zorder=4)
+ for x in np.linspace(-0.88, 0.72, 7):
+ arrow(ax, (x, 0), (x + 0.16, 0), scale=9, lw=1.3)
+ ax.text(0.0, 0.07, "current sheet at $z=0$: current $I$ per unit width along $x$,"
+ " infinite in $x$ and $y$",
+ ha="center", va="bottom", color=SRC, fontsize=9)
+
+ # depth axis
+ arrow(ax, (-1.18, 0.0), (-1.18, -0.95), color=INK, lw=1.2, scale=11)
+ ax.text(-1.23, -0.48, "depth $z$", ha="right", va="center", color=INK,
+ fontsize=9.5, rotation=90)
+
+ # spreading front, then the three receivers
+ for z in np.linspace(-0.10, -0.94, 8): # the spreading front, left margin only
+ ax.plot([-0.97, -0.45], [z, z], color=FLD, lw=0.9, ls=":", alpha=0.55, zorder=2)
+ for z, lab in ((-0.22, "100 m"), (-0.45, "300 m"), (-0.82, "1000 m")):
+ ax.add_patch(Polygon([(-0.06, z), (0.06, z), (0.0, z + 0.055)], closed=True,
+ fc=FLD, ec=FLD, zorder=6))
+ ax.text(0.11, z, f"receiver at {lab}", ha="left", va="center", color=FLD,
+ fontsize=9)
+ arrow(ax, (-0.36, z), (-0.11, z), color=FLD, lw=1.5, scale=10)
+ ax.text(-0.40, -0.22, "$E_x$", ha="right", va="center", color=FLD, fontsize=10)
+
+ bare(ax, (-1.45, 1.35), (-1.12, 0.26))
+ caption(fig, "The field diffuses downward. The solution is even in $z$, so reading $z$"
+ " as depth idealises away the air above.")
+ fig.tight_layout(rect=(0, 0.08, 1, 1))
+ fig.savefig("setup_sheet.svg", transparent=True, metadata={"Date": None})
+ plt.close(fig)
+
+
+# ---------------------------------------------------------------------------
+# Part 2: an infinite wire, and the plane the fields live in
+# ---------------------------------------------------------------------------
+def wire():
+ fig, (a0, a1) = plt.subplots(1, 2, figsize=(7.4, 3.2),
+ gridspec_kw=dict(width_ratios=[0.85, 1.15]))
+
+ # left: the wire, with the plane of interest cutting it
+ a0.plot([0, 0], [-1.1, 1.1], color=SRC, lw=3, zorder=5)
+ for y in (-0.8, -0.05, 0.6):
+ arrow(a0, (0, y), (0, y + 0.26), scale=11, lw=1.6)
+ a0.text(0.1, 1.05, "wire along $z$,\ncurrent $I$", ha="left", va="top",
+ color=SRC, fontsize=9.5)
+ a0.add_patch(Polygon([(-0.95, -0.38), (0.5, -0.64), (0.95, -0.02), (-0.5, 0.24)],
+ closed=True, fc=FLD, ec=FLD, alpha=0.20, lw=1.2, zorder=2))
+ a0.text(-1.08, -0.62, "the $(x,y)$-plane", ha="left", va="top", color=FLD,
+ fontsize=9.5)
+ a0.set_title("side view", color=INK, fontsize=9.5, pad=6)
+ bare(a0, (-1.15, 1.35), (-1.35, 1.25))
+
+ # right: that plane, seen along the wire
+ for r in (0.34, 0.60, 0.86):
+ a1.add_patch(Circle((0, 0), r, fc="none", ec=FLD, lw=1.3, ls="--", alpha=0.85))
+ a1.add_patch(Arc((0, 0), 1.20, 1.20, theta1=22, theta2=68, color=FLD, lw=1.9))
+ th = np.deg2rad(68)
+ arrow(a1, (0.60 * np.cos(th), 0.60 * np.sin(th)),
+ (0.60 * np.cos(th + 0.22), 0.60 * np.sin(th + 0.22)), color=FLD, lw=1.9,
+ scale=13)
+ a1.text(0.34, 0.80, "$\\boldsymbol{H}$ circles the wire", ha="left", va="bottom",
+ color=FLD, fontsize=9.5)
+
+ a1.add_patch(Circle((0, 0), 0.085, fc="none", ec=SRC, lw=2, zorder=6))
+ a1.add_patch(Circle((0, 0), 0.026, fc=SRC, ec=SRC, zorder=6))
+ a1.text(0.0, -0.14, "$I$ out, $E_z$ in", ha="center", va="top", color=SRC, fontsize=9)
+
+ ta = np.deg2rad(-32)
+ arrow(a1, (0, 0), (0.86 * np.cos(ta), 0.86 * np.sin(ta)), color=INK, lw=1.3,
+ scale=11)
+ a1.text(0.50, -0.40, "$\\varrho$", ha="left", va="top", color=INK, fontsize=11)
+ a1.set_title("seen along the wire", color=INK, fontsize=9.5, pad=6)
+ bare(a1, (-1.1, 1.25), (-1.1, 1.1), equal=True)
+
+ caption(fig, "Nothing depends on $z$, so the fields live in the $(x,y)$-plane:"
+ " $I$ points out of it and $E_z$ into it, $\\boldsymbol{H}$ lies in it.\n"
+ "Receivers sit at three distances from the wire.", y=0.02)
+ fig.tight_layout(rect=(0, 0.15, 1, 1))
+ fig.savefig("setup_wire.svg", transparent=True, metadata={"Date": None})
+ plt.close(fig)
+
+
+# ---------------------------------------------------------------------------
+# Part 3: the same kernel, spreading over a plane, a cylinder and a sphere
+# ---------------------------------------------------------------------------
+def ladder():
+ fig, axes = plt.subplots(1, 3, figsize=(7.6, 3.0))
+ titles = ("sheet, 1-D", "wire, 2-D", "point, 3-D")
+ pees = ("$p = 1/2$", "$p = 1$", "$p = 3/2$")
+ tps = ("$t_p = ?$", "$t_p = ?$", "$t_p = ?$")
+ spreads = ("spreads one way", "spreads over a circle", "spreads over a sphere")
+
+ for k, ax in enumerate(axes):
+ if k == 0:
+ ax.plot([-0.85, 0.85], [0, 0], color=SRC, lw=3, zorder=5)
+ for z in (0.26, 0.5, 0.74):
+ for s in (1, -1):
+ ax.plot([-0.85, 0.85], [s * z, s * z], color=FLD, lw=1.2, ls="--",
+ alpha=0.85)
+ arrow(ax, (0.0, 0.0), (0.0, -0.72), color=FLD, lw=1.6)
+ ax.text(0.09, -0.38, "$r$", color=FLD, fontsize=11, va="center")
+ else:
+ fills = (0.0, 0.11)[k == 2]
+ for r, al in ((0.30, 0.9), (0.54, 0.65), (0.78, 0.4)):
+ if fills:
+ ax.add_patch(Circle((0, 0), r, fc=FLD, ec="none", alpha=fills))
+ ax.add_patch(Circle((0, 0), r, fc="none", ec=FLD, lw=1.2, ls="--",
+ alpha=al))
+ ax.add_patch(Circle((0, 0), 0.05, fc=SRC, ec=SRC, zorder=6))
+ aa = np.deg2rad(38)
+ arrow(ax, (0, 0), (0.78 * np.cos(aa), 0.78 * np.sin(aa)), color=FLD, lw=1.6)
+ ax.text(0.20, 0.44, "$r$", color=FLD, fontsize=11,
+ ha="right", va="bottom")
+
+ ax.set_title(titles[k], color=INK, fontsize=10.5, pad=8)
+ ax.text(0.0, -1.06, spreads[k], ha="center", va="top", color=INK, fontsize=8.5)
+ ax.text(0.0, -1.34, pees[k], ha="center", va="top", color=SRC, fontsize=11.5)
+ ax.text(0.0, -1.66, tps[k], ha="center", va="top", color=INK, fontsize=10.5)
+ bare(ax, (-1.0, 1.0), (-2.0, 1.0), equal=True)
+
+ caption(fig, "One kernel, $\\exp(-\\sigma\\mu r^2/4t)$. Only the prefactor $t^{-p}$"
+ " changes. Your $t_p(r, p)$ fills in the three question marks.", y=0.025)
+ fig.tight_layout(rect=(0, 0.07, 1, 1))
+ fig.savefig("setup_ladder.svg", transparent=True, metadata={"Date": None})
+ plt.close(fig)
+
+
+# ---------------------------------------------------------------------------
+# Part 4: the loop, and the two planes the lab cuts through it
+# ---------------------------------------------------------------------------
+def loop():
+ # axonometric projection: x right, y toward the viewer (down-left), z down
+ def P(x, y, z):
+ return (x - 0.52 * y, -0.30 * y - 0.92 * z)
+
+ L, zc = 1.0, 0.72 # half-side, and the depth of the (x,y) cut
+ fig, ax = plt.subplots(figsize=(6.8, 3.6))
+
+ # the vertical cut at y = 0, drawn first so the rest sits on top
+ ax.add_patch(Polygon([P(-1.75, 0, -0.05), P(1.75, 0, -0.05), P(1.75, 0, 1.05),
+ P(-1.75, 0, 1.05)], closed=True, fc=FLD, ec=FLD, lw=1.0,
+ alpha=0.13, zorder=1))
+ # the ground surface
+ ax.add_patch(Polygon([P(-1.9, -1.8, 0), P(1.9, -1.8, 0), P(1.9, 1.8, 0),
+ P(-1.9, 1.8, 0)], closed=True, fc=GND, ec=INK, lw=0.8,
+ alpha=0.13, zorder=2))
+ ax.text(*P(1.92, -1.8, 0), s=" surface", ha="left", va="center", color=INK,
+ fontsize=9)
+ # the horizontal cut at z = 10 m
+ ax.add_patch(Polygon([P(-1.7, -1.6, zc), P(1.7, -1.6, zc), P(1.7, 1.6, zc),
+ P(-1.7, 1.6, zc)], closed=True, fc=FLD, ec=FLD, lw=1.1,
+ alpha=0.20, zorder=3))
+ ax.text(*P(1.72, -1.6, zc), s=" the $(x,y)$-plane at $z=10$ m", ha="left",
+ va="center", color=FLD, fontsize=9)
+ ax.text(*P(-1.78, 0, 0.62), s="the $(x,z)$-plane\nat $y=0$ ", ha="right",
+ va="center", color=FLD, fontsize=9)
+
+ # the loop, 1 -> 2 -> 3 -> 4: clockwise seen from above, with z downward
+ verts = [(-L, -L), (L, -L), (L, L), (-L, L)]
+ for k in range(4):
+ x0, y0 = verts[k]
+ x1, y1 = verts[(k + 1) % 4]
+ ax.plot(*zip(P(x0, y0, 0), P(x1, y1, 0)), color=SRC, lw=2.6, zorder=6,
+ solid_capstyle="round")
+ arrow(ax, P(x0 + 0.40 * (x1 - x0), y0 + 0.40 * (y1 - y0), 0),
+ P(x0 + 0.60 * (x1 - x0), y0 + 0.60 * (y1 - y0), 0), scale=13, z=7)
+ for (vx, vy), k in zip(verts, "1234"):
+ px, py = P(vx, vy, 0)
+ ax.text(px + (0.10 if vx > 0 else -0.10), py + 0.06, k, ha="center",
+ va="bottom", color=SRC, fontsize=10.5, zorder=8)
+ ax.text(-2.40, 0.88, "square loop on the surface\n$L_x = L_y = 100$ m\n"
+ "$I$ switched on at $t=0$", ha="left", va="top", color=SRC,
+ fontsize=9.5, zorder=8, linespacing=1.45)
+
+ # axes, from the centre of the loop
+ for (dx, dy, dz), lab, off in (((0.70, 0, 0), "$x$", (0.10, 0.02)),
+ ((0, 0.70, 0), "$y$", (-0.10, -0.05)),
+ ((0, 0, 0.95), "$z$", (0.10, 0.0))):
+ arrow(ax, P(0, 0, 0), P(dx, dy, dz), color=INK, lw=1.3, scale=12, z=9)
+ px, py = P(dx, dy, dz)
+ ax.text(px + off[0], py + off[1], lab, color=INK, fontsize=11, zorder=9,
+ ha="center", va="center")
+
+ bare(ax, (-2.45, 3.15), (-1.45, 0.95))
+ caption(fig, "The field is cut in the vertical plane $y = 0$ and in a horizontal"
+ " plane at depth $z$, with $D = \\sqrt{4t/\\sigma\\mu}$.")
+ fig.tight_layout(rect=(0, 0.08, 1, 1))
+ fig.savefig("setup_loop.svg", transparent=True, metadata={"Date": None})
+ plt.close(fig)
+
+
+if __name__ == "__main__":
+ sheet(); wire(); ladder(); loop()
+ print("wrote setup_sheet.svg, setup_wire.svg, setup_ladder.svg, setup_loop.svg")
diff --git a/book/3_diffusion_fields/labs/figures/setup_ladder.svg b/book/3_diffusion_fields/labs/figures/setup_ladder.svg
new file mode 100644
index 0000000..bd735ad
--- /dev/null
+++ b/book/3_diffusion_fields/labs/figures/setup_ladder.svg
@@ -0,0 +1,1476 @@
+
+
+
diff --git a/book/3_diffusion_fields/labs/figures/setup_loop.svg b/book/3_diffusion_fields/labs/figures/setup_loop.svg
new file mode 100644
index 0000000..fda768e
--- /dev/null
+++ b/book/3_diffusion_fields/labs/figures/setup_loop.svg
@@ -0,0 +1,1263 @@
+
+
+
diff --git a/book/3_diffusion_fields/labs/figures/setup_sheet.svg b/book/3_diffusion_fields/labs/figures/setup_sheet.svg
new file mode 100644
index 0000000..f776855
--- /dev/null
+++ b/book/3_diffusion_fields/labs/figures/setup_sheet.svg
@@ -0,0 +1,1242 @@
+
+
+
diff --git a/book/3_diffusion_fields/labs/figures/setup_wire.svg b/book/3_diffusion_fields/labs/figures/setup_wire.svg
new file mode 100644
index 0000000..c0f7f45
--- /dev/null
+++ b/book/3_diffusion_fields/labs/figures/setup_wire.svg
@@ -0,0 +1,1146 @@
+
+
+
diff --git a/book/3_diffusion_fields/labs/fwtools.py b/book/3_diffusion_fields/labs/fwtools.py
index e354ecf..4f44274 100644
--- a/book/3_diffusion_fields/labs/fwtools.py
+++ b/book/3_diffusion_fields/labs/fwtools.py
@@ -6,8 +6,10 @@
rates, Gaussian kernels -- rather than on plotting boilerplate.
The self-check functions have the same names and behaviour as in the Lab 1
-and Lab 2 module. Everything here is one-dimensional in space: fields are
-arrays of shape (nt, nx), one row per time.
+and Lab 2 module. The one-dimensional helpers take fields of shape (nt, nx),
+one row per time. ``show_map`` and ``show_isosurface`` take two- and
+three-dimensional fields, indexed as ``np.meshgrid(..., indexing='ij')``
+builds them, so ``F[ix, iy]`` sits at ``(x[ix], y[iy])``.
Requires: numpy, matplotlib, plotly.
"""
@@ -21,6 +23,7 @@
__all__ = [
"show_profiles", "show_spacetime", "show_spacetime_slider", "show_records",
+ "show_map", "show_isosurface",
"check", "check_close", "check_abs", "check_scalar",
]
@@ -35,7 +38,8 @@ def _fmt(v, fmt):
def show_profiles(x, frames, values, *, value_name="t", value_fmt="{:.3g}",
unit="", reference=None, reference_name="initial condition",
xlabel="x [m]", ylabel="", title="", ylim=None, height=430,
- ms_per_frame=90, colors=("#00A6D6", "#EC6842", "#7B61FF")):
+ ms_per_frame=90, colors=("#00A6D6", "#EC6842", "#7B61FF"),
+ log_x=False, log_y=False):
"""Profiles over time, with a play button and a draggable progress bar.
``frames`` is either one array of shape (len(values), nx), or a dict
@@ -74,8 +78,11 @@ def show_profiles(x, frames, values, *, value_name="t", value_fmt="{:.3g}",
fig.update_layout(
title=title, height=height, margin=dict(l=60, r=20, t=50, b=30),
- xaxis_title=xlabel, yaxis_title=ylabel,
- yaxis_range=[lo - pad, hi + pad], legend=dict(x=0.99, xanchor="right", y=0.99),
+ xaxis=dict(title=xlabel, type="log" if log_x else "linear"),
+ yaxis=dict(title=ylabel, type="log" if log_y else "linear",
+ range=([np.log10(max(lo, 1e-300)), np.log10(max(hi, 1e-299))]
+ if log_y else [lo - pad, hi + pad])),
+ legend=dict(x=0.99, xanchor="right", y=0.99),
sliders=[dict(active=0, steps=steps, len=0.82, x=0.18, y=0, pad=dict(t=50, b=8),
currentvalue=dict(prefix=f"{value_name} = ", suffix=unit,
font=dict(size=13)))],
@@ -182,6 +189,147 @@ def show_records(t, records, *, t_label="t [s]", y_label="T [K]", title="",
plt.show()
+def _arrow_lines(X, Y, U, V, length):
+ """Unit arrows at (X, Y), as one broken line trace. NaN separates them."""
+ n = np.hypot(U, V)
+ n = np.where(n > 0, n, 1.0)
+ u, v = U / n * length, V / n * length
+ x0, y0 = X - 0.5 * u, Y - 0.5 * v
+ x1, y1 = X + 0.5 * u, Y + 0.5 * v
+ ang = np.arctan2(v, u)
+ b = 0.34 * length
+ lx, ly = x1 - b * np.cos(ang - 0.42), y1 - b * np.sin(ang - 0.42)
+ rx, ry = x1 - b * np.cos(ang + 0.42), y1 - b * np.sin(ang + 0.42)
+ nan = np.full_like(x0, np.nan)
+ xs = np.stack([x0, x1, nan, lx, x1, rx, nan], axis=-1).ravel()
+ ys = np.stack([y0, y1, nan, ly, y1, ry, nan], axis=-1).ravel()
+ return xs, ys
+
+
+def show_map(x, y, frames, values, *, value_name="t", value_fmt="{:.3g}", unit="",
+ x_label="x [m]", y_label="y [m]", c_label="", title="", signed=None,
+ arrows=None, arrow_every=8, zmin=None, zmax=None, per_frame_scale=False,
+ log_c=False, log_floor=1e-4, y_down=False, equal_aspect=True,
+ height=560, ms_per_frame=140):
+ """A field over a plane, with a play button and a draggable slider.
+
+ ``frames`` has shape (len(values), len(x), len(y)): one map per slider
+ position, indexed the way ``np.meshgrid(..., indexing='ij')`` builds it.
+ ``arrows`` is an optional pair ``(U, V)`` of the same shape, drawn as unit
+ arrows every ``arrow_every`` points to show direction only.
+ ``per_frame_scale=True`` rescales the colour axis on every frame, which
+ shows the shape of a field whose amplitude changes by decades.
+ ``log_c=True`` colours by the base-10 logarithm instead, keeping the shape
+ of the pattern readable where the amplitude spans several decades; values
+ below ``log_floor`` times the largest are clipped.
+ """
+ x, y = np.asarray(x, float), np.asarray(y, float)
+ frames = np.asarray(frames, float)
+ if log_c:
+ if signed:
+ raise ValueError("log_c colours a magnitude, so the field must not change sign")
+ signed = False
+ scale = np.nanmax(np.abs(frames), axis=(1, 2), keepdims=True) if per_frame_scale \
+ else np.nanmax(np.abs(frames))
+ frames = np.log10(np.maximum(np.abs(frames) / scale, log_floor))
+ per_frame_scale = False
+ zmin, zmax = np.log10(log_floor), 0.0
+ if signed is None:
+ signed = np.nanmin(frames) < 0 < np.nanmax(frames)
+ if per_frame_scale:
+ lim = [float(np.nanmax(np.abs(f))) or 1.0 for f in frames]
+ rng = [(-L, L) if signed else (0.0, L) for L in lim]
+ else:
+ hi = float(np.nanmax(np.abs(frames))) if signed else float(np.nanmax(frames))
+ lo = -hi if signed else float(min(0.0, np.nanmin(frames)))
+ rng = [(lo if zmin is None else zmin, hi if zmax is None else zmax)] * len(values)
+ if log_c:
+ rng = [(zmin, zmax)] * len(values)
+ cs = "RdBu_r" if signed else "Inferno"
+
+ def heat(k):
+ return go.Heatmap(z=frames[k].T, x=x, y=y, colorscale=cs,
+ zmin=rng[k][0], zmax=rng[k][1], zsmooth="best",
+ colorbar=dict(title=c_label))
+
+ quiver = None
+ if arrows is not None:
+ U, V = (np.asarray(a, float) for a in arrows)
+ sx, sy = x[::arrow_every], y[::arrow_every]
+ Xa, Ya = np.meshgrid(sx, sy, indexing="ij")
+ step = float(min(sx[1] - sx[0], sy[1] - sy[0])) * 0.85
+ quiver = [_arrow_lines(Xa, Ya, U[k][::arrow_every, ::arrow_every],
+ V[k][::arrow_every, ::arrow_every], step)
+ for k in range(len(values))]
+
+ fig = go.Figure(heat(0))
+ if quiver is not None:
+ # cyan, not white: Inferno runs to near-white at the top, so white arrows
+ # vanish over the bright core exactly where the pattern is worth reading
+ fig.add_trace(go.Scatter(x=quiver[0][0], y=quiver[0][1], mode="lines",
+ line=dict(color="#00A6D6", width=1.4),
+ hoverinfo="skip", showlegend=False))
+ data_k = (lambda k: [heat(k), go.Scatter(x=quiver[k][0], y=quiver[k][1])]) if quiver \
+ else (lambda k: [heat(k)])
+ fig.frames = [go.Frame(name=str(k), data=data_k(k),
+ traces=list(range(len(fig.data)))) for k in range(len(values))]
+
+ still = dict(mode="immediate", frame=dict(duration=0, redraw=True),
+ transition=dict(duration=0))
+ steps = [dict(method="animate", label=_fmt(v, value_fmt), args=[[str(k)], still])
+ for k, v in enumerate(values)]
+ fig.update_layout(
+ title=title, height=height, margin=dict(l=60, r=20, t=50, b=30),
+ xaxis=dict(title=x_label, constrain="domain"),
+ yaxis=dict(title=y_label, autorange="reversed" if y_down else True,
+ scaleanchor="x" if equal_aspect else None, scaleratio=1),
+ sliders=[dict(active=0, steps=steps, len=0.80, x=0.18, y=0, pad=dict(t=50, b=8),
+ currentvalue=dict(prefix=f"{value_name} = ", suffix=unit,
+ font=dict(size=13)))],
+ updatemenus=[dict(type="buttons", direction="left", showactive=False,
+ x=0.0, y=0, xanchor="left", yanchor="top", pad=dict(t=55, b=8),
+ buttons=[
+ dict(label="\u25b6 play", method="animate",
+ args=[None, dict(mode="immediate", fromcurrent=False,
+ frame=dict(duration=ms_per_frame, redraw=True),
+ transition=dict(duration=0))]),
+ dict(label="\u275a\u275a", method="animate",
+ args=[[None], still])])])
+ fig.show()
+
+
+def show_isosurface(x, y, z, V, *, levels=(0.1, 0.3, 0.6), opacity=0.3, c_label="",
+ title="", wire=None, wire_name="source", x_label="x [m]",
+ y_label="y [m]", z_label="z [m]", z_down=True, height=600):
+ """Nested translucent surfaces of constant ``V``, free to rotate and zoom.
+
+ ``V`` has shape (len(x), len(y), len(z)), built with ``indexing='ij'``.
+ ``levels`` are fractions of the largest value in ``V``. Transparency keeps
+ the inner surfaces visible through the outer ones. ``wire`` is an optional
+ (xs, ys, zs) polyline, e.g. the loop that generated the field.
+ """
+ x, y, z = (np.asarray(a, float) for a in (x, y, z))
+ V = np.asarray(V, float)
+ X, Y, Z = np.meshgrid(x, y, z, indexing="ij")
+ vmax = float(np.nanmax(V))
+ lo, hi = sorted(float(f) * vmax for f in (min(levels), max(levels)))
+
+ fig = go.Figure(go.Isosurface(
+ x=X.ravel(), y=Y.ravel(), z=Z.ravel(), value=V.ravel(),
+ isomin=lo, isomax=hi, surface_count=len(levels), opacity=opacity,
+ colorscale="Inferno", caps=dict(x_show=False, y_show=False, z_show=False),
+ colorbar=dict(title=c_label)))
+ if wire is not None:
+ fig.add_trace(go.Scatter3d(x=wire[0], y=wire[1], z=wire[2], mode="lines",
+ line=dict(color="#00A6D6", width=6), name=wire_name))
+ fig.update_layout(
+ title=title, height=height, margin=dict(l=0, r=0, t=50, b=0),
+ scene=dict(xaxis_title=x_label, yaxis_title=y_label, zaxis_title=z_label,
+ zaxis=dict(autorange="reversed" if z_down else True),
+ aspectmode="data"))
+ fig.show()
+
+
# --------------------------------------------------------------------------
# Self-checks (same behaviour as the Lab 1 and Lab 2 module)
# --------------------------------------------------------------------------
diff --git a/book/3_diffusion_fields/labs/week05_diffusion.md b/book/3_diffusion_fields/labs/week05_diffusion.md
index 52a6600..7506c98 100644
--- a/book/3_diffusion_fields/labs/week05_diffusion.md
+++ b/book/3_diffusion_fields/labs/week05_diffusion.md
@@ -28,7 +28,7 @@ Everything you write here is plain NumPy. The module `fwtools` only draws and ch
- **Decaying modes.** Build an initial temperature from a Fourier series and explain why mode $n$ decays at $\kappa(n\pi/L)^2$, so fine structure goes first.
- **Boundaries.** Ends held at a fixed temperature let heat leave; insulated ends conserve it and drive the bar to its mean.
-- **Superposition.** Evolve any initial temperature with the Green's function, and turn the $t\propto x^2$ it implies into a measurement: of $\kappa$ in a bar, of depth in the Earth.
+- **Superposition.** Evolve any initial temperature with the Green's function, and turn the $t\propto x^2$ it implies into a measurement of $\kappa$ in a bar. The same kernel measures depth in the Earth in [Lab 6](./week06_em_diffusion.md).
---
@@ -41,7 +41,6 @@ import sys, pathlib
import numpy as np
import matplotlib.pyplot as plt
-from scipy.constants import mu_0, epsilon_0
# --- Live Code housekeeping, not part of the physics -----------------------
try:
@@ -639,150 +638,3 @@ So the two knobs answer two different physical questions, one about the far fiel
one about the source at early time. This is Task 6 read backwards. There we asked when a real
boundary makes itself felt; here we ask how far away to put a fictitious one so that it never does.
:::
-
----
-
-## Part 4 — The same kernel in the Earth
-
-A horizontal current sheet is switched on at $t=0$ with current $I$ per unit width. In the diffusive approximation the electric field a distance $z$ from the sheet is
-
-$$ E_x(z,t) = -I\sqrt{\frac{\mu}{4\pi\sigma t}}\exp\!\left(-\frac{\sigma\mu z^2}{4t}\right), \qquad t>0 . $$
-
-The exponential is the kernel of Part 2 with $\kappa$ replaced by $1/(\sigma\mu)$, so the field soaks into the ground as heat spreads along a wire. The prefactors differ by the constant $I/\sigma$, which sets the amplitude and not the shape. Rock is non-magnetic, so $\mu=\mu_0$.
-
-The sheet sits in a uniform conductor and the solution is even in $z$, so reading $z$ as depth means ignoring the air above: an idealisation, not a consequence of the derivation. The approximation also neglects $\varepsilon\,\partial_t\boldsymbol{E}$, which is justified when $t$ is much longer than the charge relaxation time $\tau_r = \varepsilon/\sigma$. This is the one-dimensional idealisation behind time-domain electromagnetic sounding; a real system uses a finite loop and records $\partial_t B$.
-
-### Task 8 — when does the field arrive, and in what ground?
-
-Write the time at which $\lvert E_x\rvert$ at depth $z$ peaks, from $\partial_t\lvert E_x\rvert = 0$. Differentiating $\ln\lvert E_x\rvert$ is easier and has the same roots. **Predict** first whether a deeper receiver peaks earlier or later, and by what power of $z$.
-
-Then the same receiver in three grounds: seawater at 4 S/m, saturated sediments at $10^{-2}$ S/m, crystalline basement at $10^{-4}$ S/m. **Predict**, before running, how they rank by arrival time at 300 m and by what factor $t_p$ changes from basement to seawater.
-
-The second blank is the margin the diffusive approximation has left, $t_p/\tau_r$ with $\tau_r=\varepsilon/\sigma$. Write it in closed form, in terms of $z$ and $\sigma$, and watch what $\sigma$ does to it.
-
-```{code-cell} ipython3
-sigma, eps_r, I_s = 0.01, 10.0, 1.0 # sediments [S/m], relative permittivity, current [A/m]
-tau_r = eps_r * epsilon_0 / sigma # charge relaxation time [s]
-t_em = np.logspace(-6, -1, 3001) # [s], from about 100 relaxation times
-_depths = [100.0, 300.0, 1000.0] # [m]
-
-def E_step(z, t, s=sigma):
- return -I_s * np.sqrt(mu_0 / (4*np.pi*s*t)) * np.exp(-s*mu_0*z**2 / (4*t))
-
-# Task 8 -- two blanks. `s` is the conductivity [S/m]; it defaults to the sediments above.
-def t_peak_em(z, s=sigma):
- """Time at which |E_x| at depth z peaks [s]."""
- return ___
-
-def peak_over_relax(z, s):
- """How many charge relaxation times the peak waits: t_p / tau_r. Dimensionless."""
- return ___
-
-# --- given, 1: the field soaking downward. Drag the bar, or press play ---
-_zs = np.linspace(1.0, 2000.0, 400) # depth [m]
-_tp = np.logspace(-6, -1.3, 60) # [s]
-_prof = np.array([np.abs(E_step(_zs, tv)) / np.abs(E_step(1.0, tv)) for tv in _tp])
-fw.show_profiles(_zs, _prof, _tp, value_name="t", value_fmt="{:.2g}", unit=" s",
- xlabel="depth z [m]", ylabel="|Ex| / its value at the surface",
- title="The field diffuses into the Earth", ylim=(0.0, 1.05))
-
-# --- given, 2: depth against time, with your arrival-time curve over it ---
-fw.show_spacetime(_zs, _tp, _prof, t_label="t [s]", x_label="depth z [m]",
- c_label=r"$|E_x|$ / surface value", signed=False,
- curve=(t_peak_em(_zs), _zs), curve_label=r"your $t_p(z)$", deeper_down=True,
- title=r"Deeper means later: $t_p \propto z^2$")
-
-# --- given, 3: what a receiver records at three depths ---
-fw.show_records(t_em, {f"z = {z:g} m": np.abs(E_step(z, t_em)) / np.max(np.abs(E_step(z, t_em)))
- for z in _depths},
- peaks={f"z = {z:g} m": (t_peak_em(z), 1.0) for z in _depths},
- t_label="t [s]", y_label=r"$|E_x|$ / peak value",
- title=r"Step response at depth, $\sigma$ = 0.01 S/m")
-for z in _depths:
- print(f"z = {z:6.0f} m: t_p = {t_peak_em(z):.3g} s, t_p / tau_r = {t_peak_em(z)/tau_r:.2g},"
- f" |E| at the peak = {abs(E_step(z, t_peak_em(z))):.3g} V/m")
-
-# --- self-check (leave this alone) ---
-for z in _depths:
- fw.check_scalar(f"t_p at {z:g} m against the peak of the record", t_peak_em(z),
- t_em[np.argmax(np.abs(E_step(z, t_em)))], rtol=0.01, unit=" s")
-fw.check_close("t_p / tau_r against t_p and tau_r built separately, over five decades of sigma",
- [peak_over_relax(300.0, sv) for sv in (1e-4, 1e-3, 1e-2, 1e-1, 1.0)],
- [t_peak_em(300.0, sv) / (eps_r * epsilon_0 / sv)
- for sv in (1e-4, 1e-3, 1e-2, 1e-1, 1.0)], rtol=1e-6,
- hint="tau_r = eps_r eps_0 / sigma shrinks as t_p grows, so the ratio carries sigma twice")
-fw.check("a tenfold drop in conductivity costs a hundredfold in t_p / tau_r",
- np.isclose(peak_over_relax(300.0, 1e-3) / peak_over_relax(300.0, 1e-4), 100.0, rtol=1e-9),
- "the margin goes as sigma^2, not as sigma")
-```
-
-```{code-cell} ipython3
-# --- given: one receiver at 300 m, three grounds -----------------------------
-media = {"seawater, 4 S/m": 4.0,
- "saturated sediments, 1e-2 S/m": 0.01,
- "crystalline basement, 1e-4 S/m": 1e-4}
-z_s = 300.0 # the receiver depth [m]
-t_s = np.logspace(-8, 1, 900) # [s]
-
-fw.show_records(t_s, {k: np.abs(E_step(z_s, t_s, sv)) / np.max(np.abs(E_step(z_s, t_s, sv)))
- for k, sv in media.items()},
- peaks={k: (t_peak_em(z_s, sv), 1.0) for k, sv in media.items()},
- t_label="t [s]", y_label="|E| / peak value",
- title=f"The same receiver at {z_s:.0f} m, in three grounds")
-
-print(f"at z = {z_s:.0f} m\n")
-print(f"{'ground':<30}{'t_p [s]':>11}{'tau_r [s]':>11}{'t_p / tau_r':>13}")
-for k, sv in media.items():
- print(f"{k:<30}{t_peak_em(z_s, sv):>11.3g}{eps_r*epsilon_0/sv:>11.3g}"
- f"{peak_over_relax(z_s, sv):>13.3g}")
-
-print("\na survey looking 500 m down, able to record for 10 ms:")
-for k, sv in media.items():
- _tp = t_peak_em(500.0, sv)
- print(f" {k:<30} t_p = {_tp:.3g} s {'inside' if _tp < 1e-2 else 'beyond'} the window")
-
-# --- self-check (leave this alone) ---
-_tf = np.logspace(-8, 1, 20001)
-for k, sv in media.items():
- fw.check_scalar(f"t_p at {z_s:.0f} m in {k.split(',')[0]}", t_peak_em(z_s, sv),
- _tf[np.argmax(np.abs(E_step(z_s, _tf, sv)))], rtol=0.01, unit=" s")
-fw.check("basement to seawater is a factor 4e4 in arrival time",
- np.isclose(t_peak_em(z_s, 4.0) / t_peak_em(z_s, 1e-4), 4e4, rtol=1e-9),
- "t_p is linear in sigma, and the two conductivities differ by 4e4")
-```
-
-:::{admonition} Solution — Task 8
-:class: dropdown
-
-```python
-def t_peak_em(z, s=sigma):
- return s * mu_0 * z**2 / 2
-
-def peak_over_relax(z, s):
- return s**2 * mu_0 * z**2 / (2 * eps_r * epsilon_0)
-```
-:::
-
-:::{admonition} Depth maps onto arrival time, as $t_p\propto z^2$
-:class: important dropdown
-
-A layer ten times deeper answers a hundred times later: 63 µs at 100 m, 6.3 ms at 1 km, for $\sigma=0.01$ S/m. On the map your curve $t_p(z)$ runs along the shoulder of the diffusing field, because $|E_x|$ has fallen to $1/\sqrt e$ of its surface value exactly where $t = \sigma\mu z^2/2$. Time-domain electromagnetic soundings rest on this: the later part of a record senses the deeper Earth. In a more conductive Earth the arrival is later, since $t_p\propto\sigma$. The ratio $t_p/\tau_r$ is 7000 at 100 m and grows as $z^2$, so neglecting $\varepsilon\,\partial_t\boldsymbol{E}$ is well justified at these depths. The blurring of Task 5 applies here too: deep structure arrives late, and therefore blurred.
-
-The peak amplitude is 0.24 V/m at 100 m and 0.024 V/m at 1 km. A factor of ten in depth costs a factor of ten in signal and a factor of a hundred in time, which is the $1/\lvert x\rvert$ law of Task 4 with depth in place of distance along the wire.
-
-The chapter's three-dimensional impulse response peaks at $t_p=\sigma\mu r^2/6$, three times earlier than the $\sigma\mu z^2/2$ you found here, and its amplitude falls as $1/r^3$. Both follow from the same differentiation: a prefactor $t^{-p}$ gives $t_p=\sigma\mu r^2/(4p)$, with $p=1/2$ for the sheet and $p=3/2$ for the point.
-:::
-
-:::{admonition} Conductivity sets the clock, and it sets the margin twice over
-:class: important dropdown
-
-At 300 m the peak arrives after 0.23 s in seawater, 0.57 ms in sediments and 5.7 µs in basement. That is a factor $4\times10^{4}$ from basement to seawater, exactly the ratio of the two conductivities, because $t_p=\sigma\mu z^2/2$ is linear in $\sigma$. A conductive ground holds the field up; a resistive one lets it through almost at once.
-
-The margin behaves differently. $\tau_r=\varepsilon/\sigma$ grows as the ground becomes resistive while $t_p$ shrinks, so
-
-$$ \frac{t_p}{\tau_r} = \frac{\sigma^2\mu z^2}{2\varepsilon} $$
-
-carries $\sigma$ twice and falls a hundredfold for every tenfold drop in conductivity: $10^{10}$ in seawater, $6.4\times10^{4}$ in sediments, and **6.4 in basement**. The diffusive approximation is comfortable in the first two and running out in the third, and it is the *resistive* ground that breaks it, not the conductive one. Dropping $\varepsilon\,\partial_t\boldsymbol{E}$ was never a statement about depth. It is a statement about how fast charge can relax, and in a resistive ground it cannot relax quickly.
-
-The survey question resolves the same way from the other end. Looking 500 m down with a 10 ms record: seawater peaks at 0.63 s, far beyond the window, so the target is never reached in time. Sediments peak at 1.6 ms, comfortably inside. Basement peaks at 16 µs, inside the window but so early that a real transmitter is still switching off, and the useful signal is buried in the turn-off transient. Conductivity, not instrument quality, decides which depths a time-domain sounding can reach.
-:::
diff --git a/book/3_diffusion_fields/labs/week06_em_diffusion.md b/book/3_diffusion_fields/labs/week06_em_diffusion.md
new file mode 100644
index 0000000..97d650a
--- /dev/null
+++ b/book/3_diffusion_fields/labs/week06_em_diffusion.md
@@ -0,0 +1,799 @@
+---
+jupytext:
+ text_representation:
+ extension: .md
+ format_name: myst
+ format_version: 0.13
+kernelspec:
+ display_name: Python 3 (ipykernel)
+ language: python
+ name: python3
+mystnb:
+ # Workbook page: the task cells contain `___` blanks by design, so it must
+ # not be executed at build time. Readers run it themselves with Live Code.
+ execution_mode: 'off'
+---
+
+# Lab 6: Diffusive electromagnetic fields
+
+:::{admonition} Computer lab
+:class: note
+
+A practical companion to the lecture on [electromagnetic diffusive field equations with a source](../em_diffusion.md). Each task states a physical question, gives the formula it rests on, and ends with a picture you can read the answer off.
+
+Everything you write is plain NumPy, and every line that implements an equation names it in a comment. The module `fwtools` only draws: `fw.show_map` (a field over a plane, with a play button), `fw.show_isosurface` (a rotatable 3-D view), `fw.show_profiles`, `fw.show_spacetime` and `fw.show_records`. Every slider runs in the browser.
+:::
+
+## Learning objectives
+
+- **One kernel, three geometries.** A sheet, a wire and a point all produce $t^{-p}\exp(-\sigma\mu r^2/4t)$. Only $p$ changes, and it fixes when the field arrives. A loop is a superposition of these along four sides, so it has no single kernel, but its arrival still obeys the same $t_p$ with $p=5/2$.
+- **Arrival time measures distance.** $t_p=\sigma\mu r^2/4p$ turns a recorded time into a depth, and the conductivity of the ground sets the clock.
+- **The magnetic field becomes a potential field.** $(\nabla\times\boldsymbol H)_z=\sigma E_z$ while the field diffuses, and both sides die together.
+- **A source forgets its shape.** Once the diffusion distance passes the size of the loop, the field no longer knows the loop was square.
+
+---
+
+## Part 0 — Setup
+
+Run this once. It contains no physics: it fetches two packages the browser lacks and locates `fwtools`.
+
+```{code-cell} ipython3
+import sys, pathlib
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.constants import mu_0, epsilon_0
+from scipy.special import erf
+
+# --- Live Code housekeeping, not part of the physics -----------------------
+try:
+ import plotly.io as pio
+except ModuleNotFoundError:
+ print("Fetching plotly. A few seconds, and only the first time...")
+ import micropip
+ await micropip.install("plotly")
+ import plotly.io as pio
+
+try:
+ import nbformat # noqa: F401
+except ModuleNotFoundError:
+ import micropip, types
+ try:
+ await micropip.install("nbformat")
+ except Exception:
+ _nb = types.ModuleType("nbformat")
+ _nb.__version__ = "5.10.4"
+ sys.modules["nbformat"] = _nb
+
+for _p in (".", "book/3_diffusion_fields/labs"):
+ if (pathlib.Path(_p) / "fwtools.py").exists():
+ sys.path.insert(0, _p)
+ break
+try:
+ import fwtools as fw
+except ModuleNotFoundError:
+ from pyodide.http import pyfetch
+ _r = await pyfetch("fwtools.py")
+ pathlib.Path("fwtools.py").write_bytes(await _r.bytes())
+ import fwtools as fw
+
+pio.renderers.default = "plotly_mimetype+notebook"
+mu = mu_0 # rock and water are non-magnetic
+# ---------------------------------------------------------------------------
+print("ready")
+```
+
+---
+
+## Part 1 — One dimension: a current sheet
+
+A horizontal current sheet carrying $I$ per unit width is switched on at $t=0$. The field a distance $z$ from it is {eq}`eq:ExtH`,
+
+$$ E_x(z,t) = -I\sqrt{\frac{\mu}{4\pi\sigma t}}\exp\left(-\frac{\sigma\mu z^2}{4t}\right),\qquad t>0 . $$
+
+The exponential is the kernel of the heat equation with $\kappa$ replaced by $1/(\sigma\mu)$, so the field soaks into the ground as heat spreads along a wire. The prefactor falls as $t^{-1/2}$.
+
+Every field in this lab has that form, $t^{-p}\exp(-\sigma\mu r^2/4t)$, with a different $p$. It is worth finding the peak time once, for any $p$.
+
+```{figure} figures/setup_sheet.svg
+:name: fig-lab6-sheet
+:width: 85%
+
+The model of this part. An infinite sheet of current switched on at $t=0$, in a uniform conductor, with receivers at three depths.
+```
+
+### Task 1 — when does the field arrive, and in what ground?
+
+Setting $\partial_t\ln|E|=0$ is easier than differentiating $|E|$ itself and gives the same answer. With a prefactor $t^{-p}$,
+
+$$ -\frac{p}{t} + \frac{\sigma\mu r^2}{4t^2} = 0 . $$
+
+Solve it for $t$ and write the result once, as `t_peak(r, p, s)`.
+
+**Predict** before running, and sketch it: on one logarithmic time axis, draw $|E_x|(t)$ at 100 m and at 1000 m, each scaled to its own peak. Are the two curves the same shape, or does the deeper one spread out? Shifted by what factor? Then rank the same receiver in three grounds, seawater at 4 S/m, saturated sediments at $10^{-2}$ S/m and crystalline basement at $10^{-4}$ S/m, by arrival time.
+
+The second blank is the margin the diffusive approximation still has. Dropping $\varepsilon\,\partial_t\boldsymbol E$ is safe only while the peak arrives long after the charge relaxation time $\tau_r=\varepsilon/\sigma$. Write $t_p/\tau_r$ for the sheet in closed form, in terms of $z$ and $\sigma$, and watch what $\sigma$ does to it.
+
+Dividing by the value at the shallowest depth cancels the common $t^{-1/2}$ and leaves the bare kernel $\exp(-\sigma\mu z^2/4t)$: how far the field has reached, not how strong it is. That ratio only climbs and saturates, so look for the shoulder at $1/\sqrt e$, not a maximum. A receiver keeps the $t^{-1/2}$, and there the arriving exponential and the decaying prefactor compete. That competition is what $t_p$ solves.
+
+```{code-cell} ipython3
+sigma_sheet, eps_r, I_sheet = 0.01, 10.0, 1.0 # sediments [S/m], relative permittivity, [A/m]
+tau_r = eps_r * epsilon_0 / sigma_sheet # charge relaxation time [s]
+
+def E_sheet(z, t, s=sigma_sheet):
+ """eq:ExtH -- the field of a current sheet switched on at t = 0 [V/m]."""
+ return -I_sheet * np.sqrt(mu / (4*np.pi*s*t)) * np.exp(-s*mu*z**2 / (4*t))
+
+def t_for_D(D, s):
+ """When the diffusion distance sqrt(4t/(s*mu)) reaches D [s]. Used only to
+ frame the figures, so that every axis follows the parameters you set."""
+ return s * mu * D**2 / 4
+
+# Task 1 -- two blanks.
+def t_peak(r, p, s):
+ """Time at which t**(-p) * exp(-s*mu*r**2/(4t)) is largest [s].
+ Solve -p/t + s*mu*r**2/(4*t**2) = 0 for t."""
+ return ___
+
+def peak_over_relax(z, s):
+ """How many charge relaxation times the sheet's peak waits, t_p / tau_r.
+ The sheet has p = 1/2, and tau_r = eps_r*epsilon_0/s."""
+ return ___
+
+# --- given: the field soaking downward. Drag the bar, or press play ---------
+z_grid = np.linspace(1.0, 2000.0, 400) # depth [m]
+t_grid = np.logspace(np.log10(t_for_D(z_grid[-1], sigma_sheet)) - 4,
+ np.log10(t_for_D(z_grid[-1], sigma_sheet)), 60) # [s]
+soak = np.array([np.abs(E_sheet(z_grid, tv)) / np.abs(E_sheet(z_grid[0], tv))
+ for tv in t_grid])
+fw.show_profiles(z_grid, soak, t_grid, value_name="t", value_fmt="{:.2g}", unit=" s", log_x=True,
+ xlabel="depth z [m]", ylabel="|Ex| / its value at the surface",
+ title="The field diffuses into the Earth", ylim=(0.0, 1.05))
+
+# --- given: depth against time, with your arrival-time curve over it --------
+fw.show_spacetime(z_grid, t_grid, soak, t_label="t [s]", x_label="depth z [m]",
+ c_label=r"$|E_x|$ / surface value", signed=False,
+ curve=(t_peak(z_grid, 0.5, sigma_sheet), z_grid),
+ curve_label=r"your $t_p(z)$", deeper_down=True,
+ title=r"$|E_x|$ against depth and time, with your $t_p(z)$ drawn over it")
+
+# --- given: what a receiver records at three depths ------------------------
+depths = [100.0, 300.0, 1000.0] # [m]
+t_rec = np.logspace(np.log10(t_for_D(min(depths), sigma_sheet)) - 2,
+ np.log10(t_for_D(max(depths), sigma_sheet)) + 2, 3001) # [s]
+fw.show_records(t_rec, {f"z = {z:g} m": np.abs(E_sheet(z, t_rec)) / np.max(np.abs(E_sheet(z, t_rec)))
+ for z in depths},
+ peaks={f"z = {z:g} m": (t_peak(z, 0.5, sigma_sheet), 1.0) for z in depths},
+ t_label="t [s]", y_label=r"$|E_x|$ / peak value",
+ title=f"Step response at depth, sigma = {sigma_sheet:g} S/m")
+
+_zc = depths[1]
+fw.check_scalar(f"t_p at {_zc:g} m against the peak of the record",
+ t_peak(_zc, 0.5, sigma_sheet),
+ t_rec[np.argmax(np.abs(E_sheet(_zc, t_rec)))], rtol=0.01, unit=" s")
+```
+
+The same expression evaluated at three conductivities. $t_p$ is linear in $\sigma$ while $\tau_r=\varepsilon/\sigma$ is inverse in it, so those two move opposite ways and their ratio, the margin, carries $\sigma$ twice. A margin above about 100 is comfortable.
+
+```{code-cell} ipython3
+# --- given: one receiver at 300 m, three grounds ---------------------------
+media = {"seawater, 4 S/m": 4.0,
+ "saturated sediments, 1e-2 S/m": 0.01,
+ "crystalline basement, 1e-4 S/m": 1e-4}
+z_s = 300.0 # the receiver depth [m]
+_tp_all = [t_peak(z_s, 0.5, sv) for sv in media.values()]
+t_s = np.logspace(np.log10(min(_tp_all)) - 3, np.log10(max(_tp_all)) + 2, 900)
+
+fw.show_records(t_s, {k: np.abs(E_sheet(z_s, t_s, sv)) / np.max(np.abs(E_sheet(z_s, t_s, sv)))
+ for k, sv in media.items()},
+ peaks={k: (t_peak(z_s, 0.5, sv), 1.0) for k, sv in media.items()},
+ t_label="t [s]", y_label="|E| / peak value",
+ title=f"The same receiver at {z_s:.0f} m, in three grounds")
+
+print(f"at z = {z_s:.0f} m\n")
+print(f"{'ground':<30}{'t_p [s]':>11}{'tau_r [s]':>11}{'t_p / tau_r':>13}")
+for k, sv in media.items():
+ print(f"{k:<30}{t_peak(z_s, 0.5, sv):>11.3g}{eps_r*epsilon_0/sv:>11.3g}"
+ f"{peak_over_relax(z_s, sv):>13.3g}")
+
+fw.check("a tenfold drop in conductivity costs a hundredfold in margin",
+ np.isclose(peak_over_relax(z_s, 1e-3) / peak_over_relax(z_s, 1e-4), 100.0, rtol=1e-9),
+ "the margin goes as sigma squared, not as sigma")
+fw.check_scalar("the margin itself, against t_p and tau_r worked out separately",
+ peak_over_relax(z_s, sigma_sheet),
+ t_peak(z_s, 0.5, sigma_sheet) / (eps_r*epsilon_0/sigma_sheet), rtol=1e-6)
+```
+
+:::{admonition} Solution — Task 1
+:class: dropdown
+
+```python
+def t_peak(r, p, s):
+ return s * mu * r**2 / (4 * p)
+
+def peak_over_relax(z, s):
+ return s**2 * mu * z**2 / (2 * eps_r * epsilon_0)
+```
+:::
+
+:::{admonition} Depth maps onto arrival time, and conductivity sets the clock
+:class: important dropdown
+
+A layer ten times deeper answers a hundred times later: 63 µs at 100 m, 6.3 ms at 1 km, for $\sigma=0.01$ S/m. Your curve runs along the shoulder of the diffusing field, because $|E_x|$ has fallen to $1/\sqrt e$ of its surface value exactly where $t=\sigma\mu z^2/2$. Time-domain electromagnetic soundings rest on this: the later part of a record senses the deeper Earth.
+
+At 300 m the peak arrives after 0.23 s in seawater, 0.57 ms in sediments and 5.7 µs in basement, a factor $4\times10^4$ from basement to seawater. That is exactly the ratio of the two conductivities, because $t_p=\sigma\mu z^2/2$ is linear in $\sigma$. A conductive ground holds the field up; a resistive one lets it through at once.
+
+The margin behaves differently. $\tau_r=\varepsilon/\sigma$ grows as the ground becomes resistive while $t_p$ shrinks, so
+
+$$ \frac{t_p}{\tau_r} = \frac{\sigma^2\mu z^2}{2\varepsilon} $$
+
+carries $\sigma$ twice: $10^{10}$ in seawater, $6.4\times10^4$ in sediments, and **6.4** in basement. The diffusive approximation is comfortable in the first two and running out in the third, and it is the *resistive* ground that breaks it. Dropping $\varepsilon\,\partial_t\boldsymbol E$ was never a statement about depth. It is a statement about how fast charge can relax.
+:::
+
+---
+
+## Part 2 — Two dimensions: an infinite wire
+
+A long wire along $z$ carries a current switched on at $t=0$. Nothing depends on $z$, so the field lives in the $(x,y)$-plane: the electric field is $E_z$, perpendicular to that plane, and the magnetic field lies in it. This is the Transverse Electric mode of the lecture, and the chapter gives
+
+$$ E_z(\varrho,t) = -\frac{\mu I}{4\pi t}\exp\left(-\frac{\sigma\mu\varrho^2}{4t}\right),\qquad \varrho=\sqrt{x^2+y^2}. $$
+
+Negative, as the sheet's field was: the current is switched on, and the field it drives in the ground opposes it.
+
+The kernel is the same as the sheet's. The prefactor is now $t^{-1}$, so $p=1$.
+
+```{figure} figures/setup_wire.svg
+:name: fig-lab6-wire
+:width: 92%
+
+The model of this part. The wire is a point in the plane you will plot, the electric field points out of that plane, and the magnetic field circles the wire within it.
+```
+
+### Task 2 — the field of the wire, and when it arrives
+
+Write $E_z$. Then **predict**, before running: at the same distance, does the wire's field peak earlier or later than the sheet's, and by what factor?
+
+The colour is logarithmic and rescaled frame by frame. The kernel falls off so steeply that one linear scale would leave a bright spot and an empty plane, and the level drops as $1/t$ across the run, so a fixed scale would lose the late frames.
+
+```{code-cell} ipython3
+sigma_wire, I_wire = 3.0, 1.0 # seawater [S/m], wire current [A]
+
+# Task 2 -- one blank.
+def E_wire(rho, t, s=sigma_wire):
+ """E_z of an infinitely long wire with a step switch-on current [V/m].
+ Prefactor t**(-1), times the kernel exp(-s*mu*rho**2/(4t)). Opposes the
+ current, so it is negative."""
+ return ___
+
+# --- given: the field spreading out from the wire, on a logarithmic colour scale ---
+n_w, half_w = 91, 50.0
+x_w = np.linspace(-half_w, half_w, n_w) # [m]
+X_w, Y_w = np.meshgrid(x_w, x_w, indexing="ij") # F[ix, iy] sits at (x[ix], y[iy])
+rho_w = np.hypot(X_w, Y_w)
+rho_w[rho_w == 0] = x_w[1] - x_w[0] # keep the wire out of the formula
+# from a fifth of the box to twice it: past that the picture is uniform colour
+t_w = np.logspace(np.log10(t_for_D(half_w / 5, sigma_wire)),
+ np.log10(t_for_D(half_w * 2, sigma_wire)), 16) # [s]
+
+E_frames = np.array([E_wire(rho_w, tv) for tv in t_w])
+fw.show_map(x_w, x_w, E_frames, t_w, value_name="t", value_fmt="{:.2g}", unit=" s",
+ x_label="x [m]", y_label="y [m]", c_label="log10(Ez / max)",
+ log_c=True, per_frame_scale=True, title="The electric field spreads out from the wire")
+
+# --- given: three receivers, with your arrival times marked ----------------
+dists = [0.36, 15.0, 30.0] # [m]
+d_mid = dists[1]
+t_r = np.logspace(np.log10(t_for_D(min(dists), sigma_wire)) - 2,
+ np.log10(t_for_D(max(dists), sigma_wire)) + 3, 3001) # [s]
+fw.show_records(t_r, {f"{d:g} m": np.abs(E_wire(d, t_r)) / np.max(np.abs(E_wire(d, t_r))) for d in dists},
+ peaks={f"{d:g} m": (t_peak(d, 1, sigma_wire), 1.0) for d in dists},
+ t_label="t [s]", y_label=r"$|E_z|$ / peak value",
+ title=f"The wire, sigma = {sigma_wire:g} S/m:"
+ f" close in, the field is already fading")
+
+print(f"the sheet and the wire at the same {d_mid:g} m, in the same ground:")
+print(f" sheet, p = 1/2 : t_p = {t_peak(d_mid, 0.5, sigma_wire)*1e3:.3f} ms")
+print(f" wire, p = 1 : t_p = {t_peak(d_mid, 1.0, sigma_wire)*1e3:.3f} ms")
+
+_tpm = t_peak(d_mid, 1, sigma_wire)
+fw.check_scalar(f"t_p at {d_mid:g} m against the peak of the record", _tpm,
+ t_r[np.argmax(np.abs(E_wire(d_mid, t_r)))], rtol=0.01, unit=" s")
+fw.check_scalar("the peak value, against -mu I / (4 pi t_p) times 1/e", E_wire(d_mid, _tpm),
+ -mu*I_wire/(4*np.pi*_tpm)*np.exp(-1.0), rtol=0.01, unit=" V/m")
+```
+
+:::{admonition} Solution — Task 2
+:class: dropdown
+
+```python
+def E_wire(rho, t, s=sigma_wire):
+ return -mu * I_wire / (4*np.pi*t) * np.exp(-s*mu*rho**2 / (4*t))
+```
+:::
+
+### Task 3 — the magnetic field climbs to its static value
+
+The wire's magnetic field circles it. Long after switch-on it must settle at the static result of the curl chapter, $\boldsymbol H = I\hat{\boldsymbol\varphi}/(2\pi\varrho)$. The chapter shows that it approaches that value through the same exponential that shapes $E_z$, with no power of $t$ in front:
+
+$$ \frac{|\boldsymbol H|}{|\boldsymbol H|_{\text{static}}} = \exp\left(-\frac{\sigma\mu\varrho^2}{4t}\right),\qquad |\boldsymbol H|_{\text{static}} = \frac{I}{2\pi\varrho}. $$
+
+Write that ratio. It runs from 0 at switch-on to 1 at late time, so the picture shows how far the field has got, not how strong it is.
+
+Then move the conductivity slider on the second figure. **Predict** first: in which ground does the magnetic field reach its static value soonest?
+
+```{code-cell} ipython3
+# Task 3 -- one blank.
+def H_over_static(rho, t, s=sigma_wire):
+ """|H| divided by its static value. Dimensionless, 0 at switch-on and 1 at late time."""
+ return ___
+
+# --- given: the magnetic field filling in, with its direction ---------------
+H_frames = np.array([H_over_static(rho_w, tv) for tv in t_w])
+Hx_dir = np.broadcast_to(-Y_w / rho_w, (len(t_w),) + rho_w.shape) # the unit vector phi-hat
+Hy_dir = np.broadcast_to( X_w / rho_w, (len(t_w),) + rho_w.shape)
+fw.show_map(x_w, x_w, H_frames, t_w, arrows=(Hx_dir, Hy_dir), arrow_every=10,
+ value_name="t", value_fmt="{:.2g}", unit=" s", zmin=0.0, zmax=1.0,
+ x_label="x [m]", y_label="y [m]", c_label="|H| / static value",
+ title="The magnetic field is already circular. It is only filling in")
+
+# --- given: how far it still has to go, over five conductivities -----------
+sigmas = np.array([0.01, 0.1, 1.0, 3.0, 10.0]) # [S/m]
+# wide enough to show the 1/t tail for every curve, and no wider: far below this
+# the exponential is zero to machine precision and nothing is left to plot
+_a_lo = t_for_D(min(dists[1:]), sigmas.min())
+_a_hi = t_for_D(max(dists[1:]), sigmas.max())
+t_h = np.logspace(np.log10(_a_hi / 500), np.log10(_a_lo * 1e7), 400) # [s]
+gap = {f"{d:g} m": np.array([1.0 - H_over_static(d, t_h, sv) for sv in sigmas])
+ for d in dists[1:]}
+fw.show_profiles(t_h, gap, sigmas, value_name="σ", value_fmt="{:.3g}", unit=" S/m",
+ xlabel="t [s]", ylabel="1 − |H| / static value", log_x=True, log_y=True,
+ title="How far the magnetic field still has to go, against time")
+
+# "early" and "late" written in terms of the problem, so the test survives a change
+# of sigma_wire or of the distance: the field has spread a fifth of d_mid, then a hundred times it.
+_early, _late = t_for_D(d_mid / 5, sigma_wire), t_for_D(d_mid * 100, sigma_wire)
+fw.check("the ratio climbs from 0 at switch-on to 1 at late time",
+ H_over_static(d_mid, _early) < 1e-6 and abs(H_over_static(d_mid, _late) - 1.0) < 1e-3,
+ "it is exp(-a/t), not 1 - exp(-a/t): at t -> 0 there is no field yet")
+# the shape alone does not pin the clock: a wrong factor in the exponent still
+# climbs from 0 to 1. This fixes when, using the t_p you already wrote.
+fw.check_scalar("|H| is at 1/e of static exactly at the wire's own peak time",
+ H_over_static(d_mid, t_peak(d_mid, 1, sigma_wire)), np.exp(-1.0), rtol=1e-6)
+```
+
+:::{admonition} Solution — Task 3
+:class: dropdown
+
+```python
+def H_over_static(rho, t, s=sigma_wire):
+ return np.exp(-s*mu*rho**2 / (4*t))
+```
+:::
+
+### Task 4 — while it diffuses, the magnetic field has curl
+
+The chapter's two-dimensional equation that connects the fields to the source is $-\partial_xH_y+\partial_yH_x+\sigma E_z=0$ away from the wire, which is Ampère's law, $\nabla\times\boldsymbol H=\sigma\boldsymbol E$, in two dimensions. The curl of $\boldsymbol H$ is not an abstraction here: it is the electric field, times $\sigma$. As $E_z$ dies the curl dies with it, and the magnetic field becomes a potential field.
+
+Write the $z$-component of the curl. The two derivatives of each component are already computed for you.
+
+With no source this says $\partial_xH_y-\partial_yH_x=\sigma E_z$, and both sides are negative because $E_z$ is. Where the two pictures agree, the diffusing magnetic field is demonstrably not curl-free. Near the wire they will not agree, because a finite grid cannot follow $1/\varrho$.
+
+```{code-cell} ipython3
+# One instant, chosen so the field has spread about as far as the ring the check uses.
+# That makes the accuracy below independent of sigma_wire.
+r_ring = (0.4*half_w, 0.9*half_w) # the comparison ring [m]
+t_c = t_for_D(np.mean(r_ring), sigma_wire) # [s]
+dx = x_w[1] - x_w[0]
+H_mag = I_wire / (2*np.pi*rho_w) * H_over_static(rho_w, t_c)
+Hx_c, Hy_c = -Y_w/rho_w * H_mag, X_w/rho_w * H_mag # [A/m]
+
+dHx_dx, dHx_dy = np.gradient(Hx_c, dx, dx) # indexing='ij', so [0] is d/dx
+dHy_dx, dHy_dy = np.gradient(Hy_c, dx, dx)
+
+# Task 4 -- one blank: the z-component of the curl of H.
+curl_z = ___
+
+# --- given: the curl beside the electric field it equals -------------------
+shown = rho_w > 7*dx # H goes as 1/rho, so the grid cannot follow it at the wire
+fw.show_map(x_w, x_w, np.array([np.where(shown, curl_z, np.nan),
+ np.where(shown, sigma_wire*E_wire(rho_w, t_c), np.nan)]),
+ ["curl of H", "sigma times E_z"], value_name="showing", value_fmt="{}",
+ x_label="x [m]", y_label="y [m]", c_label="A/m²",
+ title=f"Two ways to the same picture, at t = {t_c*1e3:g} ms")
+
+# The comparison is made on the outer ring, where the grid resolves 1/rho. Push t_c
+# up and it will fail: the true curl is on its way to zero, and round-off takes over.
+ring = (rho_w > r_ring[0]) & (rho_w < r_ring[1])
+fw.check_close("the curl of H against sigma E_z, away from the wire",
+ curl_z[ring], sigma_wire*E_wire(rho_w, t_c)[ring], rtol=0.05,
+ hint="Ampere with no source: d_x H_y - d_y H_x = sigma E_z")
+```
+
+:::{admonition} Solution — Task 4
+:class: dropdown
+
+```python
+curl_z = dHy_dx - dHx_dy
+```
+:::
+
+:::{admonition} The wire answers twice as fast as the sheet, and the curl dies with the electric field
+:class: important dropdown
+
+At 15 m in seawater the sheet peaks at 0.42 ms and the wire at 0.21 ms, a factor of two, because $t_p=\sigma\mu r^2/4p$ and $p$ went from $1/2$ to $1$. A steeper prefactor pulls the peak earlier: the field is being drained faster, so the moment at which arrival stops winning over decay comes sooner.
+
+The two magnetic pictures say different things. The map never changes shape, because the direction $\hat{\boldsymbol\varphi}$ is circular from the first instant and only the amount fills in. What takes time is the *amount*, and the second figure shows how far there is still to go. Expanding the exponential for $t\gg\sigma\mu\varrho^2/4$,
+
+$$ 1 - \exp\left(-\frac{\sigma\mu\varrho^2}{4t}\right) \approx \frac{\sigma\mu\varrho^2}{4t}, $$
+
+so the gap closes as $1/t$, which is the straight line of slope $-1$ on the logarithmic axes. A conductive ground is slower, and doubling the distance costs a factor of four in waiting.
+
+The curl map is the electric field map, rescaled. That is {eq}`eq:TEEy` read as a picture, and it is the sharpest statement of what makes a diffusive field different from a static one: the static magnetic field of the curl chapter is curl-free outside the wire, and this one is not, for exactly as long as there is an electric field to drive it. Let $t$ run on and both die, which is why the chapter says the magnetic field becomes a potential field at infinite time.
+:::
+
+---
+
+## Part 3 — Three dimensions, and one law for all of them
+
+A point source in three dimensions has the impulse response {eq}`eq:G3D`,
+
+$$ G(R,t) = \frac{(\sigma\mu)^{1/2}}{[4\pi t]^{3/2}}\exp\left(-\frac{\sigma\mu R^2}{4t}\right). $$
+
+Same kernel again. The prefactor is $t^{-3/2}$, so $p=3/2$.
+
+```{figure} figures/setup_ladder.svg
+:name: fig-lab6-ladder
+:width: 95%
+
+The three models of this lab side by side. The source spreads over a plane, a cylinder or a sphere, which is the only thing that changes between them.
+```
+
+### Task 5 — the ladder
+
+Write $G$. Your `t_peak` already covers it: the three geometries differ only in $p$, so at the same distance their peak times are in a fixed ratio. **Predict** that ratio before you run the last figure. Then one the formula does not hand you: which of the three records is the broadest on a logarithmic time axis, and which has the slowest tail? Say why before you look.
+
+```{code-cell} ipython3
+sigma_pt = 0.01 # sediments [S/m]
+
+# Task 5 -- one blank.
+def G_point(R, t, s=sigma_pt):
+ """eq:G3D -- the impulse response in three dimensions [1/(m s)]."""
+ return ___
+
+# --- given: a rotatable view of the field at one instant -------------------
+# The box is sized from the field, so changing t_show rescales both together:
+# the shell keeps its size on screen and the axis numbers change instead.
+t_show = 1.0e-3 # [s]
+D_show = np.sqrt(4*t_show / (sigma_pt*mu)) # how far it has spread [m]
+g_xy = np.linspace(-1.8*D_show, 1.8*D_show, 32) # the box follows the field
+g_z = np.linspace(0.0, 1.8*D_show, 32) # depth [m]
+GX, GY, GZ = np.meshgrid(g_xy, g_xy, g_z, indexing="ij")
+G_vol = G_point(np.sqrt(GX**2 + GY**2 + GZ**2), t_show)
+fw.show_isosurface(g_xy, g_xy, g_z, G_vol, levels=(0.08, 0.25, 0.6), opacity=0.32,
+ c_label="G", z_label="depth z [m]",
+ title=f"The point source at t = {t_show*1e3:g} ms. Drag to rotate")
+
+# --- given: the three geometries at the same distance, same ground ---------
+r0 = 300.0 # [m]
+t_l = np.logspace(np.log10(t_for_D(r0, sigma_pt)) - 2,
+ np.log10(t_for_D(r0, sigma_pt)) + 2, 4001) # [s]
+ladder = {"sheet, p = 1/2": np.abs(E_sheet(r0, t_l, sigma_pt)),
+ "wire, p = 1 ": np.abs(E_wire(r0, t_l, sigma_pt)),
+ "point, p = 3/2": G_point(r0, t_l)}
+fw.show_records(t_l, {k: v / np.max(v) for k, v in ladder.items()},
+ peaks={k: (t_peak(r0, p, sigma_pt), 1.0)
+ for k, p in zip(ladder, (0.5, 1.0, 1.5))},
+ t_label="t [s]", y_label="field / its own peak",
+ title=f"The same {r0:.0f} m in the same ground, three geometries")
+
+print(f"{'geometry':<12}{'p':>6}{'t_p [ms]':>12}{'ratio':>8}")
+for name, p in zip(("sheet", "wire", "point"), (0.5, 1.0, 1.5)):
+ tp = t_peak(r0, p, sigma_pt)
+ print(f"{name:<12}{p:>6}{tp*1e3:>12.3f}{tp/t_peak(r0, 1.5, sigma_pt):>8.1f}")
+
+fw.check_scalar("the point source peaks at sigma mu R^2 / 6",
+ t_peak(r0, 1.5, sigma_pt), t_l[np.argmax(G_point(r0, t_l))],
+ rtol=0.01, unit=" s")
+# integrate out to ten diffusion distances, where the Gaussian has nothing left
+_R = np.linspace(1e-3, 10*D_show, 60000) # a radial grid [m]
+fw.check_scalar("the integral of G over all space, which fixes the constant in front",
+ np.sum(G_point(_R, t_show) * 4*np.pi*_R**2) * (_R[1] - _R[0]),
+ 1.0 / (sigma_pt*mu), rtol=0.01, unit=" m^2/s")
+```
+
+:::{admonition} Solution — Task 5
+:class: dropdown
+
+```python
+def G_point(R, t, s=sigma_pt):
+ return np.sqrt(s*mu) / (4*np.pi*t)**1.5 * np.exp(-s*mu*R**2 / (4*t))
+```
+:::
+
+:::{admonition} The ladder, and why the sphere answers first
+:class: important dropdown
+
+$$ t_p = \frac{\sigma\mu r^2}{4p}:\qquad \frac{\sigma\mu r^2}{2},\quad \frac{\sigma\mu r^2}{4},\quad \frac{\sigma\mu r^2}{6} $$
+
+for $p=\tfrac12,1,\tfrac32$. The three peak times are in the ratio **6 : 3 : 2**, whatever the distance and whatever the ground, because the kernel is identical and only the prefactor differs. At 300 m in sediments that is 0.57 ms, 0.28 ms and 0.19 ms.
+
+The reason is worth saying in words. The exponential is the arrival: it rises from nothing as the diffusing field reaches $r$, and it does not care about the geometry. The prefactor is the decay: it is how fast the energy is spreading out, over a line, a cylinder or a sphere. The peak is where the two balance, so a geometry that spreads faster peaks earlier. The sphere spreads fastest, and answers first.
+
+The three-dimensional amplitude falls as $1/R^3$: double the distance and the peak is eight times smaller, against four for the wire and two for the sheet. Volume is expensive.
+:::
+
+---
+
+## Part 4 — A loop on the surface
+
+A real transmitter is a loop, not an infinite wire. A square loop of side $L_x=L_y=100$ m lies on the surface and its current is switched on at $t=0$. Integrating the three-dimensional Green's function along the four segments gives {eq}`eq:Exloop` and {eq}`eq:Eyloop`:
+
+$$ E_y(\boldsymbol r,t) = -\frac{\mu}{8\pi t}\left\{\exp\left(-\frac{x_m^2+z^2}{D^2}\right) - \exp\left(-\frac{x_p^2+z^2}{D^2}\right)\right\}\left[\mathrm{erf}\left(\frac{y_p}{D}\right) - \mathrm{erf}\left(\frac{y_m}{D}\right)\right] $$
+
+with $x_p=x+L_x/2$, $x_m=x-L_x/2$, and the **diffusion distance**
+
+$$ D = \sqrt{\frac{4t}{\sigma\mu}} , $$
+
+which is how far the field has spread by time $t$. The two exponentials are the two segments that run along $y$, at $x=\pm L_x/2$. They carry opposite currents, so they enter with opposite signs, and that difference is the only part you have to write.
+
+```{figure} figures/setup_loop.svg
+:name: fig-lab6-loop
+:width: 88%
+
+The model of this part, and the two planes you will cut through it. The current runs clockwise from vertex 1 to vertex 4, seen from above with $z$ pointing down.
+```
+
+### Task 6 — the loop forgets its shape
+
+Fill in the brace of $E_y$. Each of the three figures that follow then asks you one question, posed just above it.
+
+Integrating the point response along the four sides leaves a Gaussian across each pair of parallel wires and an error function along them. In the vertical plane through the middle of the loop you see the two segments that run along $y$, entering with opposite signs because they carry opposite currents.
+
+```{code-cell} ipython3
+sigma_loop, Lx, Ly = 0.3, 100.0, 100.0 # [S/m], loop sides [m]
+
+def D_of(t, s=sigma_loop):
+ """The diffusion distance [m]: how far the field has spread by time t."""
+ return np.sqrt(4*t / (s*mu))
+
+def Ex_loop(x, y, z, t, s=sigma_loop):
+ """eq:Exloop -- given, as the pattern to follow [V/m]."""
+ D = D_of(t, s)
+ xp, xm = x + Lx/2, x - Lx/2 # to the two segments that run along y
+ yp, ym = y + Ly/2, y - Ly/2 # to the two segments that run along x
+ return (-mu / (8*np.pi*t)
+ * (np.exp(-(yp**2+z**2)/D**2) - np.exp(-(ym**2+z**2)/D**2))
+ * (erf(xp/D) - erf(xm/D)))
+
+# Task 6 -- one blank. eq:Eyloop is eq:Exloop with x and y exchanged, except
+# that the two exponentials come in the opposite order, because the two
+# segments at x = +-Lx/2 carry opposite currents.
+def Ey_loop(x, y, z, t, s=sigma_loop):
+ """eq:Eyloop [V/m]."""
+ D = D_of(t, s)
+ xp, xm = x + Lx/2, x - Lx/2
+ yp, ym = y + Ly/2, y - Ly/2
+ return ___
+
+def E_loop(x, y, z, t, s=sigma_loop):
+ """Both horizontal components, as a pair."""
+ return Ex_loop(x, y, z, t, s), Ey_loop(x, y, z, t, s)
+
+# --- given: the vertical plane through the middle of the loop --------------
+n_L = 91
+x_L = np.linspace(-3*Lx, 3*Lx, n_L) # three loop widths either way [m]
+z_L = np.linspace(0.0, 3*Lx, n_L)
+XZ_x, XZ_z = np.meshgrid(x_L, z_L, indexing="ij")
+# from when the field has spread a tenth of the loop to when it has spread ten times it
+t_L = np.logspace(np.log10(t_for_D(Lx/10, sigma_loop)),
+ np.log10(t_for_D(Lx*10, sigma_loop)), 16) # [s]
+
+xz = np.array([E_loop(XZ_x, 0.0, XZ_z, tv)[1] for tv in t_L])
+fw.show_map(x_L, z_L, xz, t_L*1e3, value_name="t", value_fmt="{:.3g}", unit=" ms",
+ x_label="x [m]", y_label="depth z [m]", c_label="Ey [V/m]",
+ per_frame_scale=True, y_down=True, equal_aspect=False,
+ title="Two lobes of opposite sign, over the two wires, moving apart and down")
+```
+
+The same field, finely sampled, from the lecturer's own script. It runs over 200 instants instead of the 16 on your slider.
+
+```{video} figures/ELoopwsxz.mp4
+:width: 80%
+```
+
+The diffusion distance $D=\sqrt{4t/\sigma\mu}$ is the only length the field itself carries. Holding it against $L_x$ is what decides whether the pattern still remembers the shape of its source.
+
+**Question 1.** Drag the slider, which carries $D$. At which $D$ does the square stop looking square? Compare it with $L_x$.
+
+```{code-cell} ipython3
+# --- given: the plane at 10 m depth, seen from above. Arrows give direction ---
+z_xy = 10.0 # the depth of this slice [m]
+# its own narrower grid: the question is about the shape of the loop, and at
+# three loop widths either way the square is only a dozen pixels across
+x_xy = np.linspace(-1.5*Lx, 1.5*Lx, n_L) # [m]
+XY_x, XY_y = np.meshgrid(x_xy, x_xy, indexing="ij")
+EE = [E_loop(XY_x, XY_y, z_xy, tv) for tv in t_L]
+mag = np.array([np.hypot(ex, ey) for ex, ey in EE])
+# the slider carries D, not t, because the question is about D against Lx
+fw.show_map(x_xy, x_xy, mag, D_of(t_L),
+ arrows=(np.array([e[0] for e in EE]), np.array([e[1] for e in EE])),
+ arrow_every=10, log_c=True, per_frame_scale=True, y_down=True,
+ value_name="D", value_fmt="{:.0f}", unit=f" m (Lx = {Lx:g} m)",
+ x_label="x [m]", y_label="y [m]", c_label="log10(|E| / max)",
+ title=f"At z = {z_xy:g} m, seen from above")
+
+print(f"{'t [ms]':>9}{'D [m]':>9}{'D / Lx':>9}")
+for tv in t_L[::3]:
+ print(f"{tv*1e3:>9.3g}{D_of(tv):>9.1f}{D_of(tv)/Lx:>9.2f}")
+```
+
+Now time is held fixed and depth is stepped instead, with one colour scale shared by every frame so they can be compared directly.
+
+**Question 2.** Before you run it, write down with a neighbour what each of you expects to change as $z$ grows: the brightness, the width of the pattern, the direction of the arrows, or the position of the maximum. You will disagree about at least one of those. Run the cell, then settle it from {eq}`eq:Eyloop` rather than from the picture.
+
+```{code-cell} ipython3
+# --- given: the same instant, at a depth you choose -------------------------
+t_fix = 3e-4 # [s]
+# depths spaced against the diffusion distance, so the frames always span the drop
+# and stay inside the colour floor of 1e-4 that show_map clips at
+z_try = np.round(D_of(t_fix) * np.array([0.0, 0.6, 1.2, 1.8, 2.4])) # depth [m]
+deep = np.array([np.hypot(*E_loop(XY_x, XY_y, zv, t_fix)) for zv in z_try])
+fw.show_map(x_xy, x_xy, deep, z_try, value_name="z", value_fmt="{:.0f}", unit=" m",
+ log_c=True, per_frame_scale=False, y_down=True,
+ x_label="x [m]", y_label="y [m]", c_label="log10(|E| / max)",
+ title=f"The same instant, t = {t_fix*1e3:g} ms, at {len(z_try)} depths")
+
+print(f"at t = {t_fix*1e3:g} ms the diffusion distance is D = {D_of(t_fix):.0f} m\n")
+print(f"{'depth [m]':>11}{'|E| / its surface value':>26}")
+for zv in z_try:
+ f = np.exp(-zv**2 / D_of(t_fix)**2)
+ print(f"{zv:>11.0f}{f:>26.4f}")
+```
+
+Nested surfaces of constant $|\boldsymbol E|$ on a cube around the loop, so the interior stays visible. The box is sized from $D$ and keeps framing the field at whatever instant you set.
+
+**Question 3.** Rotate it, then set `t_3d = 3e-4` and run again. What has happened to the loop relative to the field it is driving?
+
+```{code-cell} ipython3
+# --- given: the loop in three dimensions. Drag to rotate --------------------
+t_3d = 3e-5 # [s]
+# wide enough for the loop, and for the field once it has spread past it
+half3 = max(1.8*Lx/2, 1.5*D_of(t_3d)) # [m]
+g3 = np.linspace(-half3, half3, 36)
+gz3 = np.linspace(0.0, half3, 36) # depth [m]
+L3x, L3y, L3z = np.meshgrid(g3, g3, gz3, indexing="ij")
+E3x, E3y = E_loop(L3x, L3y, L3z, t_3d)
+loop_wire = ([-Lx/2, Lx/2, Lx/2, -Lx/2, -Lx/2],
+ [-Ly/2, -Ly/2, Ly/2, Ly/2, -Ly/2], [0, 0, 0, 0, 0])
+fw.show_isosurface(g3, g3, gz3, np.hypot(E3x, E3y), levels=(0.06, 0.18, 0.45), opacity=0.3,
+ wire=loop_wire, wire_name="the loop", c_label="|E| [V/m]",
+ z_label="depth z [m]",
+ title=f"t = {t_3d*1e3:g} ms: D = {D_of(t_3d):.0f} m, "
+ f"D / Lx = {D_of(t_3d)/Lx:.2f}")
+
+fw.check("the field opposes the source current over the x = +Lx/2 wire",
+ E_loop(Lx/2, 0.0, z_xy, t_3d)[1] < 0 < E_loop(-Lx/2, 0.0, z_xy, t_3d)[1],
+ "the current there runs in +y, so an opposing field has E_y < 0")
+# A square loop is unchanged by exchanging x and y, and the circulation reverses
+# under that exchange, so E_y(x, y) = -E_x(y, x). This catches an erf bracket left
+# running over the wrong coordinate, which the sign check above does not.
+_q = np.array([-1.6, -0.4, 0.0, 0.7, 2.4]) * Lx
+_QX, _QY = np.meshgrid(_q, _q, indexing="ij")
+fw.check_close("E_y(x, y) against -E_x(y, x), which a square loop must satisfy",
+ Ey_loop(_QX, _QY, z_xy, t_3d), -Ex_loop(_QY, _QX, z_xy, t_3d), rtol=1e-9,
+ hint="the erf bracket of E_y must run over y, as its exponentials run over x")
+```
+
+:::{admonition} Solution — Task 6
+:class: dropdown
+
+```python
+def Ey_loop(x, y, z, t, s=sigma_loop):
+ D = D_of(t, s)
+ xp, xm = x + Lx/2, x - Lx/2
+ yp, ym = y + Ly/2, y - Ly/2
+ return (-mu / (8*np.pi*t)
+ * (np.exp(-(xm**2+z**2)/D**2) - np.exp(-(xp**2+z**2)/D**2))
+ * (erf(yp/D) - erf(ym/D)))
+```
+:::
+
+:::{admonition} What time does to the pattern, and what depth does
+:class: important dropdown
+
+**The square survives while $D$ is small.** At $t=0.01$ ms, $D=10$ m against a 100 m loop and the four wires are separate bright lines. By $t=0.1$ ms, $D=33$ m, and the pattern is already round. The loop stops looking square at about $D\approx L_x/3$, which is $t\approx\sigma\mu L_x^2/36$, here 0.1 ms. After that no measurement can tell this loop from a round one of the same area. A source forgets its shape once the field has diffused further than the source is big, and this is the same blurring that merged two Gaussians in the last lab.
+
+**Depth only dims it.** Put $z$ into {eq}`eq:Eyloop` and the depth appears as $\exp[-(x_m^2+z^2)/D^2]$, which factors into $\exp(-z^2/D^2)$ times the surface expression, identically for both components. So
+
+$$ \boldsymbol E(x,y,z,t) = \exp\left(-\frac{z^2}{D^2}\right)\boldsymbol E(x,y,0,t). $$
+
+Going deeper multiplies the whole picture by one number. It does not blur it, does not rotate it, does not change which way the arrows point. The frames of the depth figure are the same image at different brightnesses, and the printed column is that one factor.
+
+That factor is the whole of depth sounding. At $t=0.3$ ms, $D=56$ m, so 50 m down you still have 46% of the surface field and 200 m down you have $3.5\times10^{-6}$ of it. You can only see about as deep as $D$, and $D$ grows as $\sqrt t$. To look twice as deep you must wait four times as long, which is Part 1's $t_p\propto z^2$ arriving from the other direction.
+
+**The arrows oppose the current.** Over the wire at $x=+L_x/2$ the current runs in $+\hat{\boldsymbol y}$ and $E_y$ is negative. The switched-on current drives a field that opposes it, which is Lenz's law seen in a conductor.
+:::
+
+---
+
+## Part 5 — Designing a sounding
+
+You now have everything a time-domain electromagnetic survey rests on. This part asks you to put it together. Answer the four questions first, from the formulas and not from a new calculation, then run the cell.
+
+1. A target lies 400 m down in ground of conductivity $10^{-2}$ S/m. Roughly when does the response from that depth arrive? Which $p$ did you put into `t_peak`, and does the choice matter at the factor-of-two level?
+2. Your recording window is 20 µs to 10 ms. Which of the three grounds of Part 1 puts a 400 m target inside that window?
+3. You have two loops, 50 m and 400 m on a side. The target is 400 m down. Which do you choose, and what decides it: the loop size, or the time you record for?
+4. The field at the target is $\exp(-z^2/D^2)$ of its surface value. Below which amplitude would you say the target is not being illuminated at all, and what does that make the deepest usable depth?
+
+Question 3 is the one worth arguing about, because the intuitive answer is wrong and nothing on the page has told you otherwise yet. Commit to a loop before you run the second cell, which measures both.
+
+Everything here is $\exp(-z^2/D^2)$ at the target depth, with $D$ built from each ground's conductivity, against time. A target is worth looking for where that curve is still high inside the window you can actually record in.
+
+```{code-cell} ipython3
+# --- given: the window, the ground and the depth, all on one picture -------
+z_target = 400.0 # [m]
+t_win = (2e-5, 1e-2) # the recording window [s]
+_tp_t = [t_peak(z_target, 0.5, sv) for sv in media.values()]
+t_d = np.logspace(min(np.log10(min(_tp_t)), np.log10(t_win[0])) - 2,
+ max(np.log10(max(_tp_t)), np.log10(t_win[1])) + 2, 2001)
+
+fig, ax = plt.subplots(figsize=(8.2, 3.6))
+for k, (name, sv) in enumerate(media.items()):
+ ax.plot(t_d, np.exp(-z_target**2 * sv * mu / (4*t_d)), color=f"C{k}", lw=2, label=name)
+ ax.axvline(t_peak(z_target, 0.5, sv), color=f"C{k}", ls=":", lw=1.2)
+ax.axvspan(*t_win, color="0.85", zorder=0, label="recording window")
+ax.axhline(1/np.e, color="k", lw=0.8, ls="--")
+ax.text(0.01, 1/np.e + 0.02, "z = D", fontsize=8, transform=ax.get_yaxis_transform())
+ax.set_xscale("log"); ax.set_xlabel("t [s]")
+ax.set_ylabel(f"illumination at {z_target:.0f} m, exp(-z²/D²)")
+ax.set_title(f"Which ground puts a {z_target:.0f} m target inside the window?")
+ax.grid(alpha=0.3); ax.legend(fontsize=8, loc="upper left")
+plt.show()
+
+print(f"a target at {z_target:.0f} m, window {t_win[0]*1e6:.0f} us to {t_win[1]*1e3:.0f} ms\n")
+print(f"{'ground':<30}{'t_p [s]':>11}{'in window?':>12}"
+ f"{f'illum. at {t_win[1]*1e3:g} ms':>18}")
+for name, sv in media.items():
+ tp = t_peak(z_target, 0.5, sv)
+ lit = np.exp(-z_target**2 * sv * mu / (4*t_win[1]))
+ print(f"{name:<30}{tp:>11.3g}{'yes' if t_win[0] < tp < t_win[1] else 'no':>12}{lit:>18.3g}")
+```
+
+Now the loop itself, the one you built in Part 4, at two sizes. Both are placed over the same target in the same ground, so the only thing different is the side length.
+
+```{code-cell} ipython3
+# --- given: a 50 m loop and a 400 m loop over the same target --------------
+_t_scan = np.logspace(-5, -1, 400) # [s]
+_Lx0, _Ly0 = Lx, Ly # keep your Part 4 loop
+for _L in (50.0, 400.0):
+ Lx = Ly = _L # E_loop reads these
+ _E = np.array([np.hypot(*E_loop(100.0, 0.0, z_target, tv, s=0.01)) for tv in _t_scan])
+ print(f"L = {_L:3.0f} m: peak |E| = {_E.max():.2e} V/m at t = {_t_scan[_E.argmax()]*1e3:5.2f} ms,"
+ f" D there = {np.sqrt(4*_t_scan[_E.argmax()]/(0.01*mu)):.0f} m")
+Lx, Ly = _Lx0, _Ly0
+
+print(f"\nthe loop's own arrival, t_p with p = 5/2 at R = sqrt(100^2 + {z_target:.0f}^2):"
+ f" {t_peak(np.hypot(100.0, z_target), 2.5, 0.01)*1e3:.2f} ms")
+```
+
+:::{admonition} Answers
+:class: dropdown
+
+**1.** $t_p=\sigma\mu z^2/4p$. The sheet's $p=1/2$ gives 1.0 ms, the wire 0.50 ms, the point source 0.34 ms, and the loop's $p=5/2$ gives 0.20 ms. A survey uses a loop, so 0.2 ms is the honest number, but any of them sets the scale: the spread is a factor of five across geometries that differ completely. What matters is $t\sim\sigma\mu z^2$, not the number in the denominator.
+
+**2.** Sediments, at 1.0 ms, sit comfortably inside, and the printed table says so. Seawater peaks at 0.40 s, forty times beyond the end of the window: the window closes long before the response from 400 m arrives. Basement peaks at 10 µs, which is before the window opens, so its response from 400 m is already over by the time the first sample is taken. Only the middle ground is workable, and that is typical.
+
+**3.** The time you record for. Depth of investigation is set by $D=\sqrt{4t/\sigma\mu}$, which contains no loop size at all, and the cell measures what an eightfold change in $L$ actually does to the arrival: the peak moves from 0.215 ms to 0.271 ms and the diffusion distance at the peak from 262 m to 294 m. That is 26% and 12%, against a factor of 8 in the loop. The 50 m loop lands on $t_p$ with $p=5/2$ at $R=\sqrt{100^2+400^2}$, 0.214 ms, because at that size it is already a point as far as the target is concerned; the 400 m loop is the later of the two because the target now sits under the middle of it rather than off to one side. What the big loop buys is amplitude, a factor of 34 at the target, and that decides whether you can detect the response rather than where the response comes from. A bigger loop is worth having for exactly that reason, and for no other. You buy depth with late time and a quiet receiver.
+
+**4.** Any threshold in the range 1% to 10% gives nearly the same rule, because the Gaussian is so steep. Taking $1/e$ gives $z_{\max}=D=\sqrt{4t/\sigma\mu}$ exactly, which is the standard definition of the diffusion depth. Read the two ends of the window off the printed column. At the 1.0 ms peak in sediments $D=564$ m, so the 400 m target sits at $\exp(-0.50)=0.60$ of the surface field; by the 10 ms end of the window $D$ has grown to 1784 m and the target is at 0.95, which is what the table prints. In seawater 10 ms gives only $D=89$ m and the target is at $2\times10^{-9}$, which is nothing at all. Illumination is not what limits you in sediments; signal strength is.
+:::
+
+:::{admonition} What this lab was about
+:class: important dropdown
+
+One kernel, $\exp(-\sigma\mu r^2/4t)$, appeared in every part. It came from the heat equation last week with $\kappa$ replaced by $1/(\sigma\mu)$, and the only thing that changed from a sheet to a wire to a point was the power of $t$ in front, which is a statement about geometry and not about electromagnetism.
+
+From that one kernel came three results a geophysicist uses directly. Arrival time measures distance, $t_p=\sigma\mu r^2/4p$. Depth of investigation is the diffusion distance, $z\approx D=\sqrt{4t/\sigma\mu}$, and it is bought with time rather than with equipment. And a source is blurred into a point once $D$ exceeds its size, which is why the shape of a transmitter stops mattering, and also why deep structure is always seen blurred.
+
+The magnetic field made the fourth point. For as long as the electric field is there, $(\nabla\times\boldsymbol H)_z=\sigma E_z$ is not zero, and the magnetic field is not a potential field. It becomes one only in the limit, which is where the magnetostatics of the earlier chapters lives.
+:::
diff --git a/book/4_mechanical_waves/acoustic_waves.md b/book/4_mechanical_waves/acoustic_waves.md
index 48a52a0..16ca9fa 100644
--- a/book/4_mechanical_waves/acoustic_waves.md
+++ b/book/4_mechanical_waves/acoustic_waves.md
@@ -48,7 +48,7 @@ $$
in which $c$ can be seen as the velocity of the wave, for which we have: $c = \sqrt{K/ \rho}$.
-This is the basic wave equation for one-dimensional acoustic waves. Equally well, we can derive from the same two equations {eq}`eq:eqdeform` and {eq}`eq:eqmotion` (see also *EXERCISES*) that the particle velocity also satisfies the wave equation:
+This is the basic wave equation for one-dimensional acoustic waves. Equally well, we can derive from the same two equations {eq}`eq:eqdeform` and {eq}`eq:eqmotion` (see also the [*EXERCISES*](./exercises.md)) that the particle velocity also satisfies the wave equation:
$$
\frac{\partial^2 v_x}{\partial x^2}
@@ -57,49 +57,8 @@ $$ (eq:waveqacv)
where $c$ is the same as in the wave equation of the pressure.
-The solution to the wave equation for the pressure is:
+Actually, the movie we showed before was generated using the above equations (but then discretized), so the equation of deformation {eq}`eq:eqdeform` and the equation of motion {eq}`eq:eqmotion`. It is shown again below.
-$$
-p(x,t) = s(t \pm x/c)
-$$ (eq:pxt)
-
-where $s(t)$ is some function. Note the dependency on space and time via $(t \pm x/c)$, which denotes that it is a travelling wave, and behaves as a direct wave as discussed in [](./snells_law.md). The sign in the argument depends on the direction the wave is travelling.
-
-Often, seismic responses are analyzed in terms of frequencies, i.e., Fourier spectra. The definition used here is:
-
-$$
-G(\omega) = \int_{-\infty}^{+\infty} g(t) \exp(- i \omega t) \, dt
-$$ (eq:fourier1)
-
-where the radial frequency $\omega$ is used, so $\omega = 2 \pi f$ where $f$ is the (ordinary) frequency. Using this convention, it is easy to show that differentiation with respect to time is equivalent to multiplication with $i\omega$ in the Fourier domain. When we transform the solution of the wave equation ({eq}`eq:pxt`) to the Fourier domain, we obtain:
-
-$$
-P(x,\omega) = S(\omega) \exp (\pm i\omega x/c).
-$$ (eq:Pxw)
-
-Note here that the time delay $x/c$ (in the time domain) becomes a *linear* phase in the Fourier domain, so the phase $\pm \omega x/c$ is linear as a function of $\omega$.
-
-In the above, we gave an expression for the pressure, but one can also derive an equivalent expression for the particle velocity $v_x$. For that purpose, we can use the equation of motion as expressed in {eq}`eq:eqmotion`, but then in its Fourier-transformed version, which is:
-
-$$
-V_x(x,\omega) = -\frac{1}{i\omega \rho} \frac{\partial P(x,\omega)}{\partial x}.
-$$
-
-When we substitute the solution for the pressure from above ({eq}`eq:Pxw`), we get for the negative sign:
-
-$$
-\begin{aligned}
-V_x(x,\omega) & = -\frac{1}{i \omega \rho}
- S(\omega) \frac{- i \omega}{c} \exp(-i\omega x/c)
-\\
- & = S(\omega) \frac{1}{\rho c} \exp(-i\omega x/c).
-\end{aligned}
-$$
-
-Notice that the particle velocity is a scaled version of the pressure:
-
-$$
-V_x(x,\omega) = \frac{P(x,\omega)}{\rho c}.
-$$
-
-The scaling factor is ($\rho c$), being called the *seismic impedance*.
+```{video} figures/seismicwave_1layer_vz.mp4
+:width: 80%
+```
diff --git a/book/4_mechanical_waves/data/shot_group_0.npz b/book/4_mechanical_waves/data/shot_group_0.npz
new file mode 100644
index 0000000..9cb9a49
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_0.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_1.npz b/book/4_mechanical_waves/data/shot_group_1.npz
new file mode 100644
index 0000000..710b7da
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_1.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_10.npz b/book/4_mechanical_waves/data/shot_group_10.npz
new file mode 100644
index 0000000..43711c9
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_10.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_11.npz b/book/4_mechanical_waves/data/shot_group_11.npz
new file mode 100644
index 0000000..0ef2489
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_11.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_12.npz b/book/4_mechanical_waves/data/shot_group_12.npz
new file mode 100644
index 0000000..086c0fb
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_12.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_13.npz b/book/4_mechanical_waves/data/shot_group_13.npz
new file mode 100644
index 0000000..ce23250
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_13.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_14.npz b/book/4_mechanical_waves/data/shot_group_14.npz
new file mode 100644
index 0000000..69248b9
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_14.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_15.npz b/book/4_mechanical_waves/data/shot_group_15.npz
new file mode 100644
index 0000000..c4331de
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_15.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_2.npz b/book/4_mechanical_waves/data/shot_group_2.npz
new file mode 100644
index 0000000..29c66a3
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_2.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_3.npz b/book/4_mechanical_waves/data/shot_group_3.npz
new file mode 100644
index 0000000..27d0fb7
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_3.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_4.npz b/book/4_mechanical_waves/data/shot_group_4.npz
new file mode 100644
index 0000000..0fe9ead
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_4.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_5.npz b/book/4_mechanical_waves/data/shot_group_5.npz
new file mode 100644
index 0000000..89549e9
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_5.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_6.npz b/book/4_mechanical_waves/data/shot_group_6.npz
new file mode 100644
index 0000000..2eedbfe
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_6.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_7.npz b/book/4_mechanical_waves/data/shot_group_7.npz
new file mode 100644
index 0000000..e4c046a
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_7.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_8.npz b/book/4_mechanical_waves/data/shot_group_8.npz
new file mode 100644
index 0000000..f41fb7a
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_8.npz differ
diff --git a/book/4_mechanical_waves/data/shot_group_9.npz b/book/4_mechanical_waves/data/shot_group_9.npz
new file mode 100644
index 0000000..2e2dab2
Binary files /dev/null and b/book/4_mechanical_waves/data/shot_group_9.npz differ
diff --git a/book/4_mechanical_waves/elastic_waves.md b/book/4_mechanical_waves/elastic_waves.md
index d6a60de..51b03e7 100644
--- a/book/4_mechanical_waves/elastic_waves.md
+++ b/book/4_mechanical_waves/elastic_waves.md
@@ -105,4 +105,4 @@ With respect to the (1-D) reflection coefficients for P- and S-waves at the boun
- For P-waves, we require the continuity of the tensile stress $\tau_{xx}$ and the particle velocity $v_x$ at the interface, and
- For S-waves we require the continuity of the shear stress $\tau_{zx}$ and the particle velocity $v_z$ at the interface.
-The resulting expressions for the reflection and transmission coefficient are the same as for the acoustic case, only now the velocities are for the elastic solid. (See *EXERCISES*)
+The resulting expressions for the reflection and transmission coefficient are the same as for the acoustic case, only now the velocities are for the elastic solid. (See the [*EXERCISES*](./exercises.md))
diff --git a/book/4_mechanical_waves/exercise_4a.md b/book/4_mechanical_waves/exercise_4a.md
new file mode 100644
index 0000000..d41c07a
--- /dev/null
+++ b/book/4_mechanical_waves/exercise_4a.md
@@ -0,0 +1,185 @@
+---
+jupytext:
+ text_representation:
+ extension: .md
+ format_name: myst
+ format_version: 0.13
+kernelspec:
+ display_name: Python 3 (ipykernel)
+ language: python
+ name: python3
+mystnb:
+ # The cells are to be filled in and run by the reader, so the page must not be
+ # executed at build time.
+ execution_mode: 'off'
+---
+
+# Exercise 4a: Model traveltimes for seismic wavefield from one shot
+
+In this exercise you will model the traveltimes of a seismic wavefield for a direct, refracted and reflected wave. The model contains only one boundary, giving rise to a reflection and a head wave (a critically refracted wave) when the wave speed of the second medium is higher than the wave speed of the first medium.
+
+To this aim, we will use a synthetically generated wave field for one shot position. Sixteen (0-15) wave-field recordings are present:
+
+- Take the first character of your first name, take the number of this character in the alphabet and if it is larger than 15, subtract 15.
+- Make this your group number as a string variable.
+
+```{code-cell} ipython3
+your_group_no = ""
+```
+
+```{code-cell} ipython3
+# Housekeeping, no physics. Running this page in the browser, the next lines fetch
+# your own shot gather from the book. Running the notebook on your own machine, with
+# the .npz file already beside it, they do nothing.
+import os
+import pathlib
+
+filename_npz = "shot_group_" + your_group_no + ".npz"
+
+if not pathlib.Path(filename_npz).exists():
+ from pyodide.http import pyfetch
+ # Live Code runs Python in a worker, which would resolve a plain "data/..."
+ # against itself and not against this page. The working directory mirrors the
+ # folder this page sits in on the server, so it gives the address to ask for.
+ url = os.getcwd().split("/home/pyodide/book")[-1] + "/data/" + filename_npz
+ response = await pyfetch(url)
+ if response.status != 200:
+ raise FileNotFoundError(
+ url + " -> HTTP " + str(response.status)
+ + ". Is your_group_no a number from 0 to 15?")
+ pathlib.Path(filename_npz).write_bytes(await response.bytes())
+ print("fetched", filename_npz)
+
+# The hover read-out on the figures needs ipympl, which the browser lacks.
+try:
+ import ipympl # noqa: F401
+except ModuleNotFoundError:
+ import micropip # browser only; locally, pip install ipympl
+ await micropip.install("ipympl")
+
+# Cosmetic only. ipympl stacks its toolbar vertically, which this theme stretches
+# across the page, and it names two of its icons in Font Awesome 4 while the book
+# ships Font Awesome 6.
+try:
+ from ipympl.backend_nbagg import Canvas
+ from IPython.display import display, HTML
+ Canvas.toolbar_position.default_value = "top"
+ display(HTML(r""))
+except Exception:
+ pass # the figures work regardless of how it looks
+```
+
+Now you will read this wavefield data set, that includes the geometry of the data, so the positions of the source and receivers. To that end:
+
+- Run the next code, and consider/analyze the source and receiver configuration.
+- Where is the source positioned compared to the receiver positions (beginning, end, or somewhere else)?
+
+```{code-cell} ipython3
+#%matplotlib qt
+%matplotlib ipympl
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Reading NPZ file
+filename_npz = "shot_group_" + your_group_no + ".npz"
+
+with open(filename_npz, "rb") as f:
+ npzfile = np.load(f) ; print(npzfile.files)
+ Data = npzfile['Data']
+ timevector = npzfile['tt']
+ x_receivers = npzfile['x_receivers']
+ x_source = npzfile['x_source']
+
+# Check:
+#print(Data[1,1])
+#print(timevector[1])
+#print(x_receivers[1])
+#print(x_source)
+
+plt.close("geometry") # so a re-run redraws instead of stacking
+plt.figure( "geometry", figsize = ( 12, 6 ), dpi = 60 ) # fits the canvas to the page
+plt.plot(x_receivers,'o',color='b')
+plt.plot(x_source,'x',color='r', markersize=10)
+plt.legend(['x_receivers','x_source'])
+plt.xlabel('index'); plt.ylabel('Distance (m)')
+plt.title(' Distance of source and receivers versus index of array', fontsize=14 )
+plt.show()
+```
+
+Next, you need to generate the travel-time curves for that geometry. But, of course, you need to specify the parameters of the model, i.e., the velocities of the different waves and the depth of the boundary. The depth of the boundary may be different for the refractor and the reflector, but initially you can keep those depths the same. So to that end, using the variable names as set in the next code:
+
+- Specify the P-wave velocities of the first and second medium, and the surface-wave velocity (of the first medium)
+- Specify the depths of the refractor and reflector (initially, take them the same)
+- Play with the velocities and depths such that you can see the different events in the travel-time curves, in particular that you can see both a reflection and a head wave (critically refracted wave)
+
+```{code-cell} ipython3
+# Specify velocities and depths
+c_directP =
+c_directSurf =
+c_refract =
+depth_refract =
+depth_reflect =
+```
+
+```{code-cell} ipython3
+tt_directP = np.zeros_like(x_receivers)
+tt_directSurf = np.zeros_like(x_receivers)
+tt_refract = np.zeros_like(x_receivers)
+tt_reflect = np.zeros_like(x_receivers)
+
+# direct P-wave:
+tt_directP[:] = abs(x_receivers[:]-x_source)/c_directP ; print('tt_directP = ',tt_directP[0:3])
+
+# direct Surface wave:
+tt_directSurf[:] = abs(x_receivers[:]-x_source)/c_directSurf ; print('tt_directSurf = ',tt_directSurf[0:3])
+
+# refracted P-wave:
+sin_crit = c_directP/c_refract
+cos_crit = np.sqrt(1-sin_crit*sin_crit)
+tt_refract[:] = abs(x_receivers[:]-x_source)/c_refract + 2.0*depth_refract*cos_crit/c_directP ; print('tt_refract = ',tt_refract[0:3])
+
+# find limiting indices of array when distance < x_crit (then no critical refraction):
+x_crit = 2.0 * depth_refract * sin_crit/cos_crit ; print('x_crit = ',x_crit,' m')
+# left end:
+i1 = 0
+ii_end = len(x_receivers) #; print('iiend = ',ii_end)
+for ii in range (0,len(x_receivers)):
+ if abs(x_receivers[ii]-x_source) < x_crit:
+ break
+ else:
+ i1 = ii
+#print('i1 =',i1,'; x_receivers[i1] = ',x_receivers[i1])
+#print('abs(x_receivers[i1]-x_source) = ',abs(x_receivers[i1]-x_source))
+
+# right end:
+i2 = len(x_receivers)
+for ii in range (ii_end-1,0,-1):
+ if abs(x_receivers[ii]-x_source) < x_crit:
+ break
+ else:
+ i2 = ii
+#print('i2 =',i2,'; x_receivers[i2] = ',x_receivers[i2])
+#print('abs(x_receivers[i2]-x_source) = ',abs(x_receivers[i2]-x_source))
+
+
+# reflected P-wave:
+t0 = 2*depth_reflect/c_directP
+tt_reflect[:] = np.sqrt(t0*t0 + (x_receivers[:]-x_source)*(x_receivers[:]-x_source)/(c_directP*c_directP))
+print('tt_reflect = ',tt_reflect[0:3])
+
+plt.close("traveltimes") # so a re-run redraws instead of stacking
+plt.figure( "traveltimes", figsize = ( 16, 12 ), dpi = 45 ) # fits the canvas to the page
+plt.plot(x_receivers, tt_directP, 'x',color='b')
+plt.plot(x_receivers, tt_directSurf, 'x',color='c')
+plt.plot(x_receivers[0:i1], tt_refract[0:i1], 'x',color='g')
+plt.plot(x_receivers[i2:-1], tt_refract[i2:-1], 'x',color='g')
+plt.plot(x_receivers, tt_reflect, 'x',color='r')
+plt.axis( [np.min(x_receivers), np.max(x_receivers), np.min(timevector), np.max(timevector)] )
+plt.gca().invert_yaxis()
+plt.legend(['directP','directSurf','refraction','','reflection'])
+plt.show()
+```
diff --git a/book/4_mechanical_waves/exercise_4b.md b/book/4_mechanical_waves/exercise_4b.md
new file mode 100644
index 0000000..fbfc919
--- /dev/null
+++ b/book/4_mechanical_waves/exercise_4b.md
@@ -0,0 +1,193 @@
+---
+jupytext:
+ text_representation:
+ extension: .md
+ format_name: myst
+ format_version: 0.13
+kernelspec:
+ display_name: Python 3 (ipykernel)
+ language: python
+ name: python3
+mystnb:
+ # The cells are to be filled in and run by the reader, so the page must not be
+ # executed at build time.
+ execution_mode: 'off'
+---
+
+# Exercise 4b: Interpret velocities and depths from seismic wavefield from one shot
+
+In this exercise you will interpret direct, refracted and reflected waves from a (synthetically generated) seismic wavefield recording. And will check the results via forward modelling, as you did in [Exercise 4a](./exercise_4a.md).
+
+To this aim, use the same synthetically generated wave field that you used before:
+
+- Take the first character of your first name, take the number of this character in the alphabet and if it is larger than 15, subtract 15.
+- Make this your group number as a string variable.
+
+```{code-cell} ipython3
+your_group_no = ""
+```
+
+```{code-cell} ipython3
+# Housekeeping, no physics. Running this page in the browser, the next lines fetch
+# your own shot gather from the book. Running the notebook on your own machine, with
+# the .npz file already beside it, they do nothing.
+import os
+import pathlib
+
+filename_npz = "shot_group_" + your_group_no + ".npz"
+
+if not pathlib.Path(filename_npz).exists():
+ from pyodide.http import pyfetch
+ # Live Code runs Python in a worker, which would resolve a plain "data/..."
+ # against itself and not against this page. The working directory mirrors the
+ # folder this page sits in on the server, so it gives the address to ask for.
+ url = os.getcwd().split("/home/pyodide/book")[-1] + "/data/" + filename_npz
+ response = await pyfetch(url)
+ if response.status != 200:
+ raise FileNotFoundError(
+ url + " -> HTTP " + str(response.status)
+ + ". Is your_group_no a number from 0 to 15?")
+ pathlib.Path(filename_npz).write_bytes(await response.bytes())
+ print("fetched", filename_npz)
+
+# The hover read-out on the figures needs ipympl, which the browser lacks.
+try:
+ import ipympl # noqa: F401
+except ModuleNotFoundError:
+ import micropip # browser only; locally, pip install ipympl
+ await micropip.install("ipympl")
+
+# Cosmetic only. ipympl stacks its toolbar vertically, which this theme stretches
+# across the page, and it names two of its icons in Font Awesome 4 while the book
+# ships Font Awesome 6.
+try:
+ from ipympl.backend_nbagg import Canvas
+ from IPython.display import display, HTML
+ Canvas.toolbar_position.default_value = "top"
+ display(HTML(r""))
+except Exception:
+ pass # the figures work regardless of how it looks
+```
+
+Now you will read this wavefield data set, that includes the geometry of the data and data themselves. To that end:
+
+- Run the next code, consider/analyze the source and receiver configuration again, and consider/analyze the wavefield data.
+- In the code, the maximum amplitude (`maxAmplitude`) is set; adapt this maximum amplitude and look at the effect on the plot.
+
+```{code-cell} ipython3
+#%matplotlib qt
+%matplotlib ipympl
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Reading NPZ file
+filename = "shot_group_" + your_group_no
+filename_npz = filename + ".npz"
+
+with open(filename_npz, "rb") as f:
+ npzfile = np.load(f) ; print(npzfile.files)
+ Data = npzfile['Data']
+ timevector = npzfile['tt']
+ x_receivers = npzfile['x_receivers']
+ x_source = npzfile['x_source']
+
+# Check:
+#print(Data[1,1])
+#print(timevector[1])
+#print(x_receivers[1])
+#print(x_source)
+
+# Plot geometry of shot gather
+plt.close("geometry") # so a re-run redraws instead of stacking
+plt.figure( "geometry", figsize = ( 12, 6 ), dpi = 60 ) # fits the canvas to the page
+plt.plot(x_receivers,'o',color='b')
+plt.plot(x_source,'x',color='r', markersize=10)
+plt.legend(['x_receivers','x_source'])
+plt.xlabel('index'); plt.ylabel('Distance (m)')
+plt.title(' Distance of source and receivers versus index of array', fontsize=14 )
+plt.show()
+```
+
+```{code-cell} ipython3
+# Plot extracted shot gather
+plt.close("shot gather") # so a re-run redraws instead of stacking
+plt.figure( "shot gather", figsize = ( 16, 12 ), dpi = 45 ) # fits the canvas to the page
+maxAmplitude = 0.04 * np.max( np.abs(Data) )
+extent = [ np.min(x_receivers), np.max(x_receivers), np.max(timevector), np.min(timevector) ]
+# Standard colormap:
+# plt.imshow( Data, vmin=-maxAmplitude, vmax = maxAmplitude, extent = extent, aspect = 'auto' )
+plt.imshow( Data, vmin=-maxAmplitude, vmax = maxAmplitude, extent = extent, aspect = 'auto', cmap='seismic' )
+titleplot = filename + " : original"
+plt.title(titleplot, fontsize=14 )
+plt.colorbar()
+plt.show()
+```
+
+## Identification of events
+
+The figure of the wavefield data is a simulation of a seismic wavefield recording. On this wavefield recording, determine which event(s) can be interpreted as:
+
+- Direct P-wave
+- (Critically) refracted P-wave
+- Direct surface wave
+- P-wave reflections
+
+You may have to play with `maxAmplitude` to see some events better.
+
+## Estimation of wave velocity of events
+
+Based on this identification, determine the velocity of the:
+
+- Direct P-wave
+- Head wave (Critically refracted P-wave)
+- Direct surface wave
+
+In order to estimate the velocity of a reflected wave, first consider the equation for a reflection-hyperbola, i.e.,
+
+$$
+T^2 = T_0^2 + \frac{X^2}{c^2},
+$$
+
+where $T_0$ is the time at the apex of the hyperbola, and $X$ is the horizontal distance from the apex, so the relative horizontal distance between the source and receiver. From this equation, the velocity can be made explicit:
+
+$$
+c = \frac{X}{\left( T^2 - T_0^2 \right)^{1/2}}.
+$$
+
+Based on this:
+
+- Estimate the velocity from one of the reflections you picked above.
+
+## Estimation of depth of refractor and reflector
+
+Based on the above, depth estimates can be made. Let us first estimate the depth of the **refractor**, using the cross-over distance, the distance where the direct wave and the head-wave (critically refracted) cross each other. This is given by:
+
+$$
+Z_{\mathrm{refractor}} = \frac{X_{\mathrm{cross}}}{2}
+ \left[ \frac{c_2-c_1}{c_2+c_1} \right]^{1/2},
+$$
+
+where $c_1$ is the velocity of the direct P-wave and $c_2$ the velocity of the (P-wave) refractor.
+
+Based on this:
+
+- Determine the depth of the refractor
+
+Let us now consider the **reflection** you took. At a distance $X=0$, the reflection time becomes $T_0 = 2z/c_{\mathrm{refl}}$. That reflection time can be picked from the wavefield recording, and the velocity you estimated earlier. Then the depth can be determined via $z = c_{\mathrm{refl}} T_0 /2$.
+
+Based on this:
+
+- Determine the depth of the reflector
+
+## Final check: forward model results via Exercise 4a
+
+Based on the above, you have now determined all the velocities and depths for which you can model the travel times as you did in [Exercise 4a](./exercise_4a.md). So therefore:
+
+- Run the travel-time code of Exercise 4a again, but now with the model parameters you estimated for your wavefield recording
+- Check whether the travel times match the ones from your seismic wavefield recording above.
+
+This should give you confidence in your physical interpretation of the wavefield recording.
diff --git a/book/4_mechanical_waves/exercises.md b/book/4_mechanical_waves/exercises.md
new file mode 100644
index 0000000..e09167b
--- /dev/null
+++ b/book/4_mechanical_waves/exercises.md
@@ -0,0 +1,29 @@
+# Exercises on fundamental physical laws, P- and S-waves
+
+1. From the deformation equation and the equation of motion for acoustic waves ({eq}`eq:eqdeform` and {eq}`eq:eqmotion`) show that the particle velocity $v_x$ also satisfies the wave equation, i.e., {eq}`eq:waveqacv`.
+
+2. Here we look at deriving the wave equation in a different way:
+
+ - Show that from substituting the deformation equation {eq}`eq:eqdeform` and the equation of motion {eq}`eq:eqmotion` in the commutative property for the pressure, i.e.,
+
+ $$
+ \frac{\partial}{\partial x} \left( \frac{\partial p}{\partial t} \right)
+ =
+ \frac{\partial}{\partial t} \left( \frac{\partial p}{\partial x} \right),
+ $$
+
+ the wave equation for the particle velocity is nicely obtained.
+ - In the same way show that from substituting the deformation equation {eq}`eq:eqdeform` and the equation of motion {eq}`eq:eqmotion` in the commutative property for the *particle velocity*, as expressed in {eq}`eq:commut`, the wave equation for the pressure is nicely obtained.
+
+3. From the deformation equation and the equation of motion for shear waves ({eq}`eq:deformS` and {eq}`eq:eqmotionS`) show that the particle velocity $v_z$ also satisfies the wave equation.
+
+4. Here we consider shear waves:
+
+ - How can the solution of the wave equation for the shear stress ({eq}`eq:waveqStau`) be written?
+ - Based on this solution, *derive* the solution for the particle velocity $v_z$.
+ - What is the difference between the two? And how is that called?
+
+5. Here we consider shear waves arriving at an interface:
+
+ - When we look at a reflection and transmission across an interface between two materials: we have an incident S-wave coming from above. How would you write the wave field (shear stress and particle velocity) above the interface? And how below?
+ - At the interface, we require the continuity of the shear stress and particle velocity, as given above: from this, *derive* the reflection and transmission coefficients.
diff --git a/book/4_mechanical_waves/figures/seismicwave_1layer_vz.mp4 b/book/4_mechanical_waves/figures/seismicwave_1layer_vz.mp4
new file mode 100644
index 0000000..18a2de9
Binary files /dev/null and b/book/4_mechanical_waves/figures/seismicwave_1layer_vz.mp4 differ
diff --git a/book/4_mechanical_waves/reflection_transmission.md b/book/4_mechanical_waves/reflection_transmission.md
index b64a95f..b7f348f 100644
--- a/book/4_mechanical_waves/reflection_transmission.md
+++ b/book/4_mechanical_waves/reflection_transmission.md
@@ -1,7 +1,13 @@
(subsec-wavetheory)=
# Reflection and transmission at boundaries
-In this section, we will consider what happens at a boundary between two layers with different mechanical properties. Then we can show which characteristics are responsible for, e.g., a reflection. To that end, we will focus on the basic physical and theoretical principles of reflection and transmission in this section. The basic equations describing wave motion in one direction will be used to derive an expression for the reflection and transmission coefficient at a boundary between two layers with different wave speeds and densities. In the previous section, the solution for the wave equation in the frequency domain was given, i.e.:
+In this section, we will consider what happens at a boundary between two layers with different mechanical properties. Then we can show which characteristics are responsible for, e.g., a reflection. To that end, we will focus on the basic physical and theoretical principles of reflection and transmission in this section. The basic equations describing wave motion in one direction will be used to derive an expression for the reflection and transmission coefficient at a boundary between two layers with different wave speeds and densities. So in the movie we showed before (see again below), we consider only the part of the wave that is travelling purely downwards: it was already noticed that part of the wave is reflected at the boundary and part of it is transmitted. How much this is, and what it depends on, will be worked out below.
+
+```{video} figures/seismicwave_1layer_vz.mp4
+:width: 80%
+```
+
+In the previous section, the solution for the wave equation in the frequency domain was given, i.e.:
$$
P(x,\omega) = S(\omega) \exp (\pm i\omega x/c).
diff --git a/book/4_mechanical_waves/snells_law.md b/book/4_mechanical_waves/snells_law.md
index bb0e864..4ab840f 100644
--- a/book/4_mechanical_waves/snells_law.md
+++ b/book/4_mechanical_waves/snells_law.md
@@ -116,7 +116,13 @@ $$ (eq:trefract)
We recognize this equation as the equation of a straight line when $t$ is considered as a function of the distance $x,$ the line along which we do our wavefield measurements.
-We have now derived the equations for the three rays, and we can plot their travel times as a function of distance $x$. This is done in {numref}`fig-raytimes`. This picture is important. When we measure seismic data in the field, the characteristics in this plot can most of the time be observed.
+We have now derived the equations for the three rays, and we can plot their travel times as a function of distance $x$. This is done in {numref}`fig-raytimes`. This picture is important. When we measure seismic data in the field, the characteristics in this plot can be observed most of the time.
+
+We considered travel times, but they are a result of a seismic wave field propagating in the subsurface. To get an idea of how a wave moves, we consider the movie below. A seismic source initiates a seismic wave at the top and the wave propagates outward (on a circle, for a homogeneous medium). Then later the wave encounters a boundary at which part of the wave is reflected and part of it is transmitted into the lower medium. There it travels with a higher speed; it can be seen that the wave travels faster along the boundary there. Then, part of it is diffracted back into the first medium (quite faint): that is the head wave, or critically refracted ray.
+
+```{video} figures/seismicwave_1layer_vz.mp4
+:width: 80%
+```
```{figure} figures/raytimes.png
:name: fig-raytimes
diff --git a/book/4_mechanical_waves/wave_equation_solutions.md b/book/4_mechanical_waves/wave_equation_solutions.md
new file mode 100644
index 0000000..459ad4d
--- /dev/null
+++ b/book/4_mechanical_waves/wave_equation_solutions.md
@@ -0,0 +1,48 @@
+# Solution(s) to the wave equation
+
+The solution to the wave equation for the pressure is:
+
+$$
+p(x,t) = s(t \pm x/c)
+$$ (eq:pxt)
+
+where $s(t)$ is some function. Note the dependency on space and time via $(t \pm x/c)$, which denotes that it is a travelling wave, and behaves as a direct wave as discussed in [](./snells_law.md). The sign in the argument depends on the direction the wave is travelling.
+
+Often, seismic responses are analyzed in terms of frequencies, i.e., Fourier spectra. The definition used here is:
+
+$$
+G(\omega) = \int_{-\infty}^{+\infty} g(t) \exp(- i \omega t) \, dt
+$$ (eq:fourier1)
+
+where the radial frequency $\omega$ is used, so $\omega = 2 \pi f$ where $f$ is the (ordinary) frequency. Using this convention, it is easy to show that differentiation with respect to time is equivalent to multiplication with $i\omega$ in the Fourier domain. When we transform the solution of the wave equation ({eq}`eq:pxt`) to the Fourier domain, we obtain:
+
+$$
+P(x,\omega) = S(\omega) \exp (\pm i\omega x/c).
+$$ (eq:Pxw)
+
+Note here that the time delay $x/c$ (in the time domain) becomes a *linear* phase in the Fourier domain, so the phase $\pm \omega x/c$ is linear as a function of $\omega$.
+
+In the above, we gave an expression for the pressure, but one can also derive an equivalent expression for the particle velocity $v_x$. For that purpose, we can use the equation of motion as expressed in {eq}`eq:eqmotion`, but then in its Fourier-transformed version, which is:
+
+$$
+V_x(x,\omega) = -\frac{1}{i\omega \rho} \frac{\partial P(x,\omega)}{\partial x}.
+$$
+
+When we substitute the solution for the pressure from above ({eq}`eq:Pxw`), we get for the negative sign:
+
+$$
+\begin{aligned}
+V_x(x,\omega) & = -\frac{1}{i \omega \rho}
+ S(\omega) \frac{- i \omega}{c} \exp(-i\omega x/c)
+\\
+ & = S(\omega) \frac{1}{\rho c} \exp(-i\omega x/c).
+\end{aligned}
+$$
+
+Notice that the particle velocity is a scaled version of the pressure:
+
+$$
+V_x(x,\omega) = \frac{P(x,\omega)}{\rho c}.
+$$
+
+The scaling factor is ($\rho c$), being called the *seismic impedance*.
diff --git a/book/_toc.yml b/book/_toc.yml
index 4c35deb..825c9c3 100644
--- a/book/_toc.yml
+++ b/book/_toc.yml
@@ -106,17 +106,22 @@ parts:
- file: 3_diffusion_fields/heat_equation_1d.md
- file: 3_diffusion_fields/em_diffusion.md
- file: 3_diffusion_fields/labs/week05_diffusion.md
+ - file: 3_diffusion_fields/labs/week06_em_diffusion.md
- caption: Mechanical waves
chapters:
- file: 4_mechanical_waves/intro.md
sections:
- file: 4_mechanical_waves/basic_notions.md
- file: 4_mechanical_waves/snells_law.md
+ - file: 4_mechanical_waves/exercise_4a.md
- file: 4_mechanical_waves/field_recordings.md
+ - file: 4_mechanical_waves/exercise_4b.md
- file: 4_mechanical_waves/fundamental_laws.md
- file: 4_mechanical_waves/acoustic_waves.md
+ - file: 4_mechanical_waves/wave_equation_solutions.md
- file: 4_mechanical_waves/reflection_transmission.md
- file: 4_mechanical_waves/elastic_waves.md
+ - file: 4_mechanical_waves/exercises.md
- file: 4_mechanical_waves/appendix_acoustic_1d.md
- file: 4_mechanical_waves/appendix_strain_p_s.md
- caption: Miscellaneous