diff --git a/book/0_overview/schedule.md b/book/0_overview/schedule.md
index 7552623..7847f84 100644
--- a/book/0_overview/schedule.md
+++ b/book/0_overview/schedule.md
@@ -1,3 +1,60 @@
 # Weekly schedule
 
-Click on the dropdown blocks below to find the schedule of each week's activities.
\ No newline at end of file
+Each unit runs over one week and consists of three lessons, numbered after the unit. Unit 1.4, for example, is made up of lessons 1.4.1, 1.4.2 and 1.4.3. The first two are lectures and the third is the practical, where the computer labs in this book are used.
+
+[My Timetable](https://mytimetable.tudelft.nl/) carries the authoritative schedule, with rooms and any changes.
+
+## Quarter 1
+
+| Unit | Topic | Lesson | Date | Time |
+| :--- | :--- | :--- | :--- | :--- |
+| 1.1 | Gradient and divergence | 1.1.1 | Tue 1 Sep 2026 | 08:45-10:30 |
+|  |  | 1.1.2 | Thu 3 Sep 2026 | 15:45-17:30 |
+|  |  | 1.1.3 | Fri 4 Sep 2026 | 13:45-15:30 |
+| 1.2 | Curl | 1.2.1 | Mon 7 Sep 2026 | 10:45-12:30 |
+|  |  | 1.2.2 | Tue 8 Sep 2026 | 13:45-15:30 |
+|  |  | 1.2.3 | Fri 11 Sep 2026 | 13:45-16:30 |
+| 1.3 | Potential fields: history and experiments | 1.3.1 | Mon 14 Sep 2026 | 10:45-12:30 |
+|  |  | 1.3.2 | Tue 15 Sep 2026 | 13:45-15:30 |
+|  |  | 1.3.3 | Fri 18 Sep 2026 | 13:45-15:30 |
+| 1.4 | Potential fields: gravity, magnetic field of the Earth | 1.4.1 | Mon 21 Sep 2026 | 10:45-12:30 |
+|  |  | 1.4.2 | Tue 22 Sep 2026 | 13:45-15:30 |
+|  |  | 1.4.3 | Fri 25 Sep 2026 | 13:45-16:30 |
+| 1.5 | Electric field. Diffusion fields: hot wire | 1.5.1 | Mon 28 Sep 2026 | 10:45-12:30 |
+|  |  | 1.5.2 | Tue 29 Sep 2026 | 13:45-15:30 |
+|  |  | 1.5.3 | Fri 2 Oct 2026 | 13:45-15:30 |
+| 1.6 | Diffusion fields: boundary conditions, heat in 2D and 3D | 1.6.1 | Mon 5 Oct 2026 | 10:45-12:30 |
+|  |  | 1.6.2 | Tue 6 Oct 2026 | 15:45-17:30 |
+|  |  | 1.6.3 | Fri 9 Oct 2026 | 13:45-16:30 |
+| 1.7 | Mechanical waves: strings, acoustic waves, 2D and 3D | 1.7.1 | Mon 12 Oct 2026 | 10:45-12:30 |
+|  |  | 1.7.2 | Tue 13 Oct 2026 | 13:45-15:30 |
+|  |  | 1.7.3 | Fri 16 Oct 2026 | 13:45-15:30 |
+| 1.8 | Mechanical waves: power flux | 1.8.1 | Mon 19 Oct 2026 | 10:45-12:30 |
+|  |  | 1.8.2 | Tue 20 Oct 2026 | 13:45-15:30 |
+|  |  | 1.8.3 | Fri 23 Oct 2026 | 13:45-16:30 |
+| 1.9 | Unsupervised study |  |  |  |
+| 1.10 | Midterm week: exam and discussion of solutions | 1.10.1 |  |  |
+|  |  | 1.10.2 |  |  |
+|  |  | 1.10.3 | Fri 6 Nov 2026 | 13:30-16:30 |
+
+## Quarter 2
+
+| Unit | Topic | Lesson | Date | Time |
+| :--- | :--- | :--- | :--- | :--- |
+| 2.1 | Electromagnetism: Maxwell's equations, plane waves, telegraph equation | 2.1.1 | Mon 9 Nov 2026 | 13:45-15:30 |
+|  |  | 2.1.2 | Tue 10 Nov 2026 | 13:45-15:30 |
+|  |  | 2.1.3 | Fri 13 Nov 2026 | 10:45-12:30 |
+| 2.2 | 3D waves, Poynting vector, polarisation, lossy and lossless media | 2.2.1 | Mon 16 Nov 2026 | 13:45-15:30 |
+|  |  | 2.2.2 | Tue 17 Nov 2026 | 13:45-15:30 |
+|  |  | 2.2.3 | Fri 20 Nov 2026 | 09:45-12:30 |
+| 2.3 | Reflection, transmission, refraction. Multi-layered media | 2.3.1 | Mon 23 Nov 2026 | 13:45-15:30 |
+|  |  | 2.3.2 | Tue 24 Nov 2026 | 13:45-15:30 |
+|  |  | 2.3.3 | Fri 27 Nov 2026 | 10:45-12:30 |
+| 2.4 | Trapped and surface waves. Phase and group velocity | 2.4.1 | Mon 30 Nov 2026 | 13:45-15:30 |
+|  |  | 2.4.2 | Tue 1 Dec 2026 | 08:45-10:30 |
+|  |  | 2.4.3 | Fri 4 Dec 2026 | 09:45-12:30 |
+| 2.8 | Unsupervised study |  |  |  |
+| 2.9 | Unsupervised study |  |  |  |
+| 2.10 | Resit exam |  | Wed 27 Jan 2027 | 13:30-16:30 |
+
+The exam is on Wed 16 Dec 2026, 13:30-16:30. The longer practicals, running to 16:30, end with a formative assessment.
diff --git a/book/1_gradient_divergence_curl/coord_sys.md b/book/1_gradient_divergence_curl/coord_sys.md
new file mode 100644
index 0000000..4d3ccca
--- /dev/null
+++ b/book/1_gradient_divergence_curl/coord_sys.md
@@ -0,0 +1,315 @@
+# Coordinate systems
+
+Three-dimensional space can be built from many frames of reference, but here it is built from the rectangular, cylindrical and spherical coordinate systems.
+
+- The **rectangular** (Cartesian) coordinate system is characterised by the three base vectors $\hat{\boldsymbol x},\hat{\boldsymbol y},\hat{\boldsymbol z}$ with coordinates and ranges $-\infty<x<\infty$, $-\infty<y<\infty$, $-\infty<z<\infty$, length elements $\mathrm{d}x,\mathrm{d}y,\mathrm{d}z$, and all three coordinate surfaces are planes.
+- The **cylindrical** coordinate system is characterised by the base vectors $\hat{\boldsymbol\varrho},\hat{\boldsymbol\phi},\hat{\boldsymbol z}$ with coordinates and ranges $0<\varrho<\infty$, $0<\phi<2\pi$, $-\infty<z<\infty$, length elements $\mathrm{d}\varrho,\ \varrho\,\mathrm{d}\phi,\ \mathrm{d}z$. The coordinate surfaces for constant $\phi$ and for constant $z$ are planes, while the surface of constant $\varrho$ has a cylindrical shape.
+- The **spherical** coordinate system is characterised by the base vectors $\hat{\boldsymbol r},\hat{\boldsymbol\phi},\hat{\boldsymbol\theta}$ with coordinates and ranges $0<r<\infty$, $0<\phi<2\pi$, $0<\theta<\pi$, length elements $\mathrm{d}r,\ r\sin(\theta)\,\mathrm{d}\phi,\ r\,\mathrm{d}\theta$. The coordinate surface for constant $r$ is a sphere, for constant $\phi$ a plane, and for constant $\theta$ a cone.
+
+:::{admonition} A note on the radial symbol
+:class: note
+The cylindrical radial coordinate is written here as $\varrho$ and the spherical one as $r$. Both measure a distance from an axis or from a point, but they are *not* the same quantity: $\varrho=\sqrt{x^2+y^2}$ while $r=\sqrt{x^2+y^2+z^2}$, so that $r=\sqrt{\varrho^2+z^2}$. In the figures below, and in some textbooks, the cylindrical radius is also drawn as $r$. Keeping the two apart from the start avoids a great deal of confusion later, when the gradient and the divergence are written out in both systems.
+:::
+
+Physical measurements take place in space-time, and the physical quantities that are described in classical physics are tensor functions of space and time. The time coordinate is indicated by the symbol $t$. The derivative of a differentiable function $f(\boldsymbol r,t)$ with respect to the time coordinate is given by
+
+$$
+\partial_t f(\boldsymbol r,t) = \frac{\partial f(\boldsymbol r,t)}{\partial t},
+$$
+
+where it is understood that the subscript $t$ is reserved for indicating time.
+
+## Cartesian reference frame
+
+Tensor calculus has been developed as the mathematical vehicle to describe macroscopic physical phenomena. It is a direct consequence of the postulate of continuum physics that all quantities occurring in the laws of macroscopic physics are geometric objects and therefore they are tensor quantities. Tensor calculus prescribes the rules for algebraic manipulations on tensors and how to apply the rules of differentiation and integration to tensors. The corresponding notation is known as the subscript notation. Because the subscript notation is not the generally used notation, here the vector notation is used.
+
+The Cartesian reference frame in three-dimensional Euclidean space is the set of base vectors $\{\hat{\boldsymbol x},\hat{\boldsymbol y},\hat{\boldsymbol z}\}$. The base vectors are of unit length each and in the given order they form a right-handed system. We can write the position vector $\boldsymbol r$ as the linear combination of the base vectors,
+
+$$
+\boldsymbol r = x\hat{\boldsymbol x} + y\hat{\boldsymbol y} + z\hat{\boldsymbol z}.
+$$
+
+It is customary in the Earth Sciences to let the positive $\hat{\boldsymbol z}$-axis point downward, as depicted in {numref}`fig-cartframe`.
+
+```{figure} figures/Cartframe.png
+:name: fig-cartframe
+:width: 60%
+
+The Cartesian reference frame, its base vectors, and a resulting vector to denote a point in space.
+```
+
+### Distance and the scalar product
+
+The distance between any two points with position vectors $\boldsymbol r$ and $\boldsymbol r'$ is given by the scalar product of the distance vector,
+
+$$
+d(\boldsymbol r,\boldsymbol r') = \left[(\boldsymbol r-\boldsymbol r')\cdot(\boldsymbol r-\boldsymbol r')\right]^{1/2} = |\boldsymbol r-\boldsymbol r'|,
+$$
+
+where the scalar product of two vectors $\boldsymbol a$ and $\boldsymbol b$ is given by
+
+$$
+\boldsymbol a\cdot\boldsymbol b = a_x b_x + a_y b_y + a_z b_z .
+$$
+
+The notation for the scalar product is a dot between the two vectors, for which reason it is also known as the dot product. The scalar product of two vectors multiplies the lengths of the two vectors and scales the result by the cosine of the angle $\phi$ between these vectors, hence
+
+$$
+\boldsymbol a\cdot\boldsymbol b = |\boldsymbol a||\boldsymbol b|\cos(\phi).
+$$
+
+### Translation and rotation
+
+Translation and rotation of the vector $\boldsymbol r$ are linear transformations, and the resulting vector $\boldsymbol r'$ is obtained from $\boldsymbol r$ by
+
+$$
+\boldsymbol r' = \hat{\boldsymbol e}_x(x-\lambda_x) + \hat{\boldsymbol e}_y(y-\lambda_y) + \hat{\boldsymbol e}_z(z-\lambda_z).
+$$
+
+Since distance should be invariant under rotation, only orthogonal tensors $\{\hat{\boldsymbol e}_x,\hat{\boldsymbol e}_y,\hat{\boldsymbol e}_z\}$ describe a rotation, and when these vectors have unit amplitude they span a base for the new frame of reference, in which case they are orthonormal vectors.
+
+```{figure} figures/rotCartframe.png
+:name: fig-rotcartframe
+:width: 60%
+
+The translation and rotation of the Cartesian reference frame.
+```
+
+### The vector product
+
+Any two non-parallel vectors span a plane, and the direction perpendicular to that plane is obtained by the so-called vector product of these two vectors. It is written as
+
+$$
+\boldsymbol a\times\boldsymbol b = \hat{\boldsymbol x}(a_y b_z - a_z b_y) + \hat{\boldsymbol y}(a_z b_x - a_x b_z) + \hat{\boldsymbol z}(a_x b_y - a_y b_x),
+$$
+
+and the notation is to use a cross between the vectors, for which reason it is also known as the cross product. The geometrical interpretation is that the new vector has a length equal to the area spanned by the two vectors $\boldsymbol a$ and $\boldsymbol b$, and is directed perpendicular to the plane containing the two vectors according to the right-hand rule, hence
+
+$$
+\boldsymbol c = \boldsymbol a\times\boldsymbol b = |\boldsymbol a||\boldsymbol b|\sin(\phi)\,\hat{\boldsymbol n},
+$$
+
+where $\hat{\boldsymbol n}$ is the unit vector in the direction of $\boldsymbol c$, as depicted in {numref}`fig-crossprod`. This implies that the cross-product rule does not commute, hence
+
+$$
+\boldsymbol a\times\boldsymbol b \ne \boldsymbol b\times\boldsymbol a,
+\qquad
+\boldsymbol a\times\boldsymbol b = -\,\boldsymbol b\times\boldsymbol a .
+$$
+
+```{figure} figures/crossprod.png
+:name: fig-crossprod
+:width: 35%
+
+The vector product of the two vectors $\boldsymbol a$ and $\boldsymbol b$ results in the vector $\boldsymbol c$ that is perpendicular to the plane containing the two vectors, according to the right-hand rule. The length of $\boldsymbol c$ is equal to the grey coloured area, which is the parallelogram spanned by the two vectors.
+```
+
+Because the vector product creates a coordinate system when applied to any two non-parallel vectors, it can also be used to create the base vectors of the right-handed Cartesian reference frame:
+
+$$
+\hat{\boldsymbol z} = \hat{\boldsymbol x}\times\hat{\boldsymbol y},
+\qquad
+\hat{\boldsymbol y} = \hat{\boldsymbol z}\times\hat{\boldsymbol x},
+\qquad
+\hat{\boldsymbol x} = \hat{\boldsymbol y}\times\hat{\boldsymbol z}.
+$$
+
+## Curvilinear coordinate systems
+
+Apart from the Cartesian reference frame, which is a rectangular coordinate system, two coordinate systems are often employed that use curved base vectors. One is the cylindrical coordinate system and the other is the spherical coordinate system. The point shown in {numref}`fig-cartframe` is depicted in cylindrical and spherical coordinate systems in {numref}`fig-cylsphere`.
+
+In a **cylindrical** coordinate system a point in 3D space is specified by its coordinates $(\varrho,\phi,z)$, where $\phi$ is measured from the $x$-axis perpendicular to the $(x,z)$-plane, which is then a plane of constant $\phi$. Hence, the curved base vector $\hat{\boldsymbol\phi}$ lies in the $(x,y)$-plane. The surface of constant $\varrho$ is the boundary surface of the cylinder, which is the cylindrical $(\phi,z)$-surface, and the plane of constant $z$ is the circular cross-section of the cylinder, which is the $(\varrho,\phi)$-plane.
+
+In a **spherical** coordinate system a point in 3D space is specified by its coordinates $(r,\phi,\theta)$, where $\phi$ is measured from the $x$-axis to the projection of the point onto the $(x,y)$-plane. The plane of constant $\phi$ is the $(r,\theta)$-plane. The surface of constant $r$ is the boundary surface of the ball with radius $r$, which is the spherical $(\phi,\theta)$-surface, and the surface of constant $\theta$ is the conical cross-section of the ball with radius $r$, whose intersection with the spherical surface is a circle; this is the $(r,\phi)$-surface.
+
+From these descriptions the coordinate values are easily transformed to the rectangular coordinate system from the cylindrical and spherical coordinate systems and vice versa, as shown in {numref}`tab-coord-transform`. Note that the radial vectors $\hat{\boldsymbol\varrho}$ and $\hat{\boldsymbol r}$ in the cylindrical and spherical coordinate systems are **not** the same.
+
+```{figure} figures/cylsphere.png
+:name: fig-cylsphere
+:width: 100%
+
+The Cartesian reference frame together with the cylindrical (left) and spherical (right) coordinate systems, their base vectors, and a resulting vector to denote a point in space. In the figure the cylindrical radius is drawn as $r$; in the text it is written $\varrho$.
+```
+
+```{list-table} Spatial coordinate transformations.
+:header-rows: 1
+:name: tab-coord-transform
+
+* - Cylindrical to rectangular
+  - Spherical to rectangular
+* - $x = \varrho\cos(\phi)$
+  - $x = r\cos(\phi)\sin(\theta)$
+* - $y = \varrho\sin(\phi)$
+  - $y = r\sin(\phi)\sin(\theta)$
+* - $z = z$
+  - $z = r\cos(\theta)$
+* - **Rectangular to cylindrical**
+  - **Rectangular to spherical**
+* - $\varrho = \sqrt{x^2+y^2},\quad \varrho\ge 0$
+  - $r = \sqrt{x^2+y^2+z^2},\quad r\ge 0$
+* - $\phi = \mathrm{atan}(y/x),\quad -\pi<\phi\le\pi$
+  - $\phi = \mathrm{atan}(y/x),\quad -\pi<\phi\le\pi$
+* - $z = z$
+  - $\theta = \mathrm{atan}\!\left(\dfrac{\sqrt{x^2+y^2}}{z}\right),\quad 0\le\theta\le\pi$
+```
+
+To cover the full range $-\pi<\phi\le\pi$, the two-argument arctangent $\mathrm{atan2}(y,x)$ must be used, because $\mathrm{atan}(y/x)$ alone cannot distinguish the second and third quadrants from the fourth and first. The same holds for $\theta$, where the sign of $z$ decides between the upper and lower hemisphere.
+
+The rectangular components $v_x,v_y,v_z$ of a vector function $\boldsymbol v$ at a point $\boldsymbol r$ with spherical components $v_r,v_\phi,v_\theta$ are given by
+
+$$
+\begin{aligned}
+v_x &= v_r\cos(\phi)\sin(\theta) - v_\phi\sin(\phi) + v_\theta\cos(\phi)\cos(\theta), \\
+v_y &= v_r\sin(\phi)\sin(\theta) + v_\phi\cos(\phi) + v_\theta\sin(\phi)\cos(\theta), \\
+v_z &= v_r\cos(\theta) - v_\theta\sin(\theta).
+\end{aligned}
+$$ (eq:sph2cart-components)
+
+In these relations the direction cosines occur from projections of spherical coordinates to rectangular coordinates. Projections are described by scalar products, hence the direction cosines are the results of the scalar products of the base vectors in two different coordinate systems. All direction cosines between the three coordinate systems are given in {numref}`tab-direction-cosines`. The upper-right diagonal entries are left empty and return as an exercise below.
+
+```{list-table} Direction cosines from base vector scalar products. Columns 1–3 are the Cartesian base vectors, columns 4–6 the cylindrical ones, and columns 7–9 the spherical ones; the rows follow the same order. Entries above the diagonal are left open (see the exercises).
+:header-rows: 1
+:name: tab-direction-cosines
+
+* -
+  - $\hat{\boldsymbol x}$
+  - $\hat{\boldsymbol y}$
+  - $\hat{\boldsymbol z}$
+  - $\hat{\boldsymbol\varrho}$
+  - $\hat{\boldsymbol\phi}$
+  - $\hat{\boldsymbol z}$
+  - $\hat{\boldsymbol r}$
+  - $\hat{\boldsymbol\phi}$
+  - $\hat{\boldsymbol\theta}$
+* - $\hat{\boldsymbol x}$
+  - $1$
+  - –
+  - –
+  - –
+  - –
+  - –
+  - –
+  - –
+  - –
+* - $\hat{\boldsymbol y}$
+  - $0$
+  - $1$
+  - –
+  - –
+  - –
+  - –
+  - –
+  - –
+  - –
+* - $\hat{\boldsymbol z}$
+  - $0$
+  - $0$
+  - $1$
+  - –
+  - –
+  - –
+  - –
+  - –
+  - –
+* - $\hat{\boldsymbol\varrho}$
+  - $\cos(\phi)$
+  - $\sin(\phi)$
+  - $0$
+  - $1$
+  - –
+  - –
+  - –
+  - –
+  - –
+* - $\hat{\boldsymbol\phi}$
+  - $-\sin(\phi)$
+  - $\cos(\phi)$
+  - $0$
+  - $0$
+  - $1$
+  - –
+  - –
+  - –
+  - –
+* - $\hat{\boldsymbol z}$
+  - $0$
+  - $0$
+  - $1$
+  - $0$
+  - $0$
+  - $1$
+  - –
+  - –
+  - –
+* - $\hat{\boldsymbol r}$
+  - $\cos(\phi)\sin(\theta)$
+  - $\sin(\phi)\sin(\theta)$
+  - $\cos(\theta)$
+  - $\sin(\theta)$
+  - $0$
+  - $\cos(\theta)$
+  - $1$
+  - –
+  - –
+* - $\hat{\boldsymbol\phi}$
+  - $-\sin(\phi)$
+  - $\cos(\phi)$
+  - $0$
+  - $0$
+  - $1$
+  - $0$
+  - $0$
+  - $1$
+  - –
+* - $\hat{\boldsymbol\theta}$
+  - $\cos(\phi)\cos(\theta)$
+  - $\sin(\phi)\cos(\theta)$
+  - $-\sin(\theta)$
+  - $\cos(\theta)$
+  - $0$
+  - $-\sin(\theta)$
+  - $0$
+  - $0$
+  - $1$
+```
+
+The vector $\boldsymbol v$ is expressed in spherical coordinates as
+
+$$
+\boldsymbol v = v_r\hat{\boldsymbol r} + v_\phi\hat{\boldsymbol\phi} + v_\theta\hat{\boldsymbol\theta}.
+$$
+
+## Exercises
+
+1. What is the length of the vector function $\boldsymbol v$?
+2. Show that $\hat{\boldsymbol x} = \partial\boldsymbol r/\partial x$.
+3. A two-dimensional rectangle in the plane $z=0$ is spanned by the two vectors $\boldsymbol d_x = d_x\hat{\boldsymbol x}$ and $\boldsymbol d_y = d_y\hat{\boldsymbol y}$, where $d_x$ and $d_y$ are the lengths in the $x$- and $y$-directions, respectively, and $\hat{\boldsymbol n}$ is the unit normal on the rectangle pointing in the positive $z$-direction. Give a geometrical interpretation of the product $\hat{\boldsymbol n}\cdot(\boldsymbol d_x\times\boldsymbol d_y)$.
+4. A three-dimensional rectangle, a brick, with dimensions $d_x,d_y,d_z$ is spanned by the three vectors $\boldsymbol d_x = d_x\hat{\boldsymbol x}$, $\boldsymbol d_y = d_y\hat{\boldsymbol y}$ and $\boldsymbol d_z = d_z\hat{\boldsymbol z}$. Give a geometrical interpretation of the product $\boldsymbol d_x\cdot(\boldsymbol d_y\times\boldsymbol d_z)$.
+5. Consider a smooth surface $\mathbb{S}$ with unique unit normal vector $\hat{\boldsymbol n}$. Show that any vector quantity $\boldsymbol H$ can be composed as $\boldsymbol H = (\hat{\boldsymbol n}\cdot\boldsymbol H)\hat{\boldsymbol n} + (\hat{\boldsymbol n}\times\boldsymbol H)\times\hat{\boldsymbol n}$ and give a geometric interpretation of the two terms.
+6. A smooth surface $\mathbb{S}$ has a unique unit normal vector $\hat{\boldsymbol n}$. Decompose the gradient operator into a part tangential and a part normal to the surface $\mathbb{S}$.
+7. Use the same reasoning as in the exercise above and find $\hat{\boldsymbol r} = A_r\,\partial\boldsymbol r/\partial r$, $\hat{\boldsymbol\phi} = A_\phi\,\partial\boldsymbol r/\partial\phi$, and $\hat{\boldsymbol\theta} = A_\theta\,\partial\boldsymbol r/\partial\theta$. Carry out the differentiations and find the coefficients $(A_r,A_\phi,A_\theta)$ by requiring that these vectors have unit length. These results should lead to an expression for the unit vectors of the spherical coordinate system in Cartesian coordinates.
+8. The upper right-hand side of the matrix of {numref}`tab-direction-cosines` is left empty. The $3\times3$ block in the bottom-left corner of the table relates the unit vectors of the spherical reference frame to those of the Cartesian reference frame. Suppose we write this as
+
+    $$
+    \left(\begin{array}{c} \hat{\boldsymbol r} \\ \hat{\boldsymbol\phi} \\ \hat{\boldsymbol\theta}\end{array}\right)
+    = \mathsf{R}\left(\begin{array}{c} \hat{\boldsymbol x} \\ \hat{\boldsymbol y} \\ \hat{\boldsymbol z}\end{array}\right),
+    $$ (eq:sphCar)
+
+    $$
+    \mathsf{R} = \left(\begin{array}{ccc}
+    \cos(\phi)\sin(\theta) & \sin(\phi)\sin(\theta) & \cos(\theta) \\
+    -\sin(\phi) & \cos(\phi) & 0 \\
+    \cos(\phi)\cos(\theta) & \sin(\phi)\cos(\theta) & -\sin(\theta)
+    \end{array}\right).
+    $$
+
+    Show that
+
+    $$
+    \left(\begin{array}{c} \hat{\boldsymbol x} \\ \hat{\boldsymbol y} \\ \hat{\boldsymbol z}\end{array}\right)
+    = \mathsf{R}^{t}\left(\begin{array}{c} \hat{\boldsymbol r} \\ \hat{\boldsymbol\phi} \\ \hat{\boldsymbol\theta}\end{array}\right),
+    $$
+
+    where $\mathsf{R}^{t}$ denotes the transpose of $\mathsf{R}$. Then show that $\mathsf{R}\mathsf{R}^{t} = \mathsf{I}$, where $\mathsf{I}$ is the $3\times3$ unit matrix, implying that the transpose of $\mathsf{R}$ is equal to its inverse. This follows from the fact that $\mathsf{R}$ is a matrix of orthonormal vectors and is therefore an orthonormal matrix.
+9. Give explicit expressions for the base vectors of the spherical reference frame by carrying out the matrix-vector multiplication of {eq}`eq:sphCar`. Notice that all three base vectors depend on both angles. This implies that the base vectors are not constant but depend on position! The spherical reference frame is not an inertial frame.
+10. Carry out the following differentiations: $\partial\hat{\boldsymbol r}/\partial\theta$, $\partial\hat{\boldsymbol\theta}/\partial\theta$, $\partial\hat{\boldsymbol\phi}/\partial\theta$ and $\partial\hat{\boldsymbol r}/\partial\phi$, $\partial\hat{\boldsymbol\theta}/\partial\phi$, $\partial\hat{\boldsymbol\phi}/\partial\phi$.
diff --git a/book/1_gradient_divergence_curl/curl.md b/book/1_gradient_divergence_curl/curl.md
new file mode 100644
index 0000000..d73a939
--- /dev/null
+++ b/book/1_gradient_divergence_curl/curl.md
@@ -0,0 +1,106 @@
+# Curl
+
+The curl is a peculiar operator because it turns a vector field into another vector field. Often a determinant expression is used, but we can give it in terms of even and odd permutations of the triplet $(x,y,z)$, which is easy to remember. Let $\boldsymbol v(\boldsymbol r,t)$ be a continuously differentiable space-time function. Then
+
+$$
+\nabla\times\boldsymbol v = \hat{\boldsymbol x}(\partial_y v_z - \partial_z v_y) + \hat{\boldsymbol y}(\partial_z v_x - \partial_x v_z) + \hat{\boldsymbol z}(\partial_x v_y - \partial_y v_x).
+$$ (eq:curl)
+
+Notice that the three terms have the triplet pairs $(xyz)-(xzy)$, $(yzx)-(yxz)$, and $(zxy)-(zyx)$ in (1) direction, (2) differentiation, and (3) vector component. The first in each pair is an even permutation of $(xyz)$ and the second is an odd permutation. It can be seen that in each term the direction of the resulting vector $\boldsymbol w = \nabla\times\boldsymbol v$ is perpendicular both to the direction of differentiation and to the component of $\boldsymbol v$ that is being differentiated. The curl is interpreted as the rate of circulation of the vector field $\boldsymbol v$.
+
+## The circulation integral
+
+Let us first look at the situation of an elementary area with centre point $\boldsymbol r_p = \left(\tfrac{1}{2}\mathrm{d}x,\ \tfrac{1}{2}\mathrm{d}y,\ \tfrac{1}{2}\mathrm{d}z\right)$. The circulation integral for the path from $y=0$ to $y=\mathrm{d}y$ at $x=\tfrac{1}{2}\mathrm{d}x$ and $z=0$ is given by
+
+$$
+\int_{y=0}^{\mathrm{d}y} v_y\!\left(\tfrac{1}{2}\mathrm{d}x,\,y,\,0\right)\mathrm{d}l ,
+$$
+
+and completing the circulation for all four sides of the elementary area as shown in {numref}`fig-circulation`, we find
+
+```{figure} figures/circulation.png
+:name: fig-circulation
+:width: 55%
+
+Net circulation integral on a rectangle around the point $\boldsymbol r_p$.
+```
+
+$$
+\begin{aligned}
+\oint_{\boldsymbol r}\boldsymbol\tau\cdot\boldsymbol v\,\mathrm{d}l
+&= \int_{y=0}^{\mathrm{d}y} v_y\!\left(\tfrac{1}{2}\mathrm{d}x,\,y,\,0\right)\mathrm{d}l
+ + \int_{z=0}^{\mathrm{d}z} v_z\!\left(\tfrac{1}{2}\mathrm{d}x,\,\mathrm{d}y,\,z\right)\mathrm{d}l \\
+&\quad - \int_{y=0}^{\mathrm{d}y} v_y\!\left(\tfrac{1}{2}\mathrm{d}x,\,y,\,\mathrm{d}z\right)\mathrm{d}l
+ - \int_{z=0}^{\mathrm{d}z} v_z\!\left(\tfrac{1}{2}\mathrm{d}x,\,0,\,z\right)\mathrm{d}l ,
+\end{aligned}
+$$ (eq:circint)
+
+where the symbol $\oint$ stands for a closed loop integral, $\boldsymbol\tau$ is the local unit tangent vector along the edges of $\mathbb{S}$, and $l$ is the arclength along the path.
+
+Earlier we introduced the Taylor series expansion for a function in one variable, and we give it here for each vector component of $\boldsymbol v$ around the point $\boldsymbol r_p$,
+
+$$
+\begin{aligned}
+v_x(\boldsymbol r) &= v_x(\boldsymbol r_p) + \partial_x v_x\left(x-\tfrac{1}{2}\mathrm{d}x\right) + \partial_y v_x\left(y-\tfrac{1}{2}\mathrm{d}y\right) + \partial_z v_x\left(z-\tfrac{1}{2}\mathrm{d}z\right) + \text{higher order terms}, \\
+v_y(\boldsymbol r) &= v_y(\boldsymbol r_p) + \partial_x v_y\left(x-\tfrac{1}{2}\mathrm{d}x\right) + \partial_y v_y\left(y-\tfrac{1}{2}\mathrm{d}y\right) + \partial_z v_y\left(z-\tfrac{1}{2}\mathrm{d}z\right) + \text{higher order terms}, \\
+v_z(\boldsymbol r) &= v_z(\boldsymbol r_p) + \partial_x v_z\left(x-\tfrac{1}{2}\mathrm{d}x\right) + \partial_y v_z\left(y-\tfrac{1}{2}\mathrm{d}y\right) + \partial_z v_z\left(z-\tfrac{1}{2}\mathrm{d}z\right) + \text{higher order terms},
+\end{aligned}
+$$
+
+which is the Taylor series expansion to first order. Substituting these results in {eq}`eq:circint` leads to
+
+$$
+\oint_{\boldsymbol r}\boldsymbol\tau\cdot\boldsymbol v\,\mathrm{d}l = (\partial_y v_z - \partial_z v_y)\,\mathrm{d}y\,\mathrm{d}z + \text{higher order terms}.
+$$
+
+Similarly, the net circulation integral for the elementary surface area perpendicular to the $\hat{\boldsymbol y}$-axis is $(\partial_z v_x - \partial_x v_z)\,\mathrm{d}x\,\mathrm{d}z$, and the net circulation integral for the elementary surface area perpendicular to the $\hat{\boldsymbol z}$-axis is $(\partial_x v_y - \partial_y v_x)\,\mathrm{d}x\,\mathrm{d}y$.
+
+For this reason, the physical interpretation of the curl is the net circulation per unit of surface area at the point $\boldsymbol r_p$, and it can be expressed as
+
+$$
+\hat{\boldsymbol n}\cdot(\nabla\times\boldsymbol v) = \lim_{S\rightarrow 0}\frac{\oint_{\boldsymbol r}\boldsymbol\tau\cdot\boldsymbol v\,\mathrm{d}l}{A},
+$$ (eq:circ1)
+
+in which
+
+$$
+A = \int_{\boldsymbol r\in\mathbb{S}}\mathrm{d}S
+$$
+
+is the surface area of the surface $\mathbb{S}$.
+
+:::{admonition} The curl in words
+:class: tip
+The curl of a vector field is the net circulation of that field per unit of surface area, and it points along the normal of the surface for which that circulation is largest.
+:::
+
+## Stokes' integral theorem
+
+The unit normal $\hat{\boldsymbol n}$ in {eq}`eq:circ1` is the normal vector to the surface area $\mathbb{S}$ and is oriented to the side of advance of a right-hand screw as it is turned in the direction of $\boldsymbol\tau$ around the boundary of $\mathbb{S}$. This leads to Stokes' integral theorem, given by
+
+$$
+\oint_{\boldsymbol r}\boldsymbol\tau\cdot\boldsymbol v\,\mathrm{d}l = \int_{\boldsymbol r\in\mathbb{S}}\hat{\boldsymbol n}\cdot(\nabla\times\boldsymbol v)\,\mathrm{d}S .
+$$ (eq:stokes)
+
+## Exercises
+
+1. Use Stokes' theorem to evaluate $\oint_{\boldsymbol r}\boldsymbol\tau\,\mathrm{d}l$, where $\boldsymbol\tau$ is the unit tangent along the closed boundary of the area $\mathbb{S}$. The integration runs in the direction of circulation that forms a right-handed system with the unit normal vector on $\mathbb{S}$.
+2. Show that when $\boldsymbol v(\boldsymbol r) = \boldsymbol a\,p(\boldsymbol r)$, where $\boldsymbol a$ is an arbitrary constant vector and $p(\boldsymbol r)$ is a continuously differentiable scalar function, Stokes' integral theorem gives
+
+    $$
+    \oint_{\boldsymbol r}\boldsymbol\tau\,p\,\mathrm{d}l = \int_{\boldsymbol r\in\mathbb{S}}(\hat{\boldsymbol n}\times\nabla)p\,\mathrm{d}S,
+    $$
+
+    which is Stokes' theorem for the gradient.
+3. Show that when $\boldsymbol v(\boldsymbol r) = \boldsymbol a\times\boldsymbol w(\boldsymbol r)$, where $\boldsymbol a$ is an arbitrary constant vector and $\boldsymbol w(\boldsymbol r)$ is a continuously differentiable vector function, Stokes' integral theorem gives
+
+    $$
+    \oint_{\boldsymbol r}\boldsymbol\tau\times\boldsymbol w\,\mathrm{d}l = \int_{\boldsymbol r\in\mathbb{S}}(\hat{\boldsymbol n}\times\nabla)\times\boldsymbol w\,\mathrm{d}S .
+    $$
+4. Ampère's law is given by $\nabla\times\boldsymbol H = \boldsymbol J$, which states that the electric current is equal to the curl of the magnetic field. Convert this to integral form using Stokes' theorem, by using a flat surface $\mathbb{S}$ and choosing a unit normal vector $\hat{\boldsymbol n}$ on $\mathbb{S}$. You should find
+
+    $$
+    \oint_{\boldsymbol r}\boldsymbol H\cdot\boldsymbol\tau\,\mathrm{d}l = \int_{\boldsymbol r\in\mathbb{S}}\hat{\boldsymbol n}\cdot\boldsymbol J\,\mathrm{d}S .
+    $$
+
+    Show in which direction the unit tangent vector $\boldsymbol\tau$ is circulating at the boundary of $\mathbb{S}$, and in which direction through the surface a positive current should run.
diff --git a/book/1_gradient_divergence_curl/divergence.md b/book/1_gradient_divergence_curl/divergence.md
new file mode 100644
index 0000000..1e97bf2
--- /dev/null
+++ b/book/1_gradient_divergence_curl/divergence.md
@@ -0,0 +1,250 @@
+# Divergence of a vector field
+
+Now we introduce a vector field $\boldsymbol v(x,y,z,t)$ and consider the field inside a rectangular volume with sides $\mathrm{d}x,\mathrm{d}y,\mathrm{d}z$ that we can shrink to an infinitesimal cube. Let the left side be located at position $x$ and the right side at $x+\mathrm{d}x$.
+
+If we consider $\boldsymbol v$ to be a fluid flow vector, we understand that the component $v_x$ is the flow in the $x$-direction, the component $v_y$ is the flow in the $y$-direction, and the component $v_z$ is the flow in the $z$-direction. The flow in the $x$-direction can flow in or out of the cube through faces perpendicular to the $x$-axis. Hence, the outward flow through the surface perpendicular to the $x$-axis at $x+\mathrm{d}x$ is given by $v_x(x+\mathrm{d}x,y,z,t)\,\mathrm{d}y\,\mathrm{d}z$, because $v_x(x+\mathrm{d}x,y,z,t)$ is the component of the flow perpendicular to that surface and $\mathrm{d}y\,\mathrm{d}z$ is the area of the surface through which it flows.
+
+Then it is clear that the outward flow through the left side is $-v_x(x,y,z,t)\,\mathrm{d}y\,\mathrm{d}z$, and the minus sign comes from the fact that the outward unit normal there is in the $(-x)$-direction. This means that the total outward flow in the $x$-direction is given by
+
+$$
+\left[v_x(x+\mathrm{d}x,y,z,t) - v_x(x,y,z,t)\right]\mathrm{d}y\,\mathrm{d}z = \partial_x v_x\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z = \partial_x v_x\,\mathrm{d}V .
+$$
+
+We can combine this with the outward flow in the $y$- and $z$-directions to obtain $(\partial_x v_x + \partial_y v_y + \partial_z v_z)\,\mathrm{d}V$. This can be written as
+
+$$
+\mathrm{d}\Phi = (\nabla\cdot\boldsymbol v)\,\mathrm{d}V,
+$$ (eq:flux)
+
+where $\Phi$ is the total outward flux.
+
+## Gauss' integral theorem
+
+On the other hand, the outward flow can be generally described as $(\boldsymbol v\cdot\hat{\boldsymbol n})\hat{\boldsymbol n}$, because it describes the sum of normal components directed along the outward unit normal vector on the surface. The volume of the flow through the surface per unit time is then given by
+
+$$
+\Phi = \int_{\mathbb{S}}(\boldsymbol v\cdot\hat{\boldsymbol n})\,\mathrm{d}S .
+$$ (eq:vflux)
+
+Note that the single integration symbol represents here a surface, or two-dimensional, integral, which is short-hand notation for the double integration that must be carried out. This is indicated by the integration surface $\mathbb{S}$ as a subscript to the integral, so that no confusion occurs about the kind of integration that is meant.
+
+This expression defines the **flux** $\Phi$ of the vector field $\boldsymbol v$ through the surface $\mathbb{S}$. The definition of a flux is not restricted to the flow of fluids. A flux can be computed for any vector field, such as the flux of the electric field through a surface.
+
+The statements in {eq}`eq:flux` and {eq}`eq:vflux` describe the same flow and should be equal. Let us write {eq}`eq:flux` as a total volume derivative,
+
+$$
+\frac{\mathrm{d}\Phi}{\mathrm{d}V} = \nabla\cdot\boldsymbol v,
+$$
+
+and integrate it over the volume; we find
+
+$$
+\Phi = \int_{\mathbb{D}}(\nabla\cdot\boldsymbol v)\,\mathrm{d}V .
+$$
+
+Note that the single integration symbol represents here a volume, or three-dimensional, integral, which is short-hand notation for the triple integration that must be carried out. This is indicated by the integration volume $\mathbb{D}$ as a subscript to the integral, so that no confusion occurs about the kind of integration that is meant. We find then
+
+$$
+\int_{\mathbb{D}}(\nabla\cdot\boldsymbol v)\,\mathrm{d}V = \int_{\mathbb{S}}(\boldsymbol v\cdot\hat{\boldsymbol n})\,\mathrm{d}S .
+$$ (eq:gauss)
+
+We have just derived Gauss' integral theorem, or divergence theorem, and we found the physical interpretation of the divergence of a vector field:
+
+:::{admonition} The divergence in words
+:class: tip
+The divergence of a vector field is the outward flux of the vector field per unit volume.
+:::
+
+## Incompressible flow
+
+For fluid flow under constant density and away from locations that have generated the flow, sometimes the approximation of the fluid to be incompressible can be made. Let the source or sink be located at $\boldsymbol r=\boldsymbol 0$. When the fluid is incompressible there cannot be any net outward flow, because what goes out on one side must come in on the other; otherwise the fluid would be compressible. We can state that
+
+$$
+\nabla\cdot\boldsymbol v = 0, \qquad \text{for } r\ne 0 .
+$$
+
+In a constant-density environment the flow field must be radial away from the source. Assuming the source is a volume injection at $\boldsymbol r=\boldsymbol 0$, we can write the field as
+
+$$
+\boldsymbol v = f(r)\,\boldsymbol r .
+$$
+
+The function $f(r)$ depends only on the radial distance $r=\sqrt{x^2+y^2+z^2}$. We can find $f(r)$ by requiring that the flow is divergence free,
+
+$$
+\nabla\cdot\boldsymbol v = 3f(r) + \boldsymbol r\cdot\nabla f(r) = 0, \qquad \text{for } r\ne 0 .
+$$
+
+This gives
+
+$$
+f(r) = \frac{A}{r^3},
+$$ (eq:NF)
+
+where $A$ is a constant that is determined by the source or sink.
+
+## Electric potential, electric field, and electric current
+
+The result above for fluid flow has an equivalent for the electric field and electric potential, which combines gradient and divergence. We start with the fact that for a constant current that is injected in a medium with constant resistivity $\rho$, we can write, for points away from the current injection point, that the electric field is divergence free,
+
+$$
+\nabla\cdot\boldsymbol E = 0 .
+$$ (eq:divE)
+
+The electric field can be written in terms of the gradient of an electric potential $V$ as
+
+$$
+\boldsymbol E(x,y,z) = -\nabla V(x,y,z).
+$$ (eq:EgradV)
+
+When we substitute this relation in {eq}`eq:divE`, we find the equation for the electric potential $V$ given by
+
+$$
+\nabla\cdot\nabla V = 0,
+$$ (eq:laplV)
+
+which written out in full gives
+
+$$
+\left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}\right)V = 0 .
+$$
+
+The divergence of a gradient is known as the Laplacian and is a scalar second-order differential operator. The solution to such an equation, where the potential originates from a point source in the origin, is well known and can be written as
+
+$$
+V(x,y,z) = \frac{A}{r} + B,
+$$ (eq:pot)
+
+where the distance to the source is given by $r=\sqrt{x^2+y^2+z^2}$ and the constant $B$ can be taken zero by requiring that the potential goes to zero infinitely far away from the location of the source of the potential field.
+
+To find the constant $A$ we need to find out the influence of the source we use to generate the electric potential. To do that we need the local version of Ohm's law, given by
+
+$$
+\boldsymbol J = \rho^{-1}\boldsymbol E,
+$$ (eq:bJ)
+
+in which $\boldsymbol J$ is the electric current density in A/m$^2$. Because a static electric field can be written as the gradient of the electric potential, this equation can be written as
+
+$$
+\boldsymbol J = -\rho^{-1}\nabla V .
+$$ (eq:Ohm)
+
+To find the constant $A$ we can think of a point at the origin where a current is injected. Because the spatial dependency of the potential is radial distance, the current will be distributed evenly in all directions. Therefore the surface of constant potential is a spherical surface, as can be seen from {eq}`eq:pot` as well when we find all the points where $V$ is constant. Only the current that is in the direction of the outward unit normal on any spherical surface can leave that surface, and hence we can write the total electric current leaving that surface (and hence going into the earth) as $I = 4\pi r^2\,\hat{\boldsymbol n}\cdot\boldsymbol J$, where $\hat{\boldsymbol n}$ is the outward directed unit normal on any spherical surface of constant potential and $4\pi r^2$ is the surface area of the sphere with radius $r$. We can use {eq}`eq:Ohm` to express the electric current density in terms of the electric potential and substitute the solution of the potential given in {eq}`eq:pot` to find
+
+$$
+I = 4\pi r^2\,\hat{\boldsymbol n}\cdot\boldsymbol J
+  = -4\pi r^2 A\,\hat{\boldsymbol n}\cdot\nabla\frac{1}{\rho r}
+  = 4\pi r^2 A\frac{\hat{\boldsymbol n}\cdot\boldsymbol r}{\rho r^3}
+  = \frac{4\pi A}{\rho},
+$$ (eq:Icur)
+
+because $\hat{\boldsymbol n}=\boldsymbol r/r$ and $\hat{\boldsymbol n}\cdot\hat{\boldsymbol n}=1$. Therefore
+
+$$
+A = \frac{\rho I}{4\pi},
+$$
+
+and the electric potential is given by
+
+$$
+V(x,y,z) = \frac{\rho I}{4\pi r},
+$$ (eq:dcVfull)
+
+where the point of the current injection must be avoided, so $r>0$. The total current that is injected into the ground times the total resistance equals the electric potential, which is Ohm's law. The total resistance is given by the electric resistivity $\rho$ divided by $4\pi$ times the radial distance from the current injection point.
+
+Now suppose there is a surface at $z=0$ between non-conductive air and the conductive subsurface, and the injection point is at the surface $z=0$. In that case the current can only go into the ground below the surface, hence for $z>0$, and the relevant surface area is $2\pi r^2$, because the current is now distributed over the surface area of half a sphere. Therefore, anywhere in the half-space $z>0$ the electric potential is given by
+
+$$
+V(x,y,z) = \frac{\rho I}{2\pi r},
+$$ (eq:dcVhf)
+
+where again $r>0$.
+
+You see that the electric potential depends on the value of the electric resistivity even though it does not occur in {eq}`eq:pot`. This is because the electric potential depends on the current strength, and that in turn depends on the resistivity of the ground through which this current must flow. When it is a constant it is merely a scaling parameter, but when it is a function of position it can become a complicated relation that must be found numerically.
+
+### Two electrodes at the surface
+
+{numref}`fig-dcpoth` and {numref}`fig-dcpotv` show a plot of the electric potential $V(x,y,z)$, for $z=0$ and for $y=0,\ z>0$ respectively. The arrows in the plots indicate the vector directions of the electric current. For this configuration we have the electric potential given by
+
+$$
+\begin{aligned}
+V(x,y,z) &= V(x-a/2,y,z) - V(x+a/2,y,z), \\
+V(x,y,z) &= \frac{\rho I}{2\pi}\left(\frac{1}{\sqrt{(x-a/2)^2+y^2+z^2}} - \frac{1}{\sqrt{(x+a/2)^2+y^2+z^2}}\right),
+\end{aligned}
+$$ (eq:Vpdp)
+
+where the point $x=a/2$ is the point of current injection (current goes into the ground, also known as source) and the point $x=-a/2$ is the current extraction point (current goes out of the ground, also known as sink), for which reason the potential related to that location is negative.
+
+```{figure} figures/dcpoth.png
+:name: fig-dcpoth
+:width: 75%
+
+Electric potential difference and electric current density vectors on the ground surface $z=0$, with two electrodes at $x=-a/2$ and $x=a/2$. Distances are normalised to the electrode spacing $a$.
+```
+
+```{figure} figures/dcpotv.png
+:name: fig-dcpotv
+:width: 85%
+
+Electric potential difference and electric current density vectors in the vertical cross-section $y=0,\ z>0$, with two electrodes at $x=-a/2$ and $x=a/2$. Distances are normalised to the electrode spacing $a$.
+```
+
+To make the current run in the subsurface, the two points must be connected to a current source above the ground through an electronically controlled connection with a battery or other charge-storage/current-producing device. This is because electric current can only run in closed loops. The total current running in the wire above the ground is distributed in the ground, and fractions of current run everywhere in the subsurface where the resistivity is finite.
+
+## The magnetic dipole
+
+Another example is the magnetic field of the Earth, which to first order is a dipole field. The source is therefore different from what we have seen in the fluid flow and electric potential problems. The Earth's magnetic field (to first order) is the field generated by a magnetic north pole and a magnetic south pole very close together. The dipole vector $\boldsymbol m$ points from the south pole of the dipole to the north pole, and its size is given by the strength of the dipole. The magnetic field $\boldsymbol B$ is given by
+
+$$
+\boldsymbol B = \frac{3\boldsymbol r(\boldsymbol r\cdot\boldsymbol m) - r^2\boldsymbol m}{r^5}.
+$$ (eq:magB)
+
+For this particular solution $\nabla\cdot\boldsymbol B = 0$ for all points in space, also at the source at $\boldsymbol r=\boldsymbol 0$.
+
+## Exercises
+
+1. If $(\boldsymbol v\cdot\hat{\boldsymbol n})\hat{\boldsymbol n}$ in {eq}`eq:vflux` is the fraction of $\boldsymbol v$ that leaves the volume $\mathbb{D}$ through the surface $\mathbb{S}$, what is the fraction of the flow that does not leave the volume $\mathbb{D}$?
+2. Carry out the differentiations to show that the expression for $f(r)$ in {eq}`eq:NF` is correct.
+3. The total current that can be injected into the ground must run in a cable from the source (battery and signal conditioner) to the ground. Once it is in the ground it is free to go anywhere, but the total volume integral must remain equal to the current that runs in the cable, because of the continuity of electric current. We have used the symbol $I$ to denote the total current, and in {eq}`eq:Icur` you have seen a sequence of expressions that resulted in finding the unknown coefficient $A$.
+
+    Another way of finding this result is by observing that the electric potential is the solution of {eq}`eq:laplV` under the condition that a current is injected at the origin. Hence the actual problem is obtained if you take the divergence of both sides of {eq}`eq:Ohm`. This results in $\nabla\cdot\boldsymbol E = \rho\,\nabla\cdot\boldsymbol J$. Integrate both sides of this equation over a spherical volume with fixed radius $r$ and use Gauss' theorem to show that
+
+    $$
+    -\int_{\mathbb{S}}\hat{\boldsymbol n}\cdot(\nabla V)\,\mathrm{d}S = \rho\int_{\mathbb{S}}\hat{\boldsymbol n}\cdot\boldsymbol J\,\mathrm{d}S .
+    $$ (eq:fluxintE)
+
+    The right-hand side is a constant, because it is equal to the total current $I$ that comes from the source and runs in the cable, and therefore it must run out across any spherical surface around the current injection point. Hence, we find
+
+    $$
+    \int_{\mathbb{S}}\hat{\boldsymbol n}\cdot(\nabla V)\,\mathrm{d}S = -\rho I .
+    $$
+
+    Substitute the solution proposed for $V$ of {eq}`eq:pot` with $B=0$ in this equation to verify that $A=\rho I/(4\pi)$.
+4. Evaluate the gradient of the potential expressed in {eq}`eq:Vpdp` and give the expression for the electric current density in the ground at and below the ground surface. Write a Python script that computes the electric potential and the electric current density on the ground surface and in a vertical cross-section, and reproduce the plots of {numref}`fig-dcpoth` and {numref}`fig-dcpotv`. Normalise distance to the electrode spacing $a$ and avoid the points $x=\pm a/2$. You can choose any colour map you like for the potential and choose a contrasting colour for the arrows representing the current lines and directions.
+5. The electric field associated with the electric potential given in {eq}`eq:Vpdp` can be evaluated by taking the gradient of the potential, because of {eq}`eq:EgradV`. Give an argument why the flux integral of the electric field $\int_{\mathbb{S}}\hat{\boldsymbol n}\cdot\boldsymbol E\,\mathrm{d}S = 0$ for every closed and piecewise smooth surface that does not include the current injection and extraction points $x=\pm a/2$.
+6. Verify that the magnetic field expressed in {eq}`eq:magB` is divergence free for all points in space.
+7. Show that the divergence of a vector field in cylindrical and in spherical coordinates is given by
+
+    $$
+    \begin{aligned}
+    \nabla\cdot\boldsymbol v(\varrho,\phi,z) &= \frac{1}{\varrho}\left[\partial_\varrho(\varrho v_\varrho) + \partial_\phi v_\phi\right] + \partial_z v_z, \\
+    \nabla\cdot\boldsymbol v(r,\phi,\theta) &= \frac{1}{r^2}\partial_r(r^2 v_r) + \frac{1}{r\sin(\theta)}\left[\partial_\theta(\sin(\theta)v_\theta) + \partial_\phi v_\phi\right].
+    \end{aligned}
+    $$
+
+    Please remember that in cylindrical coordinates $\varrho=\sqrt{x^2+y^2}$ and in spherical coordinates $r=\sqrt{x^2+y^2+z^2}$!
+8. Consider a general flow field $\boldsymbol v(\boldsymbol r) = \left(v_x(y,z),\,v_y(x,z),\,v_z(x,y)\right)$ flowing in an open space containing a closed surface $\mathbb{S}$. Evaluate the flux integral $\int_{\mathbb{S}}\hat{\boldsymbol n}\cdot\boldsymbol v\,\mathrm{d}S$.
+9. Show that when $\boldsymbol v(\boldsymbol r) = \boldsymbol a\,p(\boldsymbol r)$, where $\boldsymbol a$ is an arbitrary constant vector and $p(\boldsymbol r)$ is a continuously differentiable scalar function, Gauss' integral theorem gives
+
+    $$
+    \int_{\mathbb{S}}p\,\hat{\boldsymbol n}\,\mathrm{d}S = \int_{\mathbb{D}}\nabla p\,\mathrm{d}V,
+    $$
+
+    which is Gauss' theorem for the gradient.
+10. Show that when $\boldsymbol v(\boldsymbol r) = \boldsymbol a\times\boldsymbol w(\boldsymbol r)$, where $\boldsymbol a$ is an arbitrary constant vector and $\boldsymbol w(\boldsymbol r)$ is a continuously differentiable vector function, Gauss' integral theorem gives
+
+    $$
+    \int_{\mathbb{S}}\hat{\boldsymbol n}\times\boldsymbol w\,\mathrm{d}S = \int_{\mathbb{D}}\nabla\times\boldsymbol w\,\mathrm{d}V,
+    $$
+
+    which is Gauss' theorem for the curl.
diff --git a/book/1_gradient_divergence_curl/figures/Cartframe.png b/book/1_gradient_divergence_curl/figures/Cartframe.png
new file mode 100644
index 0000000..b5fdc45
Binary files /dev/null and b/book/1_gradient_divergence_curl/figures/Cartframe.png differ
diff --git a/book/1_gradient_divergence_curl/figures/circulation.png b/book/1_gradient_divergence_curl/figures/circulation.png
new file mode 100644
index 0000000..86cc280
Binary files /dev/null and b/book/1_gradient_divergence_curl/figures/circulation.png differ
diff --git a/book/1_gradient_divergence_curl/figures/crossprod.png b/book/1_gradient_divergence_curl/figures/crossprod.png
new file mode 100644
index 0000000..dcb2e77
Binary files /dev/null and b/book/1_gradient_divergence_curl/figures/crossprod.png differ
diff --git a/book/1_gradient_divergence_curl/figures/cylsphere.png b/book/1_gradient_divergence_curl/figures/cylsphere.png
new file mode 100644
index 0000000..2837e05
Binary files /dev/null and b/book/1_gradient_divergence_curl/figures/cylsphere.png differ
diff --git a/book/1_gradient_divergence_curl/figures/dcpoth.png b/book/1_gradient_divergence_curl/figures/dcpoth.png
new file mode 100644
index 0000000..48f0a9c
Binary files /dev/null and b/book/1_gradient_divergence_curl/figures/dcpoth.png differ
diff --git a/book/1_gradient_divergence_curl/figures/dcpotv.png b/book/1_gradient_divergence_curl/figures/dcpotv.png
new file mode 100644
index 0000000..1a10d17
Binary files /dev/null and b/book/1_gradient_divergence_curl/figures/dcpotv.png differ
diff --git a/book/1_gradient_divergence_curl/figures/rotCartframe.png b/book/1_gradient_divergence_curl/figures/rotCartframe.png
new file mode 100644
index 0000000..3bdd9f5
Binary files /dev/null and b/book/1_gradient_divergence_curl/figures/rotCartframe.png differ
diff --git a/book/1_gradient_divergence_curl/gradient.md b/book/1_gradient_divergence_curl/gradient.md
new file mode 100644
index 0000000..3823a25
--- /dev/null
+++ b/book/1_gradient_divergence_curl/gradient.md
@@ -0,0 +1,159 @@
+# Gradient of a scalar field
+
+Let us consider a scalar field, which is a function of the three spatial coordinates and possibly of time $t$. We write it as $p(x,y,z,t)$ in Cartesian coordinates. A scalar field quantity has a value represented by a single number at every point in space and at each time instant. You can think of the temperature in the room, or the air pressure in our atmosphere, or the density of mass inside the Earth.
+
+Such fields possess iso-surfaces. An **iso-surface** is the collection of points in space where the field has a constant value. The gradient of the field quantity points perpendicular to that iso-surface. Let us investigate how that result is found. The gradient is a partial differential operator given by
+
+$$
+\nabla = \left(\begin{array}{c}
+\dfrac{\partial}{\partial x} \\[2mm]
+\dfrac{\partial}{\partial y} \\[2mm]
+\dfrac{\partial}{\partial z}
+\end{array}\right)
+= \hat{\boldsymbol x}\frac{\partial}{\partial x} + \hat{\boldsymbol y}\frac{\partial}{\partial y} + \hat{\boldsymbol z}\frac{\partial}{\partial z}
+= \hat{\boldsymbol x}\partial_x + \hat{\boldsymbol y}\partial_y + \hat{\boldsymbol z}\partial_z .
+$$
+
+We use the short-hand notation for each scalar partial derivative, e.g. $\partial_x$ for the derivative with respect to $x$. If we apply the gradient to the scalar field quantity $p(x,y,z,t)$ we obtain a vector field quantity that we can write as
+
+$$
+\nabla p(x,y,z,t) = \hat{\boldsymbol x}\,\partial_x p(x,y,z,t) + \hat{\boldsymbol y}\,\partial_y p(x,y,z,t) + \hat{\boldsymbol z}\,\partial_z p(x,y,z,t).
+$$ (eq:gradp)
+
+## The gradient of the distance function
+
+To get an idea about what this implies, let us consider the position vector $\boldsymbol r$ introduced with the Cartesian reference frame. The vector $\boldsymbol r$ is the position vector of the point $(x,y,z)$ in space, which we write as
+
+$$
+\boldsymbol r = x\hat{\boldsymbol x} + y\hat{\boldsymbol y} + z\hat{\boldsymbol z},
+$$
+
+and its length is given by
+
+$$
+r = |\boldsymbol r| = \sqrt{x^2+y^2+z^2}.
+$$
+
+The iso-surface for $r$ has the shape of a spherical surface. We now evaluate each term in the gradient. We begin with the derivative with respect to the coordinate $x$ and obtain
+
+$$
+\begin{aligned}
+\partial_x r &= \frac{1}{2}\frac{1}{\sqrt{x^2+y^2+z^2}}(2x), \\
+             &= \frac{x}{\sqrt{x^2+y^2+z^2}}, \\
+             &= \frac{x}{r}.
+\end{aligned}
+$$
+
+We find similar results for the derivatives with respect to $y$ and $z$, and put them in the vector expression such that we end up with
+
+$$
+\begin{aligned}
+\nabla r &= \hat{\boldsymbol x}\,\partial_x r(x,y,z) + \hat{\boldsymbol y}\,\partial_y r(x,y,z) + \hat{\boldsymbol z}\,\partial_z r(x,y,z), \\
+         &= \frac{x\hat{\boldsymbol x} + y\hat{\boldsymbol y} + z\hat{\boldsymbol z}}{r}, \\
+         &= \frac{\boldsymbol r}{r}.
+\end{aligned}
+$$ (eq:gradr)
+
+From the final expression, we observe that the result is the normalised distance vector. This is what we call the outward unit normal to the spherical surface. It points away from the origin of the reference frame. The physical interpretation is that the gradient of the distance to the origin of the reference frame finds the direction in which the distance increases the most. The vector is placed perpendicular to the iso-surface of the function. We investigate whether this is a general property of the gradient.
+
+Now let us take a different distance, namely relative to an arbitrary other point $\boldsymbol r'$ in space. In that case the displacement vector is given by
+
+$$
+\boldsymbol r-\boldsymbol r' = (x-x')\hat{\boldsymbol x} + (y-y')\hat{\boldsymbol y} + (z-z')\hat{\boldsymbol z},
+$$
+
+and the length of the vector is equal to the distance $d$ given by
+
+$$
+d(x-x',y-y',z-z') = |\boldsymbol r-\boldsymbol r'| = \sqrt{(x-x')^2+(y-y')^2+(z-z')^2}.
+$$
+
+Similar to the previous result, we now find
+
+$$
+\begin{aligned}
+\nabla|\boldsymbol r-\boldsymbol r'| &= \hat{\boldsymbol x}\,\partial_x d + \hat{\boldsymbol y}\,\partial_y d + \hat{\boldsymbol z}\,\partial_z d, \\
+&= \frac{(x-x')\hat{\boldsymbol x} + (y-y')\hat{\boldsymbol y} + (z-z')\hat{\boldsymbol z}}{|\boldsymbol r-\boldsymbol r'|},
+\end{aligned}
+$$
+
+which we write in vector form as
+
+$$
+\nabla|\boldsymbol r-\boldsymbol r'| = \frac{\boldsymbol r-\boldsymbol r'}{|\boldsymbol r-\boldsymbol r'|}.
+$$ (eq:gradd)
+
+We find again an outward unit vector, pointing away from the point at $\boldsymbol r'$ and perpendicular to the spherical iso-surface.
+
+## The total derivative and the direction of steepest increase
+
+Now let us consider a field quantity $p(x,y,z,t)$ and assume we analyse this function for a single moment in time. The gradient of this function is expressed in {eq}`eq:gradp`. Let $\mathrm{d}p$ be the change in $p$ from a point $\boldsymbol r$ to another point $\boldsymbol r'$, which means that $p(\boldsymbol r)=p$ and $p(\boldsymbol r')=p+\mathrm{d}p$. We do not specify where $\boldsymbol r$ and $\boldsymbol r'$ are located, and they can be on different iso-surfaces or on the same one.
+
+The change in $p$ due to a displacement in the $x$-direction from $\boldsymbol r$ to $\boldsymbol r'$ is given by $(\partial p/\partial x)\mathrm{d}x$, while keeping $y$ and $z$ constant. Similarly, the changes in $p$ due to displacements in the $y$- and $z$-directions are given by $(\partial p/\partial y)\mathrm{d}y$ and $(\partial p/\partial z)\mathrm{d}z$. Hence, the change in $p$ along the vector from $\boldsymbol r$ to $\boldsymbol r'$ is
+
+$$
+\mathrm{d}p = \partial_x p\,\mathrm{d}x + \partial_y p\,\mathrm{d}y + \partial_z p\,\mathrm{d}z.
+$$
+
+We can write this expression as a scalar product of two vectors,
+
+$$
+\mathrm{d}p = \left(\hat{\boldsymbol x}\partial_x p + \hat{\boldsymbol y}\partial_y p + \hat{\boldsymbol z}\partial_z p\right)\cdot\left(\hat{\boldsymbol x}\,\mathrm{d}x + \hat{\boldsymbol y}\,\mathrm{d}y + \hat{\boldsymbol z}\,\mathrm{d}z\right),
+$$
+
+where the symbol $\cdot$ is used to denote scalar multiplication of two vectors. We recognise this expression as
+
+$$
+\mathrm{d}p = (\nabla p)\cdot\mathrm{d}\boldsymbol r.
+$$
+
+Suppose there is an angle $\psi$ between the two vectors, then
+
+$$
+\mathrm{d}p = |\nabla p|\,|\mathrm{d}\boldsymbol r|\cos(\psi) = |\nabla p|\,\mathrm{d}r\cos(\psi),
+$$
+
+and we can write the total derivative with respect to $r$ as
+
+$$
+\frac{\mathrm{d}p}{\mathrm{d}r} = |\nabla p|\cos(\psi).
+$$ (eq:totder)
+
+The left-hand side of {eq}`eq:totder` means the rate of change along the path from $\boldsymbol r$ to $\boldsymbol r'$, as indicated by $r$. The right-hand side shows that this can never be larger than the magnitude of the gradient of $p$. We conclude that the rate of change is equal to the magnitude of the gradient only if the point $\boldsymbol r'$ is located along the resulting vector after taking the gradient of $p$, because then $\psi=0$. For all other locations the rate of change is less. We have seen this when we evaluated the gradient of the distance function $r$.
+
+Now we have found that the magnitude of the gradient of a scalar field quantity is equal to the maximum rate of change of that field quantity with respect to position. What is left is direction, which is relatively easy to understand. It is clear that when the point $\boldsymbol r'$ is on the same iso-surface as the point $\boldsymbol r$, the rate of change is zero. We investigate differentials, which means we look at points $\boldsymbol r'$ that approach the point $\boldsymbol r$. The rate of change of the function $p$ is maximal when the point $\boldsymbol r'$ moves away in the direction perpendicular to the iso-surface. If there is a part of the path from $\boldsymbol r$ to $\boldsymbol r'$ that has a component along the iso-surface, there is no change along that part of the path, which would reduce the rate of change. Hence, the gradient results in a vector that points along the unit normal of the iso-surface in the direction where the rate is positive. We can express this as
+
+$$
+\nabla p(x,y,z,t) = |\nabla p(x,y,z,t)|\,\hat{\boldsymbol n},
+$$
+
+where $\hat{\boldsymbol n}$ is the unit normal vector on the iso-surface pointing to the positive rate of change.
+
+:::{admonition} The gradient in words
+:class: tip
+The gradient of any scalar field quantity $p(x,y,z,t)$ finds the direction in which the scalar quantity increases the most, for a fixed moment in time, and its magnitude is that maximum rate of increase.
+:::
+
+## The gradient in curvilinear coordinates
+
+Earlier we introduced the spherical and cylindrical coordinate systems. We had to look at the rotation matrices between the coordinate systems to be able to move back and forth. It is now time to generalise our notation of the gradient such that it can be used in spherical and cylindrical coordinate systems as well. Let us write
+
+$$
+\nabla p(x,y,z,t) = \sum_{i=1}^{3}\frac{1}{c_i}\frac{\partial p}{\partial x_i}\hat{\boldsymbol e}_i .
+$$ (eq:gradgen)
+
+In this expression the coefficients $c_i$ are scale factors, the coordinates are $(x_1,x_2,x_3)$, and both depend on the coordinate system we want to do our analysis in. The base vectors $(\hat{\boldsymbol e}_1,\hat{\boldsymbol e}_2,\hat{\boldsymbol e}_3)$ were introduced with the coordinate systems, and we use them here again, but we make them depend on the coordinate system as well.
+
+In the Cartesian frame these would be $(x_1,x_2,x_3)=(x,y,z)$, $(c_1,c_2,c_3)=(1,1,1)$ and $(\hat{\boldsymbol e}_1,\hat{\boldsymbol e}_2,\hat{\boldsymbol e}_3)=(\hat{\boldsymbol x},\hat{\boldsymbol y},\hat{\boldsymbol z})$. Now, in spherical coordinates, we already know all the ingredients, so we can fill them in. We find $(x_1,x_2,x_3)=(r,\theta,\phi)$, $(c_1,c_2,c_3)=(1,r,r\sin(\theta))$ and $(\hat{\boldsymbol e}_1,\hat{\boldsymbol e}_2,\hat{\boldsymbol e}_3)=(\hat{\boldsymbol r},\hat{\boldsymbol\theta},\hat{\boldsymbol\phi})$. Substituting these in {eq}`eq:gradgen` results in
+
+$$
+\nabla p(r,\theta,\phi) = \frac{\partial p}{\partial r}\hat{\boldsymbol r} + \frac{1}{r}\frac{\partial p}{\partial\theta}\hat{\boldsymbol\theta} + \frac{1}{r\sin(\theta)}\frac{\partial p}{\partial\phi}\hat{\boldsymbol\phi}.
+$$ (eq:gradsph)
+
+## Exercises
+
+1. Evaluate $\nabla|\boldsymbol r|^{-1}$, $\nabla|\boldsymbol r|^{-n}$, $\nabla|\boldsymbol r-\boldsymbol r'|^{-1}$, and $\nabla\left(|\boldsymbol r-\boldsymbol a|^{-1} - |\boldsymbol r+\boldsymbol a|^{-1}\right)$, with $\boldsymbol a=(1,0,0)$. For each result, plot several iso-surfaces and plot the vectors that correspond to the gradients.
+2. Instead of taking the gradient with respect to the point $\boldsymbol r$, we can take it with respect to $\boldsymbol r'$, which we write as $\nabla' = \hat{\boldsymbol x}\partial_{x'} + \hat{\boldsymbol y}\partial_{y'} + \hat{\boldsymbol z}\partial_{z'}$. Evaluate $\nabla'|\boldsymbol r-\boldsymbol r'|^{-1}$ and express it in terms of $\nabla|\boldsymbol r-\boldsymbol r'|^{-1}$.
+3. We have seen that the outward unit normal of a spherical surface around the point $\boldsymbol r'$ is given by $\hat{\boldsymbol n} = \nabla|\boldsymbol r-\boldsymbol r'| = (\boldsymbol r-\boldsymbol r')/|\boldsymbol r-\boldsymbol r'|$. Use your understanding of the property of the gradient to explain why $\hat{\boldsymbol n}\cdot\nabla|\boldsymbol r-\boldsymbol r'|^{-1}<0$ for $\boldsymbol r\ne\boldsymbol r'$.
+4. What is the relation between $A_r,A_\phi,A_\theta$ introduced in the exercises on coordinate systems and the scale factors $c_i$ here?
+5. Determine the expression for the gradient in cylindrical coordinates. You can use the coordinates $(\varrho,\phi,z)$ to avoid confusion with the three-dimensional radius $r$ that is used in spherical coordinates and represents the distance between two points in 3D space in general. In cylindrical coordinates $r=\sqrt{\varrho^2+z^2}$.
diff --git a/book/1_gradient_divergence_curl/intro.md b/book/1_gradient_divergence_curl/intro.md
index 491c78c..c653859 100644
--- a/book/1_gradient_divergence_curl/intro.md
+++ b/book/1_gradient_divergence_curl/intro.md
@@ -1 +1,95 @@
-## Introduction
\ No newline at end of file
+# Introduction
+
+Mathematics is the language to describe physics, and physics gives the empirical content consisting of experiments. Mathematics is grammar, the model is a simplification that is meant to understand a measurement. The model states *if $a$ then $b$*, while the experiment states *if $A$ then $B$*. If $b$ matches $B$ to our satisfaction, we adopt the model; if not, we reject it. But even when adopted, we must continue to scrutinise the model, and later an experiment will come that makes us understand that the model was adopted earlier but needs revision.
+
+Scientific progress is made through doubt.
+
+AI changes the direction of knowledge: not any more from rule to reality, but from reality to rule. The data are used to understand and to generate a model.
+
+Physical objects can be described in a quantitative way only with the aid of mathematics. In classical physics, all physical objects are geometric objects. The course *Fields and Waves* combines the physics of fields and waves with the mathematical tools required to describe these phenomena. The course enables students to develop essential understanding of the physical interpretation of mathematical formulations. Examples are radar waves that are used for earth surface and subsurface observations from antennas placed on satellites, airplanes, and close to or on the ground surface; sound waves, electromagnetic diffusion fields, electric and magnetic potential fields, and the gravity field to probe the earth's interior, the surface and the atmosphere. They are used to characterise layers and objects in terms of physical and geological parameters, and to monitor dynamic processes as well. The sources for these fields can be natural or anthropic.
+
+## Dimensions and units
+
+Seven fundamental quantities, or dimensions, have been defined. They are *length* ($L$), *time* ($T$), *temperature* ($\mathcal{T}$), *mass* ($M$), *electric current* ($I$), *amount of substance* ($n$) and *luminous intensity* ($\mathcal{L}$). Other dimensions are then secondary and can be written in terms of several or all of these seven. Electric charge has fundamental dimensions $IT$ and the fundamental dimensions of electric field are given by $ML/(IT^3)$.
+
+Dimensions must be given a value to work with them numerically. For this the international agreement is to use the so-called metric system, and the present-day variant has the seven fundamental units that correspond to the fundamental dimensions. This system is known as the SI system, from the French *Système Internationale d'Unités*, known as the International System of Units. In this system, the dimension length has unit *meter* (m), the dimension time has unit *second* (s), the dimension temperature has unit *kelvin* (K), the dimension mass has unit *kilogram* (kg), the dimension electric current has unit *ampere* (A), the dimension amount of substance has unit *mole* (mol), and the dimension luminous intensity has unit *candela* (cd). These units are defined as follows:
+
+- **Meter** (m). One meter is equal to the path length travelled by light in vacuum in a time of $t = 1/299\,792\,458$ second. This defines the electromagnetic wave propagation velocity in vacuum as $c_0 = 299\,792\,458$ m/s.
+- **Second** (s). One second is equal to the duration of $9\,192\,631\,770$ periods of radiation corresponding to the transition between two hyperfine levels of the ground state of cesium 133. This is now known as the atomic clock. Atomic clocks are accurate to approximately 1 microsecond per year. Before the atomic clock the second was defined as the mean solar day divided by 86400, but because the earth's rotation around the sun is slowing down it was regarded inaccurate as a standard. The two standards differ in the order of 1 second per year. Distant fast rotating pulsars are, with their 1000 revolutions per second, a possible new replacement for the atomic clocks and will then yield a standard with an accuracy in the order of nanoseconds per year.
+- **Kelvin** (K). One kelvin is the temperature equal to $1/273.16$ of the triple point of water, defining the triple point of water as $273.16$ kelvin. Water boils at a temperature of $T = 100\,^{\circ}\mathrm{C} = 373.15$ K.
+- **Kilogram** (kg). For a long time this was the only unit still defined by a physical prototype: a cylinder of platinum and iridium alloy stored in Sèvres, France. In this sense the kilogram was an anomaly among the unit definitions.
+- **Ampere** (A). One ampere is equal to the electric current flowing in each of two infinitely long parallel wires in vacuum separated by one meter, which produces a force of 200 nanonewton per meter of length.
+- **Mole** (mol). The mole is defined as the amount of substance of a system which contains as many "elemental entities" (e.g., atoms, molecules, ions, electrons) as there are atoms in $0.012$ kg of carbon-12. It is related to the number of particles, which is the Avogadro constant $N_A = 6.022\,141\,79\times 10^{23}$ mol$^{-1}$.
+- **Candela** (cd). One candela is the luminous intensity equal to that of $1/600\,000$ square meter of a perfect radiator at the temperature of freezing platinum at a pressure of one standard atmosphere.
+
+:::{note}
+The definitions above are the ones you will meet in most textbooks, and they are the ones used throughout these notes. Since the SI revision of 2019, the base units are instead fixed by assigning exact values to seven defining constants: $\Delta\nu_{\mathrm{Cs}}$, $c_0$, $h$, $e$, $k_{\mathrm{B}}$, $N_A$ and $K_{\mathrm{cd}}$. The kilogram is now realised from the Planck constant $h$ rather than from the prototype cylinder, and the ampere from the elementary charge $e$ rather than from the force between two wires. The numerical values change by far less than any measurement we make in this course, so nothing in what follows depends on which convention you have in mind.
+:::
+
+The other units are called secondary, or derived, units and can all be expressed as combinations of these seven. The International System of Units also recommends the use of abbreviations of units in steps of three orders of magnitude. To take length as an example, it is recommended to use 10 mm over 1 cm. In these notes the SI system is used and the abbreviation recommendation is adhered to. The metric system and its scientific prefixes are given in {numref}`tab-si-prefixes`.
+
+```{list-table} Numbers in the metric system and their prefixes.
+:header-rows: 1
+:name: tab-si-prefixes
+
+* - Numerical value
+  -
+  - Prefix
+  - Symbol
+* - 1 000 000 000 000 000 000
+  - $10^{18}$
+  - exa
+  - E
+* - 1 000 000 000 000 000
+  - $10^{15}$
+  - peta
+  - P
+* - 1 000 000 000 000
+  - $10^{12}$
+  - tera
+  - T
+* - 1 000 000 000
+  - $10^{9}$
+  - giga
+  - G
+* - 1 000 000
+  - $10^{6}$
+  - mega
+  - M
+* - 1 000
+  - $10^{3}$
+  - kilo
+  - k
+* - 1
+  - $1$
+  - one
+  - –
+* - 0.001
+  - $10^{-3}$
+  - milli
+  - m
+* - 0.000 001
+  - $10^{-6}$
+  - micro
+  - $\mu$
+* - 0.000 000 001
+  - $10^{-9}$
+  - nano
+  - n
+* - 0.000 000 000 001
+  - $10^{-12}$
+  - pico
+  - p
+* - 0.000 000 000 000 001
+  - $10^{-15}$
+  - femto
+  - f
+* - 0.000 000 000 000 000 001
+  - $10^{-18}$
+  - atto
+  - a
+```
+
+## What follows
+
+The remaining pages of this introduction deal with sums, series and approximations, and with reference frames, symbols and notations for scalar, vector, and matrix quantities, together with the notion of time and temporal variations of a function. The pages after those describe the three main spatial derivative operators — gradient, divergence, and curl — and discuss their physical meaning. Later chapters describe potential, diffusive, and wave fields, respectively.
diff --git a/book/1_gradient_divergence_curl/labs/fwtools.py b/book/1_gradient_divergence_curl/labs/fwtools.py
new file mode 100644
index 0000000..3c0bcd0
--- /dev/null
+++ b/book/1_gradient_divergence_curl/labs/fwtools.py
@@ -0,0 +1,594 @@
+"""Plotting and self-check helpers for the ECTB2140 *Fields and Waves* labs.
+
+You are **not** expected to read or edit this file during a lab session. It
+exists so that your time goes on the physics -- writing fields, gradients and
+divergences -- rather than on plotting boilerplate.
+
+Everything here works on a 3-D Cartesian grid built with ``indexing='ij'``,
+which means the array axes are (x, y, z) in that order. That choice matters:
+it makes ``np.gradient(F, dx, dy, dz)`` return the derivatives in the same
+order as the coordinates, with no index gymnastics.
+
+Requires: numpy, matplotlib, plotly.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import matplotlib.pyplot as plt
+import plotly.graph_objects as go
+
+__all__ = [
+    "z0_index", "slice_z0",
+    "box_indices", "line_integral", "area_integral", "volume_integral",
+    "show_isosurfaces", "show_cones", "show_scalar_slice", "show_field_slice",
+    "check", "check_shape", "check_close", "check_abs", "check_scalar",
+]
+
+# --------------------------------------------------------------------------
+# Grid helpers (the grid itself is built in the open, on the lab page)
+# --------------------------------------------------------------------------
+
+def z0_index(Z: np.ndarray) -> int:
+    """Index k of the z = 0 plane in an indexing='ij' grid."""
+    return int(np.argmin(np.abs(Z[0, 0, :])))
+
+
+def slice_z0(F: np.ndarray, Z: np.ndarray) -> np.ndarray:
+    """The z = 0 plane of a scalar field, as a 2-D (x, y) array."""
+    return F[:, :, z0_index(Z)]
+
+
+# --------------------------------------------------------------------------
+# Integration over grid-aligned boxes, faces and edges
+#
+# These evaluate the integrals in the definitions of the divergence and the
+# curl, and in the divergence and Stokes theorems. They are quadrature
+# boilerplate: the trapezoidal weights below simply stop the end samples
+# from being counted as full cells.
+# --------------------------------------------------------------------------
+
+def _trapezoid_weights(n: int) -> np.ndarray:
+    w = np.ones(n)
+    w[0] = w[-1] = 0.5
+    return w
+
+
+def box_indices(X: np.ndarray, half_width: float):
+    """Index range (i0, i1) of the sub-cube |x|, |y|, |z| <= ``half_width``.
+
+    The same pair works on all three axes because the grid is cubic. The
+    returned indices are snapped to the nearest grid planes, so ask for a
+    half-width that is a multiple of the spacing (0.6, 1.0 and 1.4 m are
+    exact on the default 61-point grid) if you want the box you asked for.
+    """
+    axis = X[:, 0, 0]
+    i0 = int(np.argmin(np.abs(axis + half_width)))
+    i1 = int(np.argmin(np.abs(axis - half_width)))
+    return i0, i1
+
+
+def line_integral(F1: np.ndarray, ds: float) -> float:
+    """Integrate a 1-D array of samples along the line it spans.
+
+    One edge of a closed loop, for the circulation integral, the way
+    ``area_integral`` handles one face of a box for the flux. Same
+    trapezoidal rule, same refusal to integrate over masked samples.
+    """
+    F1 = np.asarray(F1, float)
+    _reject_masked(F1, "This edge passes through masked samples")
+    return float(np.sum(F1 * _trapezoid_weights(F1.shape[0])) * ds)
+
+
+def area_integral(F2: np.ndarray, da: float, db: float) -> float:
+    """Integrate a 2-D array of samples over the rectangle it spans.
+
+    Use it on one face of a box to evaluate that face's contribution to a
+    surface integral. The rule is the trapezoidal one, second order in the
+    spacing: on the fields in this lab the flux error falls from 0.17% at
+    n = 21 to 0.018% at n = 61.
+
+    Raises if any sample is NaN. A masked sample silently integrated as zero
+    returns a plausible and wrong number -- which is what happens if you put
+    a face inside a region you have masked out.
+    """
+    F2 = np.asarray(F2, float)
+    _reject_masked(F2, "This face passes through masked samples")
+    wa, wb = _trapezoid_weights(F2.shape[0]), _trapezoid_weights(F2.shape[1])
+    return float(np.sum(F2 * wa[:, None] * wb[None, :]) * da * db)
+
+
+def volume_integral(F3: np.ndarray, dx: float, dy: float, dz: float) -> float:
+    """Integrate a 3-D array of samples over the box it spans.
+
+    Trapezoidal, second order, and it raises on masked samples for the same
+    reason ``area_integral`` does.
+    """
+    F3 = np.asarray(F3, float)
+    _reject_masked(F3, "This box contains masked samples")
+    wx, wy, wz = (_trapezoid_weights(m) for m in F3.shape)
+    w = wx[:, None, None] * wy[None, :, None] * wz[None, None, :]
+    return float(np.sum(F3 * w) * dx * dy * dz)
+
+
+def _reject_masked(F, what):
+    n = int(np.count_nonzero(~np.isfinite(F)))
+    if n:
+        raise ValueError(
+            f"{what} ({n} of {F.size} are NaN or infinite). Integrating them "
+            f"as zero would return a plausible but wrong number. Move the "
+            f"surface outside the masked region, or unmask the field.")
+
+
+# --------------------------------------------------------------------------
+# 3-D views (plotly)
+# --------------------------------------------------------------------------
+
+# Directional shading. Without it plotly lights an isosurface almost flatly
+# and a nest of transparent spheres reads as a set of flat rings; the
+# specular highlight and limb darkening are what make it look like a ball.
+_LIGHTING = dict(ambient=0.35, diffuse=0.9, specular=0.5, roughness=0.4, fresnel=0.2)
+_LIGHTPOSITION = dict(x=100, y=200, z=200)
+
+
+def show_isosurfaces(X, Y, Z, F, levels, *, title="", label="", opacity=0.3,
+                     colorscale="Viridis", reversescale=False, show_caps=False,
+                     size=620, step=2, opacity_slider=True, slice_z=None):
+    """Draw one or more isosurfaces (level sets) of a scalar field F.
+
+    An isosurface is the set of points where F takes one fixed value -- the
+    3-D analogue of a contour line. Pass an **evenly spaced** list of levels
+    and look at the nesting; anything else is refused, because plotly draws
+    evenly spaced surfaces between the extremes and would quietly move them.
+
+    Parameters that matter for seeing the shape
+    -------------------------------------------
+    opacity_slider : bool
+        Adds a slider under the figure. Drag it up towards 1 and the outermost
+        surface becomes a solid, shaded ball; drag it down towards 0.1 and it
+        turns to glass so the inner surfaces show through. Sweeping it is the
+        quickest way to convince yourself these are shells and not discs.
+    label : str
+        Colorbar title. Give it the physical quantity and its unit.
+        **Plotly does not render LaTeX here.** Colorbar titles accept plain
+        text plus a small HTML subset (``<sub>``, ``<sup>``, ``<b>``), so
+        write ``"|\u2207r|  [-]"`` and ``"[m<sup>-2</sup>]"`` with Unicode
+        symbols -- a ``$...$`` label silently comes out as garbled glyphs.
+        The matplotlib helpers below are the opposite: mathtext works there.
+    reversescale : bool
+        Flip the colorscale. Needed for a signed field on ``"RdBu"``: plotly
+        runs that scale dark *red* at the low end and dark *blue* at the high
+        end, the opposite of matplotlib's ``"RdBu_r"`` used by the 2-D
+        helpers below. Without this flag a positive lobe drawn in 3-D comes
+        out blue while the same lobe in the slice beside it is red.
+    slice_z : float or None
+        If given, also draw a filled cut plane at that value of z, exposing
+        the interior. A strong depth cue, at the cost of hiding part of the
+        nesting.
+    """
+    levels = np.sort(np.atleast_1d(np.asarray(levels, dtype=float)))
+    # plotly draws surface_count EVENLY SPACED surfaces between isomin and
+    # isomax; it never sees the individual values. Unevenly spaced levels
+    # would therefore be silently redrawn at the wrong values, so refuse them
+    # rather than return a picture that lies.
+    if levels.size > 2:
+        gaps = np.diff(levels)
+        if not np.allclose(gaps, gaps[0], rtol=1e-6):
+            drawn = np.linspace(levels[0], levels[-1], levels.size)
+            raise ValueError(
+                f"levels must be evenly spaced: plotly would draw "
+                f"{np.round(drawn, 4).tolist()} instead of "
+                f"{np.round(levels, 4).tolist()}. Use an evenly spaced set, "
+                f"or call this once per level.")
+    # Subsample before handing the volume to plotly. A full 61^3 grid embeds
+    # ~12 MB of JSON per figure; every other point looks identical on screen.
+    sl = (slice(None, None, step),) * 3
+    X, Y, Z, F = X[sl], Y[sl], Z[sl], np.asarray(F)[sl]
+    trace = go.Isosurface(
+        x=X.ravel(), y=Y.ravel(), z=Z.ravel(), value=np.asarray(F).ravel(),
+        isomin=float(levels.min()), isomax=float(levels.max()),
+        surface_count=int(levels.size), opacity=opacity,
+        colorscale=colorscale, reversescale=reversescale, showscale=True,
+        colorbar=dict(title=label, len=0.7),
+        lighting=_LIGHTING, lightposition=_LIGHTPOSITION,
+        caps=dict(x_show=show_caps, y_show=show_caps, z_show=show_caps),
+    )
+    if slice_z is not None:
+        trace.slices = dict(z=dict(show=True, locations=[float(slice_z)]))
+    fig = go.Figure(trace)
+    _style_3d(fig, title, size, bottom_margin=55 if opacity_slider else 0)
+    if opacity_slider:
+        _add_opacity_slider(fig, opacity)
+    return _display(fig)
+
+
+def _add_opacity_slider(fig, current):
+    """A client-side opacity control: no kernel needed once the figure exists."""
+    values = [round(0.1 * i, 1) for i in range(1, 11)]
+    active = int(np.argmin([abs(v - current) for v in values]))
+    fig.update_layout(sliders=[dict(
+        active=active,
+        currentvalue=dict(prefix="opacity: ", font=dict(size=13)),
+        pad=dict(t=8, b=8), len=0.7, x=0.15, y=0,
+        steps=[dict(method="restyle", args=[{"opacity": v}], label=f"{v:.1f}")
+               for v in values],
+    )])
+
+
+def show_cones(X, Y, Z, Ax, Ay, Az, *, step=8, title="", label="", unit="",
+               size=620,
+               normalise=False, length=None, head=0.35, colorscale="Viridis",
+               slider=True, width=4, log_colour=None):
+    """Draw a 3-D vector field as arrows: a shaft with a barbed head.
+
+    Every arrow is built from line segments -- a shaft, plus four barbs swept
+    back from the tip. plotly's ``go.Cone`` is not used: a cone takes both its
+    size and its colour from the norm of the vector it is given, so size and
+    colour cannot be set independently, and a field whose magnitudes are all
+    close to 1 comes out with heads larger than the box.
+
+    Only every ``step``-th sample in each direction is drawn; an arrow at every
+    grid point is an unreadable haystack.
+
+    Parameters
+    ----------
+    label, unit : str
+        The plotted quantity and its unit, kept apart, e.g.
+        ``label="|<b>E</b>|", unit="V/m"``. A linear colorbar is then titled
+        ``|<b>E</b>|  [V/m]``; a log one ``log<sub>10</sub>(|<b>E</b>| /
+        (V/m))``, because the logarithm of a dimensional quantity does not
+        carry that quantity's unit. Plotly renders no LaTeX here -- see the
+        note in ``show_isosurfaces``.
+    normalise : bool
+        Draw every arrow the same length, showing direction only. Use it for
+        fields whose magnitude spans orders of magnitude, where true-to-scale
+        arrows leave a few giants and a lot of invisible dust. The magnitude is
+        not lost -- it is still in the colour.
+    length : float or None
+        Length of the longest arrow, in metres. Defaults to 0.85 of the
+        spacing between drawn arrows, so a full-length arrow almost touches
+        its neighbour.
+    head : float
+        Fraction of an arrow taken up by its head.
+    slider : bool
+        Add a size slider under the figure, scaling whole arrows (head
+        included) between 0.5x and 2x.
+    log_colour : bool or None
+        Colour by log10|A| rather than |A|. ``None`` decides automatically and
+        switches over once the magnitude spans more than a factor of 50: on a
+        linear scale a $1/r^2$ field puts all but a handful of arrows into the
+        bottom percent of the colour range, where they are indistinguishable.
+    """
+    sl = (slice(None, None, step),) * 3
+    x, y, z = X[sl].ravel(), Y[sl].ravel(), Z[sl].ravel()
+    u = np.asarray(Ax)[sl].ravel()
+    v = np.asarray(Ay)[sl].ravel()
+    w = np.asarray(Az)[sl].ravel()
+
+    finite = np.isfinite(u) & np.isfinite(v) & np.isfinite(w)
+    x, y, z, u, v, w = (a[finite] for a in (x, y, z, u, v, w))
+
+    mag = np.sqrt(u**2 + v**2 + w**2)
+    safe = np.maximum(mag, 1e-30)
+    ux, uy, uz = u / safe, v / safe, w / safe          # unit direction
+
+    spacing = float(abs(X[step, 0, 0] - X[0, 0, 0])) if X.shape[0] > step else 1.0
+    base = 0.85 * spacing if length is None else float(length)
+    rel = np.ones_like(mag) if normalise else mag / max(float(mag.max()), 1e-30)
+
+    positive = mag[mag > 0]
+    if log_colour is None:
+        log_colour = (positive.size > 0
+                      and float(positive.max()) > 50.0 * float(positive.min()))
+    name = label or "|A|"
+    if log_colour:
+        cval = np.log10(np.maximum(mag, float(positive.min())))
+        # log10 of a dimensional quantity is dimensionless: the unit belongs
+        # inside the logarithm, as a divisor, never appended in brackets.
+        clabel = (f"log<sub>10</sub>({name} / {unit})" if unit
+                  else f"log<sub>10</sub> {name}")
+    else:
+        cval = mag
+        clabel = f"{name}  [{unit}]" if unit else name
+
+    px, py, pz = _arrow_lines(x, y, z, ux, uy, uz, rel * base, head)
+    fig = go.Figure(go.Scatter3d(
+        x=px, y=py, z=pz, mode="lines", hoverinfo="skip", showlegend=False,
+        line=dict(color=np.tile(np.repeat(cval, 3), _SEGMENTS_PER_ARROW),
+                  colorscale=colorscale, width=width,
+                  cmin=float(cval.min()), cmax=float(cval.max()),
+                  showscale=True, colorbar=dict(title=clabel, len=0.7)),
+    ))
+    _style_3d(fig, title, size, bottom_margin=75 if slider else 0)
+    # Pin the box to the sampled volume; without this the arrows themselves
+    # drive the autorange and the domain silently grows.
+    fig.update_scenes(
+        xaxis=dict(range=[float(X.min()), float(X.max())], title="x [m]"),
+        yaxis=dict(range=[float(Y.min()), float(Y.max())], title="y [m]"),
+        zaxis=dict(range=[float(Z.min()), float(Z.max())], title="z [m]"),
+    )
+    if slider:
+        _add_arrow_slider(fig, x, y, z, ux, uy, uz, rel, base, head)
+    return _display(fig)
+
+
+_SEGMENTS_PER_ARROW = 5          # one shaft, four barbs
+
+
+def _arrow_lines(x, y, z, ux, uy, uz, lengths, head):
+    """One polyline per segment, all arrows in one flat pair of arrays.
+
+    Segments are separated by NaN, which plotly renders as a break. The four
+    barbs are swept back from the tip in two mutually perpendicular planes, so
+    the head reads as a head from any viewing angle.
+    """
+    n = x.size
+    tx, ty, tz = x + ux * lengths, y + uy * lengths, z + uz * lengths
+
+    d = np.stack([ux, uy, uz], axis=1)
+    # A reference direction not parallel to d, so the cross product is stable.
+    ref = np.where(np.abs(uz)[:, None] < 0.9,
+                   np.array([0.0, 0.0, 1.0]), np.array([1.0, 0.0, 0.0]))
+    p = np.cross(d, ref)
+    p /= np.maximum(np.linalg.norm(p, axis=1, keepdims=True), 1e-30)
+    q = np.cross(d, p)
+
+    barb = lengths * head
+    spread = 0.45
+    xs, ys, zs = [], [], []
+
+    def add(x0, y0, z0, x1, y1, z1):
+        for a, b, out in ((x0, x1, xs), (y0, y1, ys), (z0, z1, zs)):
+            seg = np.empty(3 * n)
+            seg[0::3], seg[1::3], seg[2::3] = a, b, np.nan
+            out.append(seg)
+
+    add(x, y, z, tx, ty, tz)                                  # the shaft
+    for side in (p, -p, q, -q):                               # the four barbs
+        add(tx, ty, tz,
+            tx - ux * barb + side[:, 0] * barb * spread,
+            ty - uy * barb + side[:, 1] * barb * spread,
+            tz - uz * barb + side[:, 2] * barb * spread)
+    return (np.concatenate(xs).astype(np.float32),
+            np.concatenate(ys).astype(np.float32),
+            np.concatenate(zs).astype(np.float32))
+
+
+def _add_arrow_slider(fig, x, y, z, ux, uy, uz, rel, base, head):
+    """One client-side control scaling whole arrows, head included."""
+    scales = [0.5, 0.75, 1.0, 1.5, 2.0]
+    steps = []
+    for sc in scales:
+        px, py, pz = _arrow_lines(x, y, z, ux, uy, uz, rel * base * sc, head)
+        steps.append(dict(method="restyle", label=f"{sc:g}x",
+                          args=[{"x": [px], "y": [py], "z": [pz]}, [0]]))
+    fig.update_layout(sliders=[dict(
+        active=scales.index(1.0), steps=steps, len=0.7, x=0.15, y=0,
+        pad=dict(t=8, b=8),
+        currentvalue=dict(prefix="arrow size: ", font=dict(size=13)),
+    )])
+
+
+def _display(fig):
+    """Show the figure and return nothing.
+
+    Jupyter renders only the value of a cell's LAST expression, so a plotting
+    call followed by a self-check would otherwise draw nothing at all. Showing
+    it here makes the call work wherever it sits; returning None keeps it from
+    being drawn a second time when it does happen to come last.
+    """
+    fig.show()
+    return None
+
+
+def _style_3d(fig, title, size, bottom_margin=0):
+    fig.update_layout(
+        title=title, width=size, height=size,
+        margin=dict(l=0, r=0, t=40 if title else 0, b=bottom_margin),
+        scene=dict(
+            xaxis_title="x [m]", yaxis_title="y [m]", zaxis_title="z [m]",
+            aspectmode="cube",          # equal aspect: never distort a field
+            camera=dict(eye=dict(x=1.6, y=1.6, z=1.1)),
+        ),
+    )
+
+
+# --------------------------------------------------------------------------
+# 2-D views of a coordinate plane (matplotlib)
+# --------------------------------------------------------------------------
+
+def _plane_slice(X, Y, Z, F, plane):
+    """Cut a 3-D field on a coordinate plane through the origin.
+
+    ``plane="z"`` gives the z = 0 plane in (x, y); ``plane="y"`` gives the
+    y = 0 plane in (x, z) -- the vertical cross-section a geophysical survey
+    is usually drawn on. Returns the two 1-D axes, the 2-D field, and the two
+    axis labels.
+    """
+    F = None if F is None else np.asarray(F)
+    if plane == "z":
+        k = z0_index(Z)
+        return (X[:, 0, k], Y[0, :, k], None if F is None else F[:, :, k],
+                "$x$ [m]", "$y$ [m]")
+    if plane == "y":
+        j = int(np.argmin(np.abs(Y[0, :, 0])))
+        return (X[:, j, 0], Z[0, j, :], None if F is None else F[:, j, :],
+                "$x$ [m]", "$z$ [m]")
+    raise ValueError(f"plane must be 'z' or 'y', not {plane!r}")
+
+def show_scalar_slice(X, Y, Z, F, *, title="", label="", cmap=None,
+                      levels=25, symmetric=False, percentile=99, ax=None,
+                      colorbar=True, vmin=None, vmax=None, plane="z"):
+    """Filled contours of a scalar field on a coordinate plane.
+
+    ``colorbar`` is drawn whether or not the axes was supplied by the caller;
+    a panel in a side-by-side comparison needs its scale just as much as a
+    standalone figure does. ``show_field_slice`` passes ``colorbar=False``
+    because it adds its own.
+
+    Pass ``vmin``/``vmax`` to pin the colour limits. Without them the limits
+    come from percentiles of *this* panel, so two panels of a comparison end
+    up on different scales and the extremes are clipped -- give both panels
+    the same explicit pair whenever the point is that they match.
+
+    ``cmap`` defaults to a diverging map when ``symmetric=True`` and a
+    sequential one otherwise, so a one-signed field never gets a colour scale
+    implying a meaningful zero crossing.
+
+    ``plane="z"`` cuts z = 0, ``plane="y"`` cuts y = 0 for a vertical section.
+    """
+    a1, b1, f2, alab, blab = _plane_slice(X, Y, Z, F, plane)
+    x2, y2 = np.meshgrid(a1, b1, indexing="ij")
+
+    if cmap is None:
+        cmap = "RdBu_r" if symmetric else "viridis"
+    if vmax is None:
+        vmax = np.nanpercentile(np.abs(f2) if symmetric else f2, percentile)
+    if vmin is None:
+        vmin = -vmax if symmetric else np.nanpercentile(f2, 100 - percentile)
+    hi, lo = float(vmax), float(vmin)
+
+    # Contour boundaries land on round numbers, and the exact answers in this
+    # course are round numbers, so a field that is constant to round-off gets
+    # its cells sorted into the bands either side of a boundary and renders as
+    # structure that is not there. Both cases in Lab 2 are of that kind: the
+    # straining flow of Task 3 is div = +-5e-15 against a boundary at 0.0,
+    # drawn as faint red lobes, which is precisely the "it looks like a source"
+    # reading the task exists to refute; the rotation of Task 6 is curl = 2
+    # +- 1e-15 against a boundary at 2.0, drawn as a patchwork. Quantising at
+    # a billionth of the plotted range is far below any band and far above
+    # double-precision noise, so it removes the artefact and nothing else.
+    span = max(abs(hi), abs(lo))
+    if span > 0.0:
+        q = 1e-9 * span
+        f2 = np.round(f2 / q) * q
+
+    lv = np.linspace(lo, hi, levels)
+
+    created = ax is None
+    if created:
+        _, ax = plt.subplots(figsize=(5.4, 4.5))
+    cf = ax.contourf(x2, y2, np.clip(f2, lo, hi), levels=lv, cmap=cmap, extend="both")
+    ax.set_aspect("equal")                 # course rule: never distort a field plot
+    ax.set_xlabel(alab)
+    ax.set_ylabel(blab)
+    ax.set_title(title)
+    if colorbar:
+        ax.figure.colorbar(cf, ax=ax, label=label)
+    return ax, cf
+
+
+def show_field_slice(X, Y, Z, Ax, Ay, *, background=None, title="", label="",
+                     cmap="RdBu_r", density=1.3, symmetric=True, ax=None,
+                     percentile=98, colorbar=True, vmin=None, vmax=None,
+                     plane="z", levels=25, stream_color="k"):
+    """Streamlines of a vector field on a coordinate plane, over an optional
+    scalar background (typically the potential that generated it).
+
+    ``plane="z"`` cuts z = 0 and expects the (x, y) components; ``plane="y"``
+    cuts y = 0 and expects the (x, z) components -- pass ``Ax, Az`` there.
+
+    Returns ``(ax, cf)``. ``cf`` is the filled-contour mappable, or ``None``
+    if no background was given; pass ``colorbar=False`` on every panel of a
+    multi-panel figure and hand ``cf`` to ``fig.colorbar(cf, ax=axes, ...)``
+    to draw a single bar spanning the lot.
+
+    ``stream_color`` is the colour of the streamlines. Black reads well on a
+    diverging map, which is pale in the middle, and disappears on a sequential
+    one, which is dark at the bottom; pass ``"w"`` over ``inferno`` or
+    ``viridis``.
+    """
+    created = ax is None
+    if created:
+        _, ax = plt.subplots(figsize=(5.8, 4.8))
+
+    cf = None
+    if background is not None:
+        _, cf = show_scalar_slice(X, Y, Z, background, cmap=cmap, symmetric=symmetric,
+                                  percentile=percentile, ax=ax, colorbar=False,
+                                  vmin=vmin, vmax=vmax, plane=plane, levels=levels)
+
+    # streamplot needs 1-D increasing axes and arrays shaped (nb, na); our
+    # indexing='ij' arrays are (na, nb), hence the transposes.
+    x1, y1, u2, alab, blab = _plane_slice(X, Y, Z, Ax, plane)
+    _, _, v2, _, _ = _plane_slice(X, Y, Z, Ay, plane)
+    u = np.nan_to_num(u2).T
+    v = np.nan_to_num(v2).T
+    ax.streamplot(x1, y1, u, v, color=stream_color, linewidth=0.7,
+                  density=density, arrowsize=0.9)
+
+    ax.set_aspect("equal")
+    ax.set_xlim(x1.min(), x1.max())
+    ax.set_ylim(y1.min(), y1.max())
+    ax.set_xlabel(alab)
+    ax.set_ylabel(blab)
+    ax.set_title(title)
+    if colorbar and cf is not None:
+        ax.figure.colorbar(cf, ax=ax, label=label)
+    return ax, cf
+
+
+# --------------------------------------------------------------------------
+# Self-checks
+# --------------------------------------------------------------------------
+
+def check(label: str, ok: bool, hint: str = "") -> None:
+    """Report a pass, or raise with a hint about what to look at."""
+    if ok:
+        print(f"  [ok] {label}")
+    else:
+        raise AssertionError(f"{label} -- {hint}" if hint else label)
+
+
+def check_shape(label: str, arr, expected) -> None:
+    arr = np.asarray(arr)
+    check(f"{label}: shape {arr.shape}", arr.shape == tuple(expected),
+          f"expected {tuple(expected)}, got {arr.shape}. Did every array come "
+          f"from the same grid?")
+
+
+def check_close(label: str, got, want, rtol=0.05, where=None) -> None:
+    """Compare two arrays where both are finite (and where `where` is True)."""
+    got, want = np.asarray(got, float), np.broadcast_to(np.asarray(want, float), np.shape(got))
+    m = np.isfinite(got) & np.isfinite(want)
+    if where is not None:
+        m = m & where
+    if not m.any():
+        raise AssertionError(f"{label} -- nothing left to compare; the mask removed every point")
+    rel = np.abs(got[m] - want[m]) / np.maximum(np.abs(want[m]), 1e-30)
+    worst = float(np.max(rel))
+    check(f"{label}: worst error {worst:.2%}", worst < rtol,
+          f"worst relative error {worst:.2%} exceeds {rtol:.0%}. Check your "
+          f"np.gradient call -- did you pass dx, dy, dz, and in that order?")
+
+
+def check_abs(label: str, got, atol, where=None, hint: str = "") -> None:
+    """Compare an array against **zero**, on an absolute scale.
+
+    ``check_close`` divides by the expected value, so it cannot be pointed at
+    a quantity whose answer is exactly zero -- a solenoidal field, say. Give
+    this one a tolerance in the field's own units instead. ``atol`` is usually
+    a small fraction of the scale the field could have had: for a divergence
+    built from ``A ~ r`` on a grid of spacing ``dx``, anything below about
+    ``1e-9`` is round-off.
+    """
+    got = np.asarray(got, float)
+    m = np.isfinite(got)
+    if where is not None:
+        m = m & where
+    if not m.any():
+        raise AssertionError(f"{label} -- nothing left to compare")
+    worst = float(np.max(np.abs(got[m])))
+    check(f"{label}: worst |value| {worst:.2e}", worst < atol,
+          hint or f"worst deviation from zero, {worst:.2e}, exceeds {atol:.1e}")
+
+
+def check_scalar(label: str, got: float, want: float, rtol: float = 0.01,
+                 unit: str = "") -> None:
+    """Compare two single numbers and report the relative discrepancy."""
+    got, want = float(got), float(want)
+    rel = abs(got - want) / max(abs(want), 1e-30)
+    check(f"{label}: {got:.4g}{unit} vs {want:.4g}{unit} ({rel:.2%} apart)",
+          rel < rtol,
+          f"these should agree to better than {rtol:.0%}. Check the sign of "
+          f"each face, and that every face uses its own outward normal.")
diff --git a/book/1_gradient_divergence_curl/labs/week01_series_grad.md b/book/1_gradient_divergence_curl/labs/week01_series_grad.md
new file mode 100644
index 0000000..6b07d9e
--- /dev/null
+++ b/book/1_gradient_divergence_curl/labs/week01_series_grad.md
@@ -0,0 +1,882 @@
+---
+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 1: Series and Gradient
+
+:::{admonition} Computer lab
+:class: note
+
+A practical companion to the lectures on series, approximation and the gradient. Each task states a physical question, gives the steps, and ends with a self-check you can run. Plotting is supplied in the module `fwtools`, so that your effort goes into the physics rather than into rendering transparent isosurfaces.
+:::
+
+## Learning objectives
+
+By the end of this lab you should be able to:
+
+- **Truncate a series and quantify what the truncation costs.** Sum a geometric series, approximate it by its leading term, determine how many terms a given accuracy requires, and identify where a Taylor series ceases to converge.
+- **Read a gradient off a figure.** Show that $\nabla r = \hat{\boldsymbol{r}}$, that $\nabla f$ is normal to the level surfaces of $f$, and that $dp/dl = \lvert\nabla p\rvert\cos\psi$, so that the magnitude of the gradient is the maximum rate of change.
+- **Convert a potential into a field, and a field into a survey.** Apply $\boldsymbol{E} = -\nabla V$ and Ohm's law $\boldsymbol{J} = -\rho^{-1}\nabla V$, and map the potential and current density of a two-electrode DC resistivity measurement.
+
+:::{admonition} Two labs
+:class: note
+
+This is the first of two labs on the operators of this chapter. **Lab 2 covers the divergence and the curl**, and follows once those have been lectured.
+:::
+
+---
+
+## Part 0 — Setup
+
+Run this once. It contains no physics: it fetches two packages the browser lacks, locates `fwtools`, and defines the Coulomb constant.
+
+```{code-cell} ipython3
+# No physics above the k_e = ... line near the bottom.
+import sys, pathlib
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.constants import epsilon_0
+
+# --- 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/1_gradient_divergence_curl/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"
+# ---------------------------------------------------------------------------
+
+k_e = 1.0 / (4.0 * np.pi * epsilon_0)   # Coulomb constant, 8.99e9 V*m/C
+Q = 1e-9                                # 1 nC test charge
+
+print(f"epsilon_0 = {epsilon_0:.4e} F/m")
+print(f"k_e       = {k_e:.4e} V*m/C")
+```
+
+---
+
+## Part 1 — Series and truncation
+
+A physical quantity is often an infinite sum, of which only the first few terms are kept. Two questions follow: what the truncation costs, and whether the sum converges at all.
+
+### Task 1 — the bouncing ball
+
+A ball leaves the ground at $z=0$ with upward velocity $v_0$. Between bounces it is in free fall,
+
+$$ z(t) = v_0 t - \tfrac{1}{2}g t^2, $$
+
+so it returns to the ground after $T_0 = 2v_0/g$ having reached a height $H = v_0^2/2g$. At each bounce it loses a fraction $\gamma$ of its energy, so $v_n = \sqrt{1-\gamma}\;v_{n-1}$, and since flight time is proportional to launch speed,
+
+$$ T_n = (1-\gamma)^{n/2}\,T_0, \qquad T_0 = \sqrt{8H/g}. $$
+
+Fill in the three physical lines; the plotting is given.
+
+```{code-cell} ipython3
+g, v0, gamma = 9.81, 5.0, 0.1
+N = 12                                    # bounces to draw
+
+H  = ___                                  # peak height of the first flight
+T0 = ___                                  # duration of the first flight
+T  = T0 * ___                             # durations of bounces 0 .. N-1
+
+# --- given: draw one parabola per bounce ---
+t_start = np.concatenate(([0.0], np.cumsum(T)[:-1]))
+plt.figure(figsize=(9, 3.4))
+for Tn, t0 in zip(T, t_start):
+    tau = np.linspace(0, Tn, 200)
+    plt.plot(t0 + tau, (g*Tn/2)*tau - g*tau**2/2, "C0")
+plt.xlabel("$t$ [s]"); plt.ylabel("$z$ [m]"); plt.grid(alpha=0.3)
+plt.title(f"bouncing ball, $\\gamma$ = {gamma}")
+plt.show()
+
+# --- self-check (leave this alone) ---
+fw.check(f"H = {H:.4f} m", np.isclose(H, v0**2/(2*g)), "H = v0^2 / 2g")
+fw.check(f"T0 = {T0:.4f} s", np.isclose(T0, 2*v0/g), "T0 = 2 v0 / g")
+fw.check("T0 = sqrt(8H/g) too", np.isclose(T0, np.sqrt(8*H/g)))
+fw.check(f"{N} bounce durations, shrinking", len(T) == N and T[-1] < T[0])
+```
+
+:::{admonition} Solution — Task 1
+:class: dropdown
+
+```python
+H  = v0**2 / (2*g)
+T0 = 2*v0 / g
+T  = T0 * (1 - gamma)**(np.arange(N)/2)
+```
+:::
+
+The ball bounces for a total time
+
+$$ T_\infty = \sum_{m=0}^{\infty} T_m = T_0\sum_{m=0}^{\infty}\left(\sqrt{1-\gamma}\right)^{m} = \frac{\sqrt{8H/g}}{1-\sqrt{1-\gamma}}, $$
+
+a geometric series with ratio $\sqrt{1-\gamma}$. The ratio is smaller than 1 for any real bounce, so the sum is **finite**: infinitely many bounces, completed in about twenty seconds. The convergence condition matters, and Task 2 examines a series that fails it. For small $\gamma$, the expansion $\sqrt{1-\gamma}\approx 1-\gamma/2$ reduces the sum to
+
+$$ T_\infty \approx \sqrt{8H/g}\;\frac{2}{\gamma}. $$
+
+Both questions are answered below by measurement: **how accurate is that approximation**, and **how many bounces must be summed** before the running total reaches $T_\infty$?
+
+```{code-cell} ipython3
+rows = {}
+print(f"{'gamma':>7} {'T_inf':>9} {'approx':>9} {'error':>7} {'n for 99%':>10}")
+for gam in (0.5, 0.2, 0.1, 0.02):
+    T_inf  = ___                          # the exact sum, from the formula above
+    T_appr = ___                          # the small-gamma approximation
+
+    # --- given: how many bounces to reach 99% of T_inf ---
+    rows[gam] = (T_inf, T_appr)
+    cum = np.cumsum(T0 * (1 - gam)**(np.arange(4000)/2))
+    n99 = int(np.argmax(cum >= 0.99*T_inf)) + 1
+    print(f"{gam:>7.2f} {T_inf:>8.3f}s {T_appr:>8.3f}s "
+          f"{abs(T_appr-T_inf)/T_inf:>6.1%} {n99:>10}")
+
+# --- self-check (leave this alone) ---
+# The closed form against a brute-force sum of 5000 bounces: one number by
+# two routes, one of which never assumed the series converges.
+_summed = np.sum(T0 * (1 - 0.1)**(np.arange(5000)/2))
+fw.check(f"your T_inf at gamma = 0.1 ({rows[0.1][0]:.3f} s) equals the "
+         f"brute-force sum ({_summed:.3f} s)",
+         np.isclose(rows[0.1][0], _summed, rtol=1e-6))
+fw.check(f"your approximation overshoots by 2.6% there "
+         f"({rows[0.1][1]/rows[0.1][0] - 1:.2%})",
+         np.isclose(rows[0.1][1]/rows[0.1][0], 1.0263, rtol=1e-3))
+```
+
+:::{admonition} Solution — Task 1, continued
+:class: dropdown
+
+```python
+    T_inf  = np.sqrt(8*H/g) / (1 - np.sqrt(1 - gam))
+    T_appr = np.sqrt(8*H/g) * 2 / gam
+```
+:::
+
+:::{admonition} What the table says
+:class: important
+
+At $\gamma = 0.5$ the leading-term approximation is 17% wrong; at $\gamma = 0.02$ it is 0.5%. Keeping only the first term is a claim about the regime, not about the algebra, and it has to be justified case by case.
+
+The term count runs the other way. The more nearly elastic the ball, the more bounces the same accuracy requires: 14 at $\gamma = 0.5$, 456 at $\gamma = 0.02$. The approximation is cheapest exactly where the summation is most expensive. The same trade-off appears in every numerical method in this course.
+:::
+
+### Task 2 — where a Taylor series stops working
+
+A function that is smooth enough, and whose series sums back to it, can be written as a Taylor series about $x=0$,
+
+$$ f(x) = f(0) + x f'(0) + \tfrac{1}{2}x^2 f''(0) + \cdots . $$
+
+In practice the series is truncated after a few terms. Compare two cases:
+
+$$ \sin x = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \cdots, \qquad\qquad \frac{1}{1+x} = 1 - x + x^2 - x^3 + \cdots $$
+
+The two expansions look equally harmless. Add terms to each and compare.
+
+```{code-cell} ipython3
+x = np.linspace(-3, 3, 600)
+
+# term m of each series, as a function of x
+def sin_term(m, x):
+    return 0.0 if m % 2 == 0 else ___     # (-1)^((m-1)/2) x^m / m!   [math.factorial]
+
+def geo_term(m, x):
+    return ___                            # term m of 1 - x + x^2 - ...
+
+# --- given: exact curve plus four truncations, side by side ---
+fig, axes = plt.subplots(1, 2, figsize=(11, 4))
+for ax, (name, exact, term) in zip(axes, [
+        (r"$\sin x$",   np.sin,            sin_term),
+        (r"$1/(1+x)$",  lambda x: 1/(1+x), geo_term)]):
+    ax.plot(x, exact(x), "k", lw=2, label="exact")
+    for M in (2, 4, 8, 16):
+        ax.plot(x, sum(term(m, x) for m in range(M + 1)), lw=1, label=f"M = {M}")
+    ax.set_ylim(-3, 3); ax.set_xlabel("$x$"); ax.set_title(name)
+    ax.grid(alpha=0.3); ax.legend(fontsize=8)
+plt.tight_layout()
+plt.show()
+
+# --- self-check (leave this alone) ---
+_s21 = sum(sin_term(m, x) for m in range(21))
+_g_in  = sum(geo_term(m, 0.5) for m in range(40))
+_g_out = sum(geo_term(m, 1.5) for m in range(40))
+fw.check("21 terms reproduce sin(x) on -3 < x < 3", np.max(np.abs(_s21 - np.sin(x))) < 1e-6)
+fw.check(f"1/(1+x) converges at x = 0.5 ({_g_in:.4f} vs {1/1.5:.4f})", np.isclose(_g_in, 1/1.5))
+fw.check(f"1/(1+x) diverges at x = 1.5 (partial sum {_g_out:.2e})", abs(_g_out) > 1e3)
+```
+
+:::{admonition} Solution — Task 2
+:class: dropdown
+
+```python
+import math
+
+def sin_term(m, x):
+    return 0.0 if m % 2 == 0 else (-1)**((m-1)//2) * x**m / math.factorial(m)
+
+def geo_term(m, x):
+    return (-x)**m
+```
+:::
+
+:::{admonition} Radius of convergence
+:class: important
+
+$\sin x$ improves everywhere as terms are added. $1/(1+x)$ improves only inside $\lvert x\rvert < 1$; outside, each extra term makes the partial sum worse without limit, and at $x = 1.5$ the 40-term partial sum is off by millions.
+
+The series has a **radius of convergence** of 1, and no amount of computing power extends it. The radius is the distance from the expansion point to the nearest singularity of the function, here from $x=0$ to the pole at $x=-1$. The failure is invisible at $x = 0$ itself, where the function is smooth and the first few terms behave well. Expanding about $x = 1$ instead gives a radius of 2, because the pole is then twice as far away. **The radius is set by where the function is singular, not by its behaviour at the expansion point.**
+
+Compare Task 1, where more terms always helped and the only question was how many. Here, beyond $\lvert x\rvert = 1$, more terms are useless. Establishing which case applies is a prerequisite to any truncation.
+:::
+
+---
+
+## Part 2 — The distance function and its gradient
+
+All remaining parts use a single cube of sample points.
+
+```{code-cell} ipython3
+n, L = 61, 2.0                          # odd n, so the origin is a sample point
+axis = np.linspace(-L, L, n)            # one axis, shared by x, y and z
+X, Y, Z = np.meshgrid(axis, axis, axis, indexing="ij")
+dx = dy = dz = axis[1] - axis[0]
+
+c = n // 2                              # index of the origin
+# Mask reused by the self-checks. `interior` drops the two outermost cells,
+# so comparisons exclude the six faces of the box, where sampling is worst
+# and np.gradient has only one-sided neighbours.
+interior = np.zeros(X.shape, dtype=bool)
+interior[2:-2, 2:-2, 2:-2] = True
+
+print(f"grid shape {X.shape},  spacing {dx:.4f} m,  {X.size:,} sample points")
+print(f"X[i,j,k] = x[i]   ->   X[-1, 0, 0] = {X[-1, 0, 0]:.1f} m")
+```
+
+:::{admonition} Grid convention
+:class: tip
+
+**Resolution.** Every derivative on this page is a centred difference, so its error falls as $\Delta x^{2}$. Measured worst-case error against the analytic answer:
+
+| $n$ | $\Delta x$ [m] | $\lvert\nabla r\rvert$ | $\nabla(1/r)$ | $\nabla\cdot\boldsymbol{E}$ |
+| ---: | ---: | ---: | ---: | ---: |
+| 21 | 0.200 | 4.1% | 12.5% | 9.1% |
+| 41 | 0.100 | 1.8% | 3.3% | 2.4% |
+| **61** | **0.067** | **0.8%** | **1.6%** | **1.1%** |
+| 81 | 0.050 | 0.5% | 1.0% | 0.6% |
+
+These are **worst cases** over the region each self-check tests, not averages, and they are what fixes the tolerances. Halving $\Delta x$ reduces the last two columns by close to the factor of four that second-order accuracy predicts (12.5 → 3.3, 9.1 → 2.4). The first column falls by only 2.3. The $|\nabla r|$ error grows towards the source, so the worst sample in the band $0.4 < r < 1.6$ m is whichever one sits nearest the inner edge, and that sample moves when `n` changes. A worst case taken over a boundary that the grid keeps redrawing does not form a smooth sequence. The convergence cell at the end of Lab 2 measures a fixed quantity instead and does recover the factor of four.
+
+$n = 61$ was chosen from this table as the coarsest grid that keeps every task under 2%; each 3-D figure it produces is about 1.5 MB. **If you change `n`, keep it at 41 or above.** The self-checks below allow 5%, and $n = 31$ already fails Task 6 at 6.7%.
+
+Two properties of this cube matter later. It is a finite window on fields that extend to infinity: the largest closed surface in Lab 2 sits only 0.6 m inside the outer face. And $z$ points **up**, as in an ordinary right-handed frame, whereas Part 4 works in the ground, where $z$ points downwards by Earth-science convention. Neither choice is more correct; stating which one is in use is what matters.
+
+The grid is built with `indexing='ij'`, so axis 0 is $x$, axis 1 is $y$, axis 2 is $z$.
+
+1. **Derivatives come back in coordinate order:** `np.gradient(f, dx, dy, dz)` returns $\partial f/\partial x$, $\partial f/\partial y$, $\partial f/\partial z$. No transposes.
+2. **Always pass the spacings.** Omit them and the derivative is silently wrong by a factor of $1/\Delta x = 15$.
+
+NumPy's default is `indexing='xy'`, which returns the $y$-derivative first. That single difference accounts for a large share of numerical field bugs.
+:::
+
+The simplest scalar field is the distance to a point:
+
+$$ r(x,y,z) = \sqrt{(x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2} $$
+
+One number at every location in space: no charge, no potential, and no units beyond metres.
+
+This is the **spherical** radial coordinate $r$, the distance from a point. The cylindrical radius $\varrho$, the distance from an axis, is a different quantity. The equations on this page use $r$, and the code calls it `r`.
+
+### Task 3 — build the distance field
+
+**The question:** what do the surfaces of constant $r$ look like, and where do they crowd together? 
+
+The source must be movable: Task 8 places two of them at different points, so write the offsets in now rather than hard-coding the origin.
+
+```{code-cell} ipython3
+# Task 3 -- distance from a source at (x0, y0, z0) to every point of the grid.
+
+def distance_to(X, Y, Z, x0=0.0, y0=0.0, z0=0.0):
+    return ___                          # root of the sum of three squares
+
+
+r = ___                                 # call it: one source, at the origin
+
+# --- self-check (leave this alone) ---
+fw.check_shape("r", r, X.shape)
+fw.check("r = 0 at the origin", np.isclose(r[c, c, c], 0.0))
+fw.check("r = 2 m at (2,0,0)", np.isclose(r[-1, c, c], 2.0))
+fw.check("r = 2 m at (0,2,0)", np.isclose(r[c, -1, c], 2.0))
+fw.check("the source can be moved off the origin",
+         np.isclose(distance_to(X, Y, Z, 1.0, 0.0, 0.0)[c, c, c], 1.0),
+         "x0, y0, z0 have to appear in the expression -- Task 8 needs them")
+```
+
+:::{admonition} Solution — Task 3
+:class: dropdown
+
+```python
+def distance_to(X, Y, Z, x0=0.0, y0=0.0, z0=0.0):
+    return np.sqrt((X - x0)**2 + (Y - y0)**2 + (Z - z0)**2)
+
+
+r = distance_to(X, Y, Z)
+```
+:::
+
+A surface on which $r$ takes one fixed value is an **isosurface**, or level set, the three-dimensional analogue of a contour line on a map. Evenly spaced values of $r$ give evenly spaced shells: the distance function has no preferred radius.
+
+```{code-cell} ipython3
+fw.show_isosurfaces(X, Y, Z, r, levels=[0.5, 1.0, 1.5], label="r  [m]",
+                    title="Isosurfaces of the distance function r")
+```
+
+
+### Task 4 — the gradient of the distance
+
+Do this one on paper first. Differentiating $r = \sqrt{x^2+y^2+z^2}$ by the chain rule,
+
+$$ \frac{\partial r}{\partial x} = \frac{x}{r}, \qquad \frac{\partial r}{\partial y} = \frac{y}{r}, \qquad \frac{\partial r}{\partial z} = \frac{z}{r} $$
+
+so, collecting the three components,
+
+$$ \nabla r \;=\; \frac{\partial r}{\partial x}\hat{\boldsymbol{x}} + \frac{\partial r}{\partial y}\hat{\boldsymbol{y}} + \frac{\partial r}{\partial z}\hat{\boldsymbol{z}} \;=\; \frac{x\,\hat{\boldsymbol{x}} + y\,\hat{\boldsymbol{y}} + z\,\hat{\boldsymbol{z}}}{r} \;=\; \hat{\boldsymbol{r}} $$
+
+The last step is the definition of the outward unit radial vector: $\hat{\boldsymbol{r}}$ is the position vector divided by its own length. So $\nabla r$ is a **unit** vector pointing **away** from the source, with both direction and magnitude known in advance.
+
+The cell below tests whether a **finite-difference gradient** on a grid reproduces that. Two measurements: the magnitude, which should be 1, and the projection $\nabla r \cdot \hat{\boldsymbol{r}}$, which recovers the full magnitude only if the gradient is purely radial, with no component along the sphere.
+
+```{code-cell} ipython3
+# The outward unit radial vector, used again later.
+rs = np.maximum(r, 1e-12)                 # 0/0 at the source is not a lesson
+rhx, rhy, rhz = X / rs, Y / rs, Z / rs
+
+# Task 4
+#   1. grad r, as three components.
+#   2. Its magnitude.
+#   3. Its projection onto r-hat.
+#   4. Draw it, rotate the figure, and compare with the spheres above.
+
+grx, gry, grz = ___                       # all three spacings, in order
+
+grad_r_mag = ___                          # the length of that vector
+
+radial_part = ___                         # its projection onto (rhx, rhy, rhz)
+
+fw.show_cones(X, Y, Z, grx, gry, grz, step=8, label="|∇r|", unit="-",
+              title="grad r -- unit vectors pointing away from the source")
+
+# --- self-check (leave this alone) ---
+band = (r > 0.4) & (r < 1.6)
+fw.check_shape("grad r (x-component)", grx, X.shape)
+fw.check_close("|grad r| = 1 everywhere", grad_r_mag, 1.0, rtol=0.05, where=band)
+fw.check_close("grad r is purely radial", radial_part, 1.0, rtol=0.05, where=band)
+# The two checks above are the same measurement for THIS field, so they pass
+# or fail together. This one is independent: it compares the three components
+# against r-hat separately, catching a gradient of the right length but the
+# wrong direction.
+fw.check(f"grad r = r-hat, componentwise (worst "
+         f"{np.nanmax(np.abs(np.stack([grx-rhx, gry-rhy, grz-rhz]))[:, band]):.3f} "
+         f"of a unit vector)",
+         np.nanmax(np.abs(np.stack([grx - rhx, gry - rhy, grz - rhz]))[:, band]) < 0.05)
+```
+
+:::{admonition} Solution — Task 4
+:class: dropdown
+
+```python
+grx, gry, grz = np.gradient(r, dx, dy, dz)
+grad_r_mag = np.sqrt(grx**2 + gry**2 + grz**2)
+radial_part = grx * rhx + gry * rhy + grz * rhz
+
+print(f"|grad r| median in 0.4 < r < 1.6 m : "
+      f"{np.median(grad_r_mag[(r > 0.4) & (r < 1.6)]):.4f}")
+```
+:::
+
+:::{admonition} What the algebra means
+:class: important
+
+$\lvert\nabla r\rvert = 1$ needs no calculus: move one metre directly away from the source and the distance to it grows by one metre, so the steepest rate of change of $r$ is 1 m/m everywhere. A gradient carries the direction of steepest increase and a length equal to that rate, here outward and 1.
+
+The radial check fixes the other half: moving along a sphere does not change $r$, so the gradient has no component there. **$\nabla f$ is normal to the level surfaces of $f$** for every scalar field, not only this one.
+
+:::
+
+### Task 5 — the rate of change in an arbitrary direction
+
+The direction of the gradient is settled: steepest increase, normal to the level surface. Its magnitude is the untested claim. It follows from
+
+$$ dp = (\nabla p)\cdot d\boldsymbol{l} = \lvert\nabla p\rvert\,\lvert d\boldsymbol{l}\rvert\cos\psi
+\qquad\Longrightarrow\qquad
+\frac{dp}{dl} = \lvert\nabla p\rvert\cos\psi, $$
+
+where $d\boldsymbol{l}$ is a small step in any chosen direction, $dl = \lvert d\boldsymbol{l}\rvert$ is its length, and $\psi$ is the angle between the step and the gradient. The step is written $d\boldsymbol{l}$ rather than $d\boldsymbol{r}$ because $r$ already denotes the distance from the origin on this page.
+
+Two testable consequences: the rate of change in any direction is $\lvert\nabla p\rvert\cos\psi$, and it never exceeds $\lvert\nabla p\rvert$, which is reached only at $\psi = 0$.
+
+Measure it. At one point, step a short distance $\varepsilon$ along many unit vectors $\hat{\boldsymbol{u}}$ and compare each measured rate with the prediction.
+
+```{code-cell} ipython3
+p_field = 1.0 / np.maximum(r, 0.25)       # any scalar field will do
+gpx, gpy, gpz = np.gradient(p_field, dx, dy, dz)
+
+ip, jp, kp = 40, 36, 34                   # one sample point, off-axis
+gvec = np.array([gpx[ip, jp, kp], gpy[ip, jp, kp], gpz[ip, jp, kp]])
+point = np.array([axis[ip], axis[jp], axis[kp]])
+
+def p_exact(q):
+    return 1.0 / np.linalg.norm(q)        # the same field, evaluated anywhere
+
+# Task 5 -- fill in the four blanks; the plotting is given.
+grad_mag = ___                            # |grad p| at the point, from gvec
+
+rng = np.random.default_rng(0)
+eps = 1e-4
+cosines, rates = [], []
+for _ in range(200):
+    u = rng.normal(size=3)
+    u = ___                               # make it a UNIT vector
+    cosines.append(___)                   # cos(psi) = u . gvec / |grad p|
+    rates.append(___)                     # centred difference of p_exact
+                                          # along u, step eps, over 2*eps
+cosines, rates = np.asarray(cosines), np.asarray(rates)
+
+# --- given: measurements against the predicted straight line ---
+plt.figure(figsize=(5.6, 4.4))
+plt.scatter(cosines, rates, s=12, alpha=0.6, label="measured")
+cs = np.linspace(-1, 1, 50)
+plt.plot(cs, grad_mag*cs, "k", lw=1.5, label=r"$|\nabla p|\cos\psi$")
+plt.xlabel(r"$\cos\psi$")
+plt.ylabel(r"$dp/dl$  [m$^{-2}$]")
+plt.legend(); plt.grid(alpha=0.3)
+plt.show()
+
+# --- self-check (leave this alone) ---
+slope = float(np.polyfit(cosines, rates, 1)[0])
+fw.check_scalar("fitted slope = |grad p|", slope, grad_mag, rtol=0.01)
+fw.check("no direction beats |grad p|", np.max(np.abs(rates)) <= grad_mag * 1.001)
+```
+
+:::{admonition} Solution — Task 5
+:class: dropdown
+
+```python
+grad_mag = float(np.linalg.norm(gvec))
+
+# ... and inside the loop:
+    u = u / np.linalg.norm(u)
+    cosines.append(float(u @ gvec) / grad_mag)
+    rates.append((p_exact(point + eps*u) - p_exact(point - eps*u)) / (2*eps))
+```
+:::
+
+:::{admonition} The magnitude, measured
+:class: important
+
+Every measured rate lies on the line. Three readings of the same figure:
+
+- **At $\cos\psi = 1$** the step is straight up the gradient and the rate equals $\lvert\nabla p\rvert$. No direction exceeds it, which is the content of *steepest*, now measured rather than asserted.
+- **At $\cos\psi = 0$** the step lies in the level surface and $p$ does not change. This is the normality result of Task 4, recovered by a second route.
+- **At $\cos\psi = -1$** the rate is $-\lvert\nabla p\rvert$, the steepest descent, which is the direction $\boldsymbol{E} = -\nabla V$ selects in Part 3.
+
+One vector carries both a direction and a rate; the cosine gives the rate along any other direction.
+:::
+
+---
+
+## Part 3 — The inverse distance
+
+The function that appears in the physics is not the distance but its reciprocal,
+
+$$ f(r) = \frac{1}{r}, \qquad\text{so}\qquad \nabla f = \frac{d}{dr}\!\left(\frac{1}{r}\right)\hat{\boldsymbol{r}} = -\frac{1}{r^{2}}\,\hat{\boldsymbol{r}} $$
+
+The isosurfaces are the same spheres, since $f$ is constant wherever $r$ is constant, but the ordering is inverted: $f$ is largest near the source and decays to zero far away. Predict the effect on the arrows, check the prediction against the formula above, then measure it.
+
+### Task 6 — the gradient of the inverse distance
+
+```{code-cell} ipython3
+# The mask keeps the singularity at r = 0 off the grid: everything within
+# 0.25 m of the source becomes NaN and is not measured.
+r_masked = np.where(r < 0.25, np.nan, r)
+f = 1.0 / r_masked
+
+# Task 6 -- two blanks. Predict the direction before you look at the figure.
+fx, fy, fz = ___                          # grad f
+f_mag = ___                               # its magnitude, to compare with 1/r^2
+
+# --- given: the numbers, then the picture ---
+for rr in (0.6, 1.0, 1.5):
+    i = int(np.argmin(np.abs(X[:, 0, 0] - rr)))
+    print(f"r = {rr:.1f} m :  |grad f| = {f_mag[i, c, c]:8.4f}   1/r^2 = {1/rr**2:8.4f}")
+
+# normalise=True draws every arrow the same length, so the figure carries
+# direction only. The magnitude moves into the colour, on a log scale,
+# because the drawn arrows span a factor of 62.
+fw.show_cones(X, Y, Z, fx, fy, fz, step=8, normalise=True,
+              label="|∇(1/r)|", unit="m<sup>-2</sup>",
+              title="grad(1/r) -- pointing back towards the source")
+
+# --- self-check (leave this alone) ---
+outside = (r > 0.5) & interior            # `interior` was built in Part 2
+fw.check_close("|grad(1/r)| = 1/r^2", f_mag, 1.0 / r_masked**2, rtol=0.05, where=outside)
+fw.check("grad(1/r) points inward at (1,0,0)", fx[-1 - 15, c, c] < 0)
+```
+
+:::{admonition} Solution — Task 6
+:class: dropdown
+
+```python
+fx, fy, fz = np.gradient(f, dx, dy, dz)
+f_mag = np.sqrt(fx**2 + fy**2 + fz**2)
+```
+:::
+
+:::{admonition} The gradient always points towards increase
+:class: important
+
+The arrows have reversed. Same spheres, same source, opposite direction:
+
+$$ \nabla r = +\hat{\boldsymbol{r}}, \qquad\qquad \nabla\!\left(\frac{1}{r}\right) = -\frac{1}{r^{2}}\,\hat{\boldsymbol{r}} $$
+
+Nothing about space changed; what changed is **which way the function climbs**. The steepness changed as well: $1/r$ climbs faster as the source is approached, so its gradient grows as $1/r^2$ instead of staying at 1.
+
+A gradient encodes nothing about sources, sinks, charges or fields. It encodes only the uphill direction and the rate along it.
+:::
+
+### Task 7 — from geometry to physics
+
+The physics enters as a single minus sign. The electric potential of a point charge $Q$ is the inverse-distance function with a constant in front,
+
+$$ V(r) = \frac{1}{4\pi\varepsilon_0}\frac{Q}{r}\quad[\text{V}], $$
+
+and the electric field is *defined* as
+
+$$ \boldsymbol{E} = -\nabla V \quad[\text{V/m}]. $$
+
+$\nabla V$ points inward, uphill towards the charge. The minus sign reverses it, so **the field points downhill**, which is the direction a positive test charge released from rest would move, losing potential energy as it goes.
+
+```{code-cell} ipython3
+V = k_e * Q / r_masked
+
+# Task 7 -- two blanks. Mind the minus sign.
+Ex, Ey, Ez = ___                          # E = -grad V
+E_mag = ___
+
+# --- given: against the analytic k_e*Q/r^2, then the picture ---
+for rr in (0.6, 1.0, 1.5):
+    i = int(np.argmin(np.abs(X[:, 0, 0] - rr)))
+    print(f"r = {rr:.1f} m :  |E| = {E_mag[i, c, c]:8.3f} V/m   "
+          f"analytic = {k_e*Q/rr**2:8.3f} V/m")
+
+fw.show_cones(X, Y, Z, Ex, Ey, Ez, step=8, normalise=True,
+              label="|<b>E</b>|", unit="V/m",
+              title="E = -grad V for a positive point charge")
+
+# --- self-check (leave this alone) ---
+fw.check_close("|E| = Q/(4 pi eps0 r^2)", E_mag, k_e * Q / r_masked**2,
+               rtol=0.05, where=outside)
+fw.check("E points outward at (1,0,0)", Ex[-1 - 15, c, c] > 0)
+```
+
+:::{admonition} Solution — Task 7
+:class: dropdown
+
+```python
+dVdx, dVdy, dVdz = np.gradient(V, dx, dy, dz)
+Ex, Ey, Ez = -dVdx, -dVdy, -dVdz
+E_mag = np.sqrt(Ex**2 + Ey**2 + Ez**2)
+```
+:::
+
+:::{admonition} Why the potential is worth defining
+:class: tip
+
+$V$ is a scalar: one number per point, with no direction to track. $\boldsymbol{E}$ is a vector: three numbers. Any operation carried out once on $V$ and then differentiated is cheaper, in arithmetic and in bookkeeping, than the same operation carried out three times on $\boldsymbol{E}$.
+
+Part 4 is the first case where this matters.
+:::
+
+---
+
+## Part 4 — Two sources: superposition
+
+A single charge is spherically symmetric. Two are not:
+
+$$ V_{\text{total}} = \frac{1}{4\pi\varepsilon_0}\left(\frac{Q_1}{r_1} + \frac{Q_2}{r_2}\right) $$
+
+**Superposition** of potentials is the addition of two numbers at every point, because $V$ is a scalar. Superposing the two fields instead requires a vector sum at every point of the cube.
+
+Since $\nabla$ is linear, $-\nabla(V_1 + V_2) = \boldsymbol{E}_1 + \boldsymbol{E}_2$ exactly. The efficient route is therefore to **add the potentials and take a single gradient at the end**, with no loss of accuracy.
+
+### Task 8 — build a dipole
+
+```{code-cell} ipython3
+# Distances to the two charges. +Q sits at x = +d/2, -Q at x = -d/2, the same
+# placement Task 9 gives the current source and sink, so the two figures can
+# be compared directly. The guard trips only if a grid point lands exactly on
+# a charge; at n = 61 none does, so nothing is masked and the full field is
+# shown. Raise it if you change the grid.
+d_sep = 1.0                               # charge separation [m]
+r_plus = np.where(distance_to(X, Y, Z, +d_sep/2, 0.0, 0.0) < 0.01, np.nan,
+                  distance_to(X, Y, Z, +d_sep/2, 0.0, 0.0))
+r_minus = np.where(distance_to(X, Y, Z, -d_sep/2, 0.0, 0.0) < 0.01, np.nan,
+                   distance_to(X, Y, Z, -d_sep/2, 0.0, 0.0))
+
+# Task 8 -- two blanks.
+V_dip = ___                               # superpose: +Q over r_plus, -Q over r_minus
+Ex_d, Ey_d, Ez_d = ___                    # ONE gradient of the sum, negated
+
+# --- given: the z = 0 plane, potential as colour, field as streamlines ---
+fw.show_field_slice(X, Y, Z, Ex_d, Ey_d, background=V_dip,
+                    title="Source and sink: potential (colour) and field lines",
+                    label="$V$ [V]")
+plt.show()
+
+# --- self-check (leave this alone) ---
+mid = np.abs(X) < 1e-9                    # the plane x = 0, halfway between them
+fw.check_shape("V_dip", V_dip, X.shape)
+fw.check("V = 0 on the mid-plane",
+         np.nanmax(np.abs(V_dip[mid])) < 1e-6 * np.nanmax(np.abs(V_dip)))
+fw.check("E on the mid-plane points from + to -", np.nanmean(Ex_d[mid]) < 0)
+```
+
+:::{admonition} Solution — Task 8
+:class: dropdown
+
+```python
+V_dip = k_e * Q / r_plus + k_e * (-Q) / r_minus
+
+dVx, dVy, dVz = np.gradient(V_dip, dx, dy, dz)
+Ex_d, Ey_d, Ez_d = -dVx, -dVy, -dVz
+```
+:::
+
+:::{admonition} The mid-plane
+:class: tip
+
+At $x = 0$ the potential is **exactly zero**, while the field is at its strongest, pointing straight from the positive charge to the negative one, here along $-\hat{\boldsymbol{x}}$ because $+Q$ sits on the right.
+
+The field is the slope of the potential, not its value: terrain at sea level can still be steep. The figure also shows the streamlines crossing the coloured contours at right angles everywhere, which is the normality result of Task 4 appearing in a field that was not constructed radially.
+:::
+
+### The far field of the dipole
+
+$V_{\text{dip}}$ is not a series; it is two exact terms. Viewed from far enough away, however, the two charges are no longer resolvable, and what survives is a **truncation**.
+
+Expand $1/r_\pm$ in powers of $d/r$ and add. The leading terms are equal and opposite, since the pair carries no net charge, and the first surviving term is
+
+$$ V \;\approx\; \frac{1}{4\pi\varepsilon_0}\frac{\boldsymbol{p}\cdot\hat{\boldsymbol{r}}}{r^{2}}, \qquad \boldsymbol{p} = Q d\,\hat{\boldsymbol{x}}, $$
+
+with the **dipole moment** $\boldsymbol{p}$ pointing from the negative charge to the positive one. Every discarded term is smaller by a further factor of $(d/r)^2$. This is the question of Task 1 asked of distance rather than of term count: at what range is one term enough?
+
+```{code-cell} ipython3
+# --- given: exact against the one-term far field, along the +x axis ---
+p_mom = Q * d_sep                              # dipole moment [C m]
+r_ff = np.logspace(np.log10(0.8), np.log10(60), 2000)
+V_ex = k_e * Q * (1/np.abs(r_ff - d_sep/2) - 1/np.abs(r_ff + d_sep/2))
+V_ff = k_e * p_mom / r_ff**2                   # p . r-hat = p on the axis
+err_ff = np.abs(V_ff - V_ex) / np.abs(V_ex)
+
+plt.figure(figsize=(5.8, 4.2))
+plt.loglog(r_ff / d_sep, err_ff, "k", lw=1.6)
+for tol, colour in ((0.10, "C1"), (0.01, "C2"), (0.001, "C3")):
+    r_ok = r_ff[np.argmax(err_ff < tol)] / d_sep
+    plt.axhline(tol, color=colour, lw=0.8, ls=":")
+    plt.plot([r_ok], [tol], "o", color=colour, ms=5)
+    print(f"  one term is good to {tol:6.1%} beyond r = {r_ok:5.1f} separations")
+plt.xlabel("$r$ / separation $d$"); plt.ylabel("relative error of the one-term form")
+plt.grid(alpha=0.3, which="both"); plt.title("How far is far?")
+plt.show()
+
+# --- self-check (leave this alone) ---
+fw.check("the far-field error falls as (d/r)^2",
+         np.isclose(np.polyfit(np.log(r_ff[r_ff > 10]), np.log(err_ff[r_ff > 10]), 1)[0],
+                    -2.0, atol=0.05))
+```
+
+:::{admonition} The same question as the bouncing ball
+:class: important
+
+Two decades of accuracy cost a factor of ten in distance: 10% at $1.6\,d$, 1% at $5\,d$, 0.1% at $16\,d$. This is the $(d/r)^2$ law, and the factor between successive rows is $\sqrt{10}\approx 3.2$.
+
+Compare Task 1, where 1% accuracy required 88 terms, and the count grew as the ball became more elastic. Here the controlling variable is a distance rather than a term count, and the requirement grows the closer the observation point. In both cases the truncation is only as good as the regime, and in both cases the regime can be established by measurement.
+
+This single term is why a compass works. A magnet has a complicated field close up; at a metre it is a dipole and nothing else, which is why the Earth's field is written as the single term used in Lab 2's Task 2.
+:::
+
+The same object in three dimensions, with positive and negative equipotential surfaces drawn transparent:
+
+```{code-cell} ipython3
+lobe = np.nanpercentile(np.abs(V_dip), 97)
+fw.show_isosurfaces(X, Y, Z, np.nan_to_num(V_dip), levels=[-lobe, -lobe/3, lobe/3, lobe],
+                    colorscale="RdBu", reversescale=True, opacity=0.3, label="V  [V]",
+                    title="Equipotential surfaces of a dipole")
+```
+
+### Task 9 — the same mathematics as a geophysical survey
+
+Task 8 was two charges in vacuum. The mathematics below is identical; the physics is not.
+
+Drive a current $I$ into the ground through one electrode and extract it through another, a distance $a$ away. Air does not conduct, so in ground of resistivity $\rho$ the current spreads through the **lower half-space only**, and each electrode contributes $\rho I/2\pi r$ rather than $\rho I / 4\pi r$. Superposition gives
+
+$$ V(x,y,z) = \frac{\rho I}{2\pi}\left(\frac{1}{\lvert\boldsymbol{r}-\boldsymbol{a}/2\rvert} - \frac{1}{\lvert\boldsymbol{r}+\boldsymbol{a}/2\rvert}\right), \qquad z \ge 0 \ \text{(down into the ground)}. $$
+
+The field follows as before, $\boldsymbol{E} = -\nabla V$, and Ohm's law in local form turns it into a **current density**:
+
+$$ \boldsymbol{J} = \rho^{-1}\boldsymbol{E} = -\rho^{-1}\nabla V \quad [\text{A}/\text{m}^2]. $$
+
+This is a DC resistivity survey, a standard near-surface geophysical measurement. Map it two ways: on the ground surface, where the electrodes are planted, and on a vertical section cut between them.
+
+:::{admonition} $\rho$ means something else here
+:class: warning
+
+In this task $\rho$ is the **electrical resistivity** in Ω·m. In Lab 2's Task 4 it is a charge density in C/m³, written $\rho_v$ to keep the two apart. The symbol is overloaded throughout the subject; the units identify which is meant.
+:::
+
+The ground is a half-space, so this task needs its own grid: $x$ and $y$ still run from $-L$ to $L$, but $z$ runs from $0$ at the surface **downwards**, following the Earth-science convention.
+
+A real electrode is a metal stake, not a mathematical point: a conductor of finite radius $r_{\text{el}}$ held at one potential over its whole surface. Flooring the distance at $r_{\text{el}}$ models it that way and keeps $1/r$ bounded. Nothing is masked, no sample is discarded, and every derivative below acts on a field that is finite everywhere.
+
+```{code-cell} ipython3
+rho, I, a_sep = 100.0, 1.0, 1.0           # ohm.m, ampere, electrode spacing [m]
+r_el = 0.12                               # electrode radius [m]
+
+axis_g = np.linspace(-2.0, 2.0, 81)       # x and y, across the survey line
+depth  = np.linspace(0.0, 2.0, 51)        # z, down into the ground
+Xg, Yg, Zg = np.meshgrid(axis_g, axis_g, depth, indexing="ij")
+dxg = axis_g[1] - axis_g[0]
+dzg = depth[1] - depth[0]
+print(f"dxg = {dxg:.3f} m,  dzg = {dzg:.3f} m   <- deliberately not equal")
+
+def dist_to(x0):
+    """Distance to an electrode at (x0, 0, 0), floored at its own radius."""
+    return np.maximum(np.sqrt((Xg - x0)**2 + Yg**2 + Zg**2), r_el)
+
+# Task 9 -- two blanks. This is Task 8 again, in different clothes.
+#   V:  the formula above, source at x = +a_sep/2, sink at x = -a_sep/2.
+#       dist_to floors the distance at the electrode radius, so there is
+#       nothing to mask and nothing to nan_to_num.
+#   J:  -grad(V)/rho. Pass dxg, dxg, dzg. On this grid z is spaced
+#       differently from x and y, and passing dxg three times costs 6.5% on
+#       the current measured in the next cell, enough to fail its check.
+
+V_dc = ___
+Jx, Jy, Jz = ___
+
+# --- given: the survey, both panels on one colour scale and one colorbar ---
+# plane="z" is the ground surface; plane="y" is the vertical section, where
+# the in-plane components are (Jx, Jz), not (Jx, Jy).
+vm = float(np.nanpercentile(np.abs(V_dc[:, :, 0]), 98))
+fig, axes = plt.subplots(2, 1, figsize=(7.2, 9.2))
+for ax_, comps, pl, ttl in ((axes[0], (Jx, Jy), "z", "a) ground surface, $z=0$"),
+                            (axes[1], (Jx, Jz), "y", "b) vertical section, $y=0$")):
+    _, cf = fw.show_field_slice(Xg, Yg, Zg, *comps, background=V_dc, ax=ax_,
+                                plane=pl, vmin=-vm, vmax=vm, colorbar=False,
+                                density=1.2, title=ttl)
+axes[1].invert_yaxis()                    # depth increases downwards
+fig.colorbar(cf, ax=axes, label="$V$ [V]", fraction=0.05, pad=0.03)
+plt.show()
+
+# --- self-check (leave this alone) ---
+mid_dc = np.abs(Xg) < 1e-9
+fw.check("V is finite everywhere -- no holes in the model",
+         np.all(np.isfinite(V_dc)) and np.all(np.isfinite(Jx)))
+fw.check("V = 0 on the mid-plane between the electrodes",
+         np.nanmax(np.abs(V_dc[mid_dc])) < 1e-6 * np.nanmax(np.abs(V_dc)))
+fw.check("current flows from the source towards the sink at the surface",
+         np.nanmean(Jx[mid_dc]) < 0)
+```
+
+:::{admonition} Solution — Task 9
+:class: dropdown
+
+```python
+V_dc = rho * I / (2*np.pi) * (1/dist_to(+a_sep/2) - 1/dist_to(-a_sep/2))
+
+gVx, gVy, gVz = np.gradient(V_dc, dxg, dxg, dzg)
+Jx, Jy, Jz = -gVx/rho, -gVy/rho, -gVz/rho
+```
+
+A presentation point worth reusing: `show_field_slice` returns `(ax, cf)`, so passing `colorbar=False` on both panels and handing the mappable `cf` to `fig.colorbar(..., ax=axes)` draws **one** bar beside the pair. Two bars carrying identical numbers are clutter, and they suggest to the reader that the scales differ.
+:::
+
+Now use the field as an instrument. All the current injected at one electrode must cross any closed surface drawn around it, since there is nowhere else for it to go. Test that.
+
+```{code-cell} ipython3
+# The five faces of a box buried in the ground around one electrode. The top
+# face is deliberately absent: it lies in the surface z = 0, where no current
+# crosses into the air, so its contribution is zero by physics.
+def buried_box_current(xc, hw=0.3):
+    i0 = int(np.argmin(np.abs(axis_g - (xc - hw))))
+    i1 = int(np.argmin(np.abs(axis_g - (xc + hw))))
+    j0 = int(np.argmin(np.abs(axis_g + hw)))
+    j1 = int(np.argmin(np.abs(axis_g - hw)))
+    k1 = int(np.argmin(np.abs(depth - hw)))
+    sx, sy, sz = slice(i0, i1+1), slice(j0, j1+1), slice(0, k1+1)
+    return (fw.area_integral(Jx[i1, sy, sz], dxg, dzg) - fw.area_integral(Jx[i0, sy, sz], dxg, dzg)
+          + fw.area_integral(Jy[sx, j1, sz], dxg, dzg) - fw.area_integral(Jy[sx, j0, sz], dxg, dzg)
+          + fw.area_integral(Jz[sx, sy, k1], dxg, dxg))
+
+for xc, name in ((+a_sep/2, "source"), (-a_sep/2, "sink")):
+    print(f"current out of a box around the {name:6s}: {buried_box_current(xc):+7.4f} A")
+print(f"                            injected: {I:+7.4f} A")
+
+# --- self-check (leave this alone) ---
+fw.check_scalar("box around the source carries I", buried_box_current(+a_sep/2), I, rtol=0.01, unit=" A")
+fw.check_scalar("box around the sink carries -I", buried_box_current(-a_sep/2), -I, rtol=0.01, unit=" A")
+```
+
+:::{admonition} Why five faces and not six?
+:class: important
+
+The box is closed by the ground surface itself. Air does not conduct, so $J_z = 0$ at $z=0$. This is a **boundary condition**, true by physics, and not a quantity to be measured.
+
+Measuring it anyway is instructive. `np.gradient` has no neighbour above $z=0$, so it falls back to a one-sided difference and reports a spurious $J_z$ averaging $+0.25$ A/m² over the top of the box, apparently current entering from the air. The outward normal on that face is $-\hat{\boldsymbol{z}}$, so the face enters the sum as $-0.088$ A and reduces the box total from $1.003$ A to $0.915$ A, an **8.5% error** on a result that is otherwise accurate to 0.3%.
+
+The rule generalises well beyond this lab: **impose a boundary condition you know exactly, rather than asking a finite-difference stencil to recover it.** Numerical derivatives are least reliable where the domain stops.
+:::
+
+---
+
+:::{admonition} End of Lab 1
+:class: note
+
+Parts 1–4 cover the gradient. **The divergence and the curl continue in Lab 2**, which is lectured next.
+:::
diff --git a/book/1_gradient_divergence_curl/labs/week02_div_curl.md b/book/1_gradient_divergence_curl/labs/week02_div_curl.md
new file mode 100644
index 0000000..ac0e3df
--- /dev/null
+++ b/book/1_gradient_divergence_curl/labs/week02_div_curl.md
@@ -0,0 +1,1480 @@
+---
+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 2: Divergence and Curl
+
+:::{admonition} Computer lab
+:class: note
+
+The second of two labs on the operators of this chapter, following Lab 1 on series and the gradient. Parts 1 and 2 take the divergence, Part 3 the curl. Each task states a physical question, gives the steps, and ends with a self-check you can run. Plotting is supplied in the module `fwtools`, so that your effort goes into the physics rather than into rendering transparent isosurfaces.
+:::
+
+## Learning objectives
+
+By the end of this lab you should be able to:
+
+- **Distinguish diverging arrows from non-zero divergence.** Compute $\nabla\cdot\boldsymbol{v}$, justify the result by flux rather than by algebra, and identify the only radial flow that is incompressible.
+- **Use the divergence theorem as a measurement.** Verify $\oint_S\boldsymbol{v}\cdot\hat{\boldsymbol{n}}\,dS = \int_{\mathcal{D}} \nabla\cdot\boldsymbol{v}\,dV$ numerically, and account for what happens when the source shrinks to a point.
+- **Measure the curl as circulation per unit area.** Compute $\nabla\times\boldsymbol{v}$, separate rotation from the shape a streamline happens to make, and verify Stokes' theorem on a vortex with a finite core.
+
+---
+
+## Part 0 — Setup
+
+Run this once. It contains no physics: it fetches two packages the browser lacks, locates `fwtools`, and defines the Coulomb constant.
+
+```{code-cell} ipython3
+# No physics above the k_e = ... line near the bottom.
+import sys, pathlib
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.constants import epsilon_0
+
+# --- 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/1_gradient_divergence_curl/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"
+# ---------------------------------------------------------------------------
+
+k_e = 1.0 / (4.0 * np.pi * epsilon_0)   # Coulomb constant, 8.99e9 V*m/C
+Q = 1e-9                                # 1 nC test charge
+
+print(f"epsilon_0 = {epsilon_0:.4e} F/m")
+print(f"k_e       = {k_e:.4e} V*m/C")
+```
+
+### Carried over from Lab 1
+
+The same cube of sample points, and the fields Lab 1 built on it. Nothing here is an
+exercise: these are Lab 1's answers, given so that this notebook runs on its own.
+
+```{code-cell} ipython3
+n, L = 61, 2.0                          # odd n, so the origin is a sample point
+axis = np.linspace(-L, L, n)            # one axis, shared by x, y and z
+X, Y, Z = np.meshgrid(axis, axis, axis, indexing="ij")
+dx = dy = dz = axis[1] - axis[0]
+
+c = n // 2                              # index of the origin
+# Mask reused by the self-checks. `interior` drops the two outermost cells,
+# so comparisons exclude the six faces of the box, where sampling is worst
+# and np.gradient has only one-sided neighbours.
+interior = np.zeros(X.shape, dtype=bool)
+interior[2:-2, 2:-2, 2:-2] = True
+
+print(f"grid shape {X.shape},  spacing {dx:.4f} m,  {X.size:,} sample points")
+print(f"X[i,j,k] = x[i]   ->   X[-1, 0, 0] = {X[-1, 0, 0]:.1f} m")
+```
+
+```{code-cell} ipython3
+# Lab 1, Tasks 3, 4, 6 and 7 -- given here, not set again.
+def distance_to(X, Y, Z, x0=0.0, y0=0.0, z0=0.0):
+    return np.sqrt((X - x0)**2 + (Y - y0)**2 + (Z - z0)**2)
+
+r  = distance_to(X, Y, Z)                 # distance from the origin
+rs = np.maximum(r, 1e-12)                 # 0/0 at the source is not a lesson
+rhx, rhy, rhz = X / rs, Y / rs, Z / rs    # the outward unit radial vector
+
+r_masked = np.where(r < 0.25, np.nan, r)  # the singularity kept off the grid
+V = k_e * Q / r_masked                    # potential of the 1 nC point charge
+_dVx, _dVy, _dVz = np.gradient(V, dx, dy, dz)
+Ex, Ey, Ez = -_dVx, -_dVy, -_dVz          # E = -grad V
+
+print(f"grid {X.shape},  spacing {dx:.4f} m")
+print(f"|E| at (1,0,0) = {np.sqrt(Ex**2+Ey**2+Ez**2)[-1-15, c, c]:.3f} V/m")
+```
+
+---
+
+## Part 1 — Divergence
+
+The gradient takes a scalar and returns a vector. The divergence takes a vector field and returns a scalar:
+
+$$ \nabla\cdot\boldsymbol{A} \;=\; \lim_{\Delta V \to 0}\frac{1}{\Delta V}\oint_S \boldsymbol{A}\cdot\hat{\boldsymbol{n}}\,dS \;=\; \frac{\partial A_x}{\partial x} + \frac{\partial A_y}{\partial y} + \frac{\partial A_z}{\partial z} $$
+
+Read the definition on the left rather than the formula on the right: **treat $\boldsymbol{A}$ as a fluid velocity**, place a small box anywhere, and measure the net outflow through its walls per unit volume.
+
+| $\nabla\cdot\boldsymbol{A}$ | Name | Picture |
+| :---: | :--- | :--- |
+| $> 0$ | **source** | a tap: more leaves than arrives |
+| $< 0$ | **sink** | a drain: more arrives than leaves |
+| $= 0$ | **solenoidal** | whatever flows in, flows out |
+
+### Task 1 — the operator, and its independence of the origin
+
+The operator is three lines. One derivative along one axis per component: `np.gradient(Ax, dx, axis=0)` returns $\partial A_x/\partial x$ and nothing else, whereas asking for all three and discarding two costs three times the memory. The cross terms are not part of a divergence.
+
+**The question is the one raised by the definition.** Flux per unit volume is measured around a point, so does the result depend on which point is called the origin? Take the outward flow $\boldsymbol{A} = \boldsymbol{r}$, whose divergence follows on paper as $1+1+1 = 3$, then shift the whole field so that it streams out of $(0.8, -0.4, 0.3)$. Predict the divergence before computing it.
+
+```{code-cell} ipython3
+# Task 1
+def divergence(Ax, Ay, Az, dx, dy, dz):
+    return (np.gradient(Ax, dx, axis=0)
+            + ___
+            + ___)
+
+# The same outward flow, seen from somewhere else.
+x0, y0, z0 = 0.8, -0.4, 0.3
+Sx, Sy, Sz = ___                          # the field r - r0, as three arrays
+div_shifted = ___                         # its divergence
+
+# --- self-check (leave this alone) ---
+fw.check_close("div of the position vector = 3",
+               divergence(X, Y, Z, dx, dy, dz), 3.0, rtol=1e-6)
+fw.check_close("...and 3 again when the source is moved",
+               div_shifted, 3.0, rtol=1e-6)
+fw.check("the shifted field really is different from the original",
+         not np.allclose(Sx, X))
+```
+
+:::{admonition} Solution — Task 1
+:class: dropdown
+
+```python
+def divergence(Ax, Ay, Az, dx, dy, dz):
+    return (np.gradient(Ax, dx, axis=0)
+            + np.gradient(Ay, dy, axis=1)
+            + np.gradient(Az, dz, axis=2))
+
+Sx, Sy, Sz = X - x0, Y - y0, Z - z0
+div_shifted = divergence(Sx, Sy, Sz, dx, dy, dz)
+```
+:::
+
+:::{admonition} Why the answer had to be 3 either way
+:class: tip
+
+Moving the source changed every arrow in the box and changed the divergence nowhere. Differentiation removes the constant: $\partial(x - x_0)/\partial x = 1$ for any $x_0$.
+
+The divergence is a **local** quantity: it is built from a limit taken around one point, so it depends on the field in a shrinking neighbourhood of that point and not on where the axes were placed. Every operator in this course has that property, and it is what makes $\nabla\cdot\boldsymbol{E} = \rho_v/\varepsilon_0$ a statement about places rather than about coordinate systems.
+:::
+
+### Task 2 — the only incompressible radial flow
+
+Water of constant density flows outward from a source at the origin. Away from that source nothing is created or destroyed, so the flow is **incompressible**:
+
+$$ \nabla\cdot\boldsymbol{v} = 0 \qquad \text{for } r \neq 0. $$
+
+Constant density and a point source force the flow to be radial, $\boldsymbol{v} = f(r)\,\boldsymbol{r}$, and incompressibility then pins $f$ down completely:
+
+$$ \nabla\cdot\boldsymbol{v} = 3f(r) + r\frac{df}{dr} = 0 \qquad\Longrightarrow\qquad f(r) = \frac{A}{r^{3}}. $$
+
+Rather than assume this, test four candidates and let the divergence select.
+
+One decision comes first, because the four candidates are not the same size. Over the test band their divergences span a factor of a thousand, and the raw numbers cannot be ranked against each other: $f = 1/r^{4}$ returns a *smaller* $\lvert\nabla\cdot\boldsymbol{v}\rvert$ than $f = \text{const}$ does, and neither field is divergence-free. "Is 0.36 small?" has no answer until it is small compared with something.
+
+The something is the size a derivative of that same field would have if nothing cancelled. A derivative is a change in $\boldsymbol{v}$ divided by the distance over which it changes, and for a radial field $f(r)\boldsymbol{r}$ the only distance available is $r$ itself. That makes $\lvert\boldsymbol{v}\rvert/r$ the yardstick, and
+
+$$ \frac{\lvert\nabla\cdot\boldsymbol{v}\rvert}{\lvert\boldsymbol{v}\rvert/r} $$
+
+a pure number, the same for a trickle and a torrent: **1 means the three terms of the divergence did not cancel at all, and 0 means they cancelled completely.** The cell prints the raw divergence beside the ratio, so you can see for yourself why the raw column is unusable.
+
+One entry is known before the code runs. For $f = \text{const}$ the field is $\boldsymbol{v} = \boldsymbol{r}$, so $\lvert\boldsymbol{v}\rvert/r = 1$ and the ratio is nothing but $\nabla\cdot\boldsymbol{r} = 3$, which Task 1 measured. That row is the check that the statistic is being formed correctly, and it is why the last self-check looks for 300%.
+
+```{code-cell} ipython3
+# The statistic, restated: |div v| / (|v|/r), median over the test band.
+# The band avoids the source, where the field is singular, and the outer
+# corners of the box, where np.gradient runs out of neighbours.
+
+r_safe = np.where(r < 0.3, np.nan, r)
+band_i = interior & (r > 0.6) & (r < 1.6)
+
+# Task 2 -- three blanks, inside the loop.
+results, raws = {}, {}
+print(f"  {'f(r)':>7}  {'|div v|':>11}  {'|v|/r':>8}  {'ratio':>9}")
+for name, f_r in [("const", np.ones_like(r_safe)),
+                  ("1/r^2", 1/r_safe**2),
+                  ("1/r^3", 1/r_safe**3),
+                  ("1/r^4", 1/r_safe**4)]:
+    vx, vy, vz = ___                      # v = f(r) * (X, Y, Z): three arrays
+    dv = ___                              # its divergence (nan_to_num each part)
+    scale = ___                           # |v|/r AT THE BAND POINTS -- index it
+                                          # with [band_i], so it comes out 1-D
+                                          # and the same length as dv[band_i]
+
+    # --- given ---
+    raws[name] = np.nanmedian(np.abs(dv[band_i]))
+    results[name] = np.nanmedian(np.abs(dv[band_i]) / scale)
+    print(f"  {name:>7}  {raws[name]:11.4f}  {np.nanmedian(scale):8.4f}"
+          f"  {results[name]:9.2%}")
+
+# --- self-check (leave this alone) ---
+fw.check(f"the raw column cannot rank these: 1/r^4 gives a smaller |div v| "
+         f"({raws['1/r^4']:.3f}) than f = const ({raws['const']:.3f}), and "
+         f"neither is divergence-free", raws["1/r^4"] < raws["const"])
+fw.check(f"scale is one value per band point ({np.shape(scale)} vs "
+         f"{np.shape(dv[band_i])})", np.shape(scale) == np.shape(dv[band_i]),
+         "index it with [band_i] -- a whole-grid array or a single median "
+         "both change the statistic being reported")
+fw.check(f"1/r^3 is the divergence-free one ({results['1/r^3']:.2%})",
+         results["1/r^3"] < 0.05)
+fw.check("...and the other three are not",
+         min(results[k] for k in ("const", "1/r^2", "1/r^4")) > 0.5)
+fw.check(f"f = const reproduces Task 1's div(r) = 3 ({results['const']:.2%})",
+         np.isclose(results["const"], 3.0, rtol=1e-3))
+```
+
+:::{admonition} Solution — Task 2
+:class: dropdown
+
+```python
+    vx, vy, vz = f_r*X, f_r*Y, f_r*Z
+    dv = divergence(*(np.nan_to_num(q) for q in (vx, vy, vz)), dx, dy, dz)
+    scale = (np.sqrt(vx**2 + vy**2 + vz**2) / r_safe)[band_i]
+```
+:::
+
+:::{admonition} Where the inverse-square law comes from
+:class: important
+
+Two candidates give almost exactly 100%, meaning their three divergence terms did not cancel at all, and the anchor gives its predicted 300%. One gives 0.66%. Only $f = A/r^{3}$ survives, as the algebra predicts, and the raw column beside it would have told you none of this. The surviving case rewrites as
+
+$$ \boldsymbol{v} = \frac{A}{r^{3}}\boldsymbol{r} = \frac{A}{r^{2}}\,\hat{\boldsymbol{r}}. $$
+
+**This is the same $1/r^{2}$ used since Lab 1's Task 6.** Here it was not assumed and no charge was mentioned; it follows from conservation away from the source together with the three-dimensionality of space. The surface of a sphere grows as $r^{2}$, so a fixed flux crossing it must thin as $1/r^{2}$.
+
+Coulomb's law, Newtonian gravity and this flow share an exponent for that one geometric reason.
+:::
+
+### Task 2, continued — a field with no source anywhere
+
+Note the restriction on that result: $\nabla\cdot\boldsymbol{v} = 0$ **for $r \neq 0$**. The origin must be excluded, because that is where the water is injected; a closed surface around it would find the tap.
+
+The next field admits no such exception. To first order the Earth's magnetic field is a **dipole**: a north and a south pole so close together that they coincide. With dipole moment $\boldsymbol{m}$,
+
+$$ \boldsymbol{B} = \frac{3\boldsymbol{r}\,(\boldsymbol{r}\cdot\boldsymbol{m}) - r^{2}\boldsymbol{m}}{r^{5}}. $$
+
+Take $\boldsymbol{m} = \hat{\boldsymbol{z}}$ on the cube set up above, where $z$ points up, and measure the divergence with the same function.
+
+The Earth's own moment points roughly geographic south, which is why the magnetic pole in the Arctic is magnetically a **south** pole and attracts the north end of a compass needle. Reversing $\boldsymbol{m}$ reverses every arrow below and leaves $\nabla\cdot\boldsymbol{B}$ unchanged.
+
+```{code-cell} ipython3
+# Task 2, continued -- fill in the three components.
+#   With m = z-hat, the dot product r . m is simply Z.
+#   Careful with the second term: it appears only in the z-component.
+
+r_dot_m = Z
+Bx = ___
+By = ___
+Bz = ___
+
+div_B = divergence(np.nan_to_num(Bx), np.nan_to_num(By), np.nan_to_num(Bz), dx, dy, dz)
+
+# --- given: the same scale-free measure as above ---
+B_mag = np.sqrt(Bx**2 + By**2 + Bz**2)
+print(f"  dipole B : median |div B| / (|B|/r) = "
+      f"{np.nanmedian(np.abs(div_B[band_i]) / (B_mag/r_safe)[band_i]):8.2%}")
+
+# --- self-check (leave this alone) ---
+fw.check("B is divergence-free",
+         np.nanmedian(np.abs(div_B[band_i]) / (B_mag/r_safe)[band_i]) < 0.05)
+fw.check("B is not simply radial (it has a north and a south)",
+         np.nanmin((Bx*X + By*Y + Bz*Z)[band_i]) < 0)
+```
+
+:::{admonition} Solution — Task 2, continued
+:class: dropdown
+
+```python
+Bx = 3*X*r_dot_m / r_safe**5
+By = 3*Y*r_dot_m / r_safe**5
+Bz = (3*Z*r_dot_m - r_safe**2) / r_safe**5
+```
+:::
+
+:::{admonition} No magnetic monopoles
+:class: important
+
+Both fields are divergence-free over the region measured, but the two statements differ.
+
+The flow required an exclusion: $\nabla\cdot\boldsymbol{v} = 0$ away from the origin, because the origin is a tap. The dipole requires none, and $\nabla\cdot\boldsymbol{B} = 0$ holds **everywhere in space, including at the source**. No point can be excluded to reveal a magnet leaking field the way the tap leaks water. This is one of Maxwell's equations: magnetic monopoles do not exist, and field lines of $\boldsymbol{B}$ never begin or end but close on themselves.
+
+Two remarks on the numbers. Both cells report the same scale-free measure, so the results are directly comparable. The dipole's 1.8% is worse than the radial flow's 0.66%, not because the physics is less secure but because $\boldsymbol{B}$ falls off as $1/r^{3}$ rather than $1/r^{2}$, leaving a centred difference more curvature to miss. Part 2 tests the same claim far below 2% by putting a closed surface around the dipole instead of differentiating it.
+
+The second check is also informative: $\boldsymbol{B}\cdot\boldsymbol{r}$ is negative somewhere, whereas the outward flow of Task 2 is never negative. The dipole points inward over part of space; it returns. That is the numerical signature of a field closing on itself.
+:::
+
+### Task 3 — three flows
+
+Three velocity fields. For each one: **sketch it, predict the sign of the divergence, then measure.** Record the predictions first; the task is about the gap between intuition and the result.
+
+| | Field $\boldsymbol{A}$ | What it looks like |
+| :---: | :--- | :--- |
+| **(a)** | $x\,\hat{\boldsymbol{x}} + y\,\hat{\boldsymbol{y}} + z\,\hat{\boldsymbol{z}}$ | outward flow in all directions |
+| **(b)** | $-y\,\hat{\boldsymbol{x}} + x\,\hat{\boldsymbol{y}}$ | fluid rotating about the $z$-axis |
+| **(c)** | $x\,\hat{\boldsymbol{x}} - y\,\hat{\boldsymbol{y}}$ | stretching along $x$, squeezing along $y$ |
+
+```{code-cell} ipython3
+# Record the predictions BEFORE running the next cell: +1 for a source,
+# -1 for a sink, 0 for solenoidal. The next cell scores them.
+predictions = {"a": ___, "b": ___, "c": ___}
+```
+
+```{code-cell} ipython3
+# Task 3 -- six blanks: three fields, three divergences.
+zero = np.zeros_like(X)
+Aa = ___                                  # (a) outward flow, as a triple
+Ab = ___                                  # (b) rotation about z
+Ac = ___                                  # (c) stretch in x, squeeze in y
+
+div_a = ___
+div_b = ___
+div_c = ___
+
+# --- given: the three side by side, one shared scale, one colorbar ---
+for name, d in [("(a) outward flow", div_a), ("(b) rotation", div_b),
+                ("(c) straining flow", div_c)]:
+    print(f"{name:20s} div = {d.mean():+.3f}")
+
+fig, axes = plt.subplots(1, 3, figsize=(16, 4.6))
+for ax_, (name, A, d) in zip(axes, [("(a) outward flow", Aa, div_a),
+                                    ("(b) rotation", Ab, div_b),
+                                    ("(c) straining flow", Ac, div_c)]):
+    fw.show_field_slice(X, Y, Z, *A[:2], background=d, ax=ax_, density=1.1,
+                        vmin=-3, vmax=3, colorbar=(ax_ is axes[-1]),
+                        label=r"$\nabla\cdot\mathbf{A}$  [s$^{-1}$]", title=name)
+plt.tight_layout()
+plt.show()
+# Examine (b) and (c) before reading the note below: both come out a uniform
+# zero, for entirely different reasons.
+
+# --- self-check (leave this alone) ---
+# (a) has a non-zero answer, so a relative test works. (b) and (c) are exactly
+# zero, and nothing can be measured relative to zero, so they get an absolute
+# tolerance instead.
+fw.check_close("(a) div = 3", div_a, 3.0, rtol=1e-6)
+fw.check_abs("(b) div = 0 (rotation)", div_b, atol=1e-9)
+fw.check_abs("(c) div = 0 (straining flow)", div_c, atol=1e-9)
+
+for key, measured in (("a", div_a), ("b", div_b), ("c", div_c)):
+    sign = int(np.sign(np.round(measured.mean(), 6)))
+    verdict = "as predicted" if predictions[key] == sign else "NOT what you predicted"
+    print(f"  ({key}) you said {predictions[key]:+d}, measured {sign:+d}  --  {verdict}")
+```
+
+:::{admonition} Solution — Task 3
+:class: dropdown
+
+```python
+Aa = (X, Y, Z)
+Ab = (-Y, X, zero)
+Ac = (X, -Y, zero)
+
+div_a = divergence(*Aa, dx, dy, dz)
+div_b = divergence(*Ab, dx, dy, dz)
+div_c = divergence(*Ac, dx, dy, dz)
+```
+:::
+
+:::{admonition} Field (c) is the trap
+:class: warning
+
+Along the $x$-axis, field (c) flows outward and resembles a source. It is not:
+
+$$ \nabla\cdot\boldsymbol{A} = \frac{\partial}{\partial x}(x) + \frac{\partial}{\partial y}(-y) = 1 - 1 = 0 $$
+
+Place a box at the origin: fluid leaves through the left and right walls and enters through the top and bottom at exactly the same rate. The parcel changes **shape**, not **volume**.
+
+Diverging arrows are not divergence. Outflow in one direction can be cancelled exactly by inflow in another. Task 5 puts a closed surface around this field and measures the cancellation directly.
+:::
+
+### Task 4 — the divergence as a charge detector
+
+Gauss's law, for a field in vacuum, says
+
+$$ \nabla\cdot\boldsymbol{E} = \frac{\rho_v}{\varepsilon_0} $$
+
+which is a strong claim: **the divergence of $\boldsymbol{E}$ at a point gives the charge density at that point and nothing else.** Where there is no charge, $\boldsymbol{E}$ is solenoidal, however widely its arrows spread.
+
+Test that pointwise, on a source a grid can hold. A point charge cannot serve: it has infinite density at one location. Take instead a charge **distributed over a finite blob**, which is what any real charged object is:
+
+$$ \rho_v(r) = \rho_{v0}\,e^{-r^{2}/a^{2}}, \qquad \rho_{v0} = 10^{-9}\ \text{C/m}^3, \qquad a = 0.5\ \text{m} $$
+
+Here $a$ is the **width of the blob**. In Lab 1's Task 9 the same letter denoted an electrode separation, the second symbol these two labs overload, after $\rho$. The code keeps them apart as `a` here and `a_sep` there; in algebra only the context distinguishes them.
+
+Integrating over a sphere of radius $r$ gives the charge it encloses:
+
+$$ Q_{\text{enc}}(r) = \int_0^{r}\!\rho_v\,4\pi r'^{2}\,dr' = 4\pi\rho_{v0}\left[\frac{a^{3}\sqrt{\pi}}{4}\operatorname{erf}\!\left(\frac{r}{a}\right) - \frac{a^{2}r}{2}e^{-r^{2}/a^{2}}\right] $$
+
+and Gauss's law, $E_r = Q_{\text{enc}}/4\pi\varepsilon_0r^{2}$, then gives the field, with the $4\pi$ cancelling:
+
+$$ E_r(r) = \frac{\rho_{v0}}{\varepsilon_0 r^{2}}\left[\frac{a^{3}\sqrt{\pi}}{4}\operatorname{erf}\!\left(\frac{r}{a}\right) - \frac{a^{2}r}{2}e^{-r^{2}/a^{2}}\right] $$
+
+One check: near the centre $Q_{\text{enc}}$ grows as $r^{3}$ while the surface grows as $r^{2}$, so $E_r \to \rho_{v0} r/3\varepsilon_0$, zero at the centre, rising linearly, and peaking at $r \approx a$.
+
+```{code-cell} ipython3
+from scipy.special import erf
+
+a, rho_v0 = 0.5, 1e-9
+
+# --- given: the charge density, and the field Gauss's law gives it ---
+# The two bracketed terms nearly cancel for r << a, so the closed form loses
+# accuracy below r ~ 1e-6 m. On this grid the only such sample is the origin,
+# where the r-hat components are zero in any case.
+rho_v = rho_v0 * np.exp(-r**2 / a**2)
+E_r = rho_v0 / (epsilon_0 * rs**2) * (
+    (a**3 * np.sqrt(np.pi) / 4) * erf(rs / a) - (a**2 * rs / 2) * np.exp(-rs**2 / a**2)
+)
+
+# Task 4 -- two blanks.
+#   E_r is a radial MAGNITUDE. Give it a direction, then differentiate.
+Ex_b, Ey_b, Ez_b = ___                    # components along (rhx, rhy, rhz)
+div_blob = ___                            # your Task 1 operator
+
+# --- given: the two pictures, forced onto one scale so they are comparable ---
+hi = float(np.nanmax(rho_v / epsilon_0))
+units = r"[V m$^{-2}$]"
+fig, axes = plt.subplots(1, 2, figsize=(12, 4.4))
+fw.show_scalar_slice(X, Y, Z, div_blob, ax=axes[0], cmap="magma", label=units,
+                     vmin=0, vmax=hi, title=r"measured $\nabla\cdot\mathbf{E}$")
+fw.show_scalar_slice(X, Y, Z, rho_v / epsilon_0, ax=axes[1], cmap="magma", label=units,
+                     vmin=0, vmax=hi, title=r"actual $\rho_v/\varepsilon_0$")
+plt.tight_layout()
+plt.show()
+
+print(f"peak of rho_v/eps0  : {np.nanmax(rho_v/epsilon_0):8.2f}")
+print(f"peak of measured div: {np.nanmax(div_blob):8.2f}")
+
+# --- self-check (leave this alone) ---
+peak = np.nanmax(rho_v / epsilon_0)
+_e = np.abs(div_blob[interior] - (rho_v / epsilon_0)[interior]) / peak
+cart_worst, cart_median = float(_e.max()), float(np.median(_e))
+fw.check(f"div E = rho_v/eps0 pointwise (worst {cart_worst:.2%} of peak)",
+         cart_worst < 0.05, "check the component construction Ex_b = E_r * rhx")
+```
+
+:::{admonition} Solution — Task 4
+:class: dropdown
+
+```python
+Ex_b, Ey_b, Ez_b = E_r * rhx, E_r * rhy, E_r * rhz
+div_blob = divergence(Ex_b, Ey_b, Ez_b, dx, dy, dz)
+```
+:::
+
+:::{admonition} What the two panels show
+:class: important
+
+The two panels show the same distribution. The location of the charge was never supplied to the code: a field was differentiated, and the charge distribution came back out.
+
+Note where the divergence vanishes: everywhere outside the blob, where the field is still large and still spreading. **Strong field, zero divergence**: the two quantities are unrelated.
+:::
+
+### The same operator, a different formula
+
+Everything so far used the Cartesian formula, because `np.gradient` differentiates along array axes. The divergence is flux per unit volume, a physical quantity that cannot depend on the choice of axes. Only the formula changes:
+
+| | Gradient $\nabla T$ | Divergence $\nabla\cdot\boldsymbol{A}$ |
+| :--- | :--- | :--- |
+| Cartesian $(x,y,z)$ | $\dfrac{\partial T}{\partial x}\hat{\boldsymbol{x}} + \dfrac{\partial T}{\partial y}\hat{\boldsymbol{y}} + \dfrac{\partial T}{\partial z}\hat{\boldsymbol{z}}$ | $\dfrac{\partial A_x}{\partial x} + \dfrac{\partial A_y}{\partial y} + \dfrac{\partial A_z}{\partial z}$ |
+| Cylindrical $(\varrho,\phi,z)$ | $\dfrac{\partial T}{\partial \varrho}\hat{\boldsymbol{\varrho}} + \dfrac{1}{\varrho}\dfrac{\partial T}{\partial \phi}\hat{\boldsymbol{\phi}} + \dfrac{\partial T}{\partial z}\hat{\boldsymbol{z}}$ | $\dfrac{1}{\varrho}\dfrac{\partial (\varrho v_\varrho)}{\partial \varrho} + \dfrac{1}{\varrho}\dfrac{\partial v_\phi}{\partial \phi} + \dfrac{\partial v_z}{\partial z}$ |
+| Spherical $(r,\phi,\theta)$ | $\dfrac{\partial T}{\partial r}\hat{\boldsymbol{r}} + \dfrac{1}{r}\dfrac{\partial T}{\partial \theta}\hat{\boldsymbol{\theta}} + \dfrac{1}{r\sin\theta}\dfrac{\partial T}{\partial \phi}\hat{\boldsymbol{\phi}}$ | $\dfrac{1}{r^{2}}\dfrac{\partial (r^{2}v_r)}{\partial r} + \dfrac{1}{r\sin\theta}\dfrac{\partial (v_\theta \sin\theta)}{\partial \theta} + \dfrac{1}{r\sin\theta}\dfrac{\partial v_\phi}{\partial \phi}$ |
+
+Cylindrical $\varrho=\sqrt{x^2+y^2}$ is the distance from the $z$-axis; spherical $r=\sqrt{x^2+y^2+z^2}$, used throughout this lab, is the distance from the origin. They are written differently precisely to keep them apart.
+
+One reading note: the spherical coordinates are named $(r,\phi,\theta)$, but the terms in the row above are listed as $r$, then $\theta$, then $\phi$, the order in which the scale factors $(1,\ r,\ r\sin\theta)$ are derived. The order of terms in a sum is immaterial.
+
+Both fields built so far are spherically symmetric, $\boldsymbol{E} = E_r(r)\,\hat{\boldsymbol{r}}$ with no $\theta$ or $\phi$ dependence, so two of the three spherical terms vanish and the divergence reduces to one ordinary derivative along one line:
+
+$$ \nabla\cdot\boldsymbol{E} \;=\; \frac{1}{r^{2}}\frac{d}{dr}\!\left(r^{2}E_r\right) $$
+
+```{code-cell} ipython3
+dr = 0.005
+r_line = np.arange(0.05, 2.0 + dr, dr)          # one radial line, not a cube
+
+# the same two fields as before, as functions of r alone
+E_R_blob = rho_v0 / (epsilon_0 * r_line**2) * (
+    (a**3 * np.sqrt(np.pi) / 4) * erf(r_line / a)
+    - (a**2 * r_line / 2) * np.exp(-r_line**2 / a**2))
+E_r_point = k_e * Q / r_line**2
+
+div_blob_sph = np.gradient(r_line**2 * E_R_blob, dr) / r_line**2
+div_point_sph = np.gradient(r_line**2 * E_r_point, dr) / r_line**2
+
+rho_v_line = rho_v0 * np.exp(-r_line**2 / a**2)
+
+# --- given: the radial profile Task 4 asserted but never drew ---
+Q_total = np.pi**1.5 * rho_v0 * a**3                  # all of the blob's charge
+plt.figure(figsize=(5.8, 4.2))
+plt.plot(r_line, E_R_blob, "k", lw=1.8, label="$E_r(r)$, exact")
+plt.plot(r_line, rho_v0 * r_line / (3 * epsilon_0), "C1--", lw=1.2,
+         label=r"small $r$:  $\rho_{v0}r/3\varepsilon_0$")
+plt.plot(r_line, Q_total / (4 * np.pi * epsilon_0 * r_line**2), "C2:", lw=1.4,
+         label=r"large $r$:  $Q/4\pi\varepsilon_0 r^2$")
+plt.axvline(a, color="C0", lw=1, alpha=0.6)
+plt.annotate("$r = a$", (a, 1.12 * E_R_blob.max()), color="C0", ha="left")
+plt.xlabel("$r$ [m]"); plt.ylabel(r"$E_r$ [V m$^{-1}$]")
+plt.ylim(0, 1.25 * E_R_blob.max()); plt.grid(alpha=0.3); plt.legend(fontsize=8)
+plt.title("The blob's field: linear inside, inverse-square outside")
+plt.show()
+
+err_sph = np.abs(div_blob_sph - rho_v_line / epsilon_0)[1:-1] / np.max(rho_v_line / epsilon_0)
+print(f"blob : {r_line.size} samples on a line vs {X.size:,} in the cube")
+print(f"       worst error {err_sph.max():.3%} of peak, median {np.median(err_sph):.4%}")
+print(f"       Cartesian, from Task 4: {cart_worst:.3%} and {cart_median:.4%}")
+print(f"point: r^2 E_r varies by {np.ptp(r_line**2 * E_r_point):.1e} over the whole line")
+print(f"       max |div E| = {np.abs(div_point_sph).max():.1e}  (round-off, not physics)")
+```
+
+:::{admonition} Why curvilinear coordinates are worth the trouble
+:class: important
+
+Same field, same operator, same answer, obtained from a few hundred samples on a line rather than a quarter of a million in a cube, and several times more accurately.
+
+For the point charge the gain is certainty rather than accuracy. $r^{2}E_r = Q/4\pi\varepsilon_0$ is a **constant**, so its derivative is analytically zero for every $r>0$: not 1% of something, but zero, in one line of algebra. What the cell prints is the precision with which double arithmetic subtracts two equal numbers, around $10^{-13}$, or exactly $0$ when the cancellation is exact. Changing `dr` moves that last digit; it does not move the algebra. Cartesian coordinates could establish only that the divergence is small.
+
+Matching the coordinates to the symmetry of the source replaces three noisy numerical derivatives with one line of algebra. That is the purpose of the second and third rows of the table.
+:::
+
+---
+
+## Part 2 — Flux, and the divergence theorem
+
+Part 1 used the differential form of Gauss's law, which compares two numbers at one point. The integral form relates a volume to the surface enclosing it:
+
+$$ \oint_S \boldsymbol{E}\cdot\hat{\boldsymbol{n}}\,dS \;=\; \int_{\mathcal{D}} \nabla\cdot\boldsymbol{E}\;dV \;=\; \frac{Q_{\text{enc}}}{\varepsilon_0} $$
+
+with $S$ the closed surface, $\hat{\boldsymbol{n}}$ its outward unit normal, and $\mathcal{D}$ the volume it encloses.
+
+The first equality is the **divergence theorem**, pure vector calculus, valid for any well-behaved field. The second is the physics. Together they state that measuring $\boldsymbol{E}$ on a closed surface gives the charge inside, and nothing about its arrangement or about any charge outside.
+
+Take $S$ to be a cube of half-width $h$ centred on the origin, faces on grid planes. On the $+x$ face the outward normal is $+\hat{\boldsymbol{x}}$, so it contributes $\int\!\!\int E_x\,dy\,dz$; on the $-x$ face the normal is $-\hat{\boldsymbol{x}}$ and the same integral enters negatively. Six faces, three pairs.
+
+### Task 5 — close the surface
+
+```{code-cell} ipython3
+# `fw.area_integral(F2, da, db)` integrates a 2-D array over the face it
+# spans; `fw.volume_integral(F3, dx, dy, dz)` does the same over a box.
+# `fw.box_indices(X, h)` gives the index range of the cube |x|,|y|,|z| <= h.
+#
+# The x pair is written for you; the pattern is one row per axis:
+#
+#     face pair   outward samples   inward samples    spacings
+#     x           Ax[i1, s, s]      Ax[i0, s, s]      dy, dz
+#     y           Ay[s, i1, s]      Ay[s, i0, s]      dx, dz
+#     z           Az[s, s, i1]      Az[s, s, i0]      dx, dy
+#
+# The axis you pin to i0/i1 is the axis whose spacing you leave out.
+# NOTE: this closes over X, dx, dy, dz from the cell above, so it is tied to
+# this grid and is not a general-purpose function.
+
+def closed_box_flux(Ax, Ay, Az, half_width):
+    """Net outward flux through the cube |x|,|y|,|z| <= half_width."""
+    i0, i1 = fw.box_indices(X, half_width)
+    s = slice(i0, i1 + 1)
+    flux_x = (fw.area_integral(Ax[i1, s, s], dy, dz)
+              - fw.area_integral(Ax[i0, s, s], dy, dz))
+    flux_y = ___
+    flux_z = ___
+    return flux_x + flux_y + flux_z
+
+
+# --- given: three routes to the same number, and the Task 3 arbiter ---
+# Finish closed_box_flux above; everything below is written for you.
+flux_1m = closed_box_flux(Ex_b, Ey_b, Ez_b, 1.0)
+
+print(f"{'h [m]':>6} {'surface':>12} {'volume':>12} {'Q_enc/eps0':>12}")
+for h in (0.6, 1.0, 1.4):
+    i0, i1 = fw.box_indices(X, h)
+    s_ = slice(i0, i1 + 1)
+    surf = closed_box_flux(Ex_b, Ey_b, Ez_b, h)
+    vol = fw.volume_integral(div_blob[s_, s_, s_], dx, dy, dz)
+    qenc = fw.volume_integral(rho_v[s_, s_, s_], dx, dy, dz) / epsilon_0
+    print(f"{h:6.1f} {surf:12.3f} {vol:12.3f} {qenc:12.3f}")
+
+# Task 3 settled by measurement rather than by argument. Both (b) and (c)
+# appear to throw fluid outwards somewhere; a closed surface is the arbiter,
+# and it differentiates nothing.
+print(f"\nflux of (b), the rotation       : {closed_box_flux(-Y, X, zero, 1.0):+.2e}")
+print(f"flux of (c), the straining flow : {closed_box_flux(X, -Y, zero, 1.0):+.2e}")
+
+# --- self-check (leave this alone) ---
+i0, i1 = fw.box_indices(X, 1.0)
+s = slice(i0, i1 + 1)
+fw.check_scalar("closed-surface flux = Q_enc/eps0", flux_1m,
+                fw.volume_integral(rho_v[s, s, s], dx, dy, dz) / epsilon_0,
+                rtol=0.01, unit=" V*m")
+fw.check_scalar("divergence theorem: surface = volume", flux_1m,
+                fw.volume_integral(div_blob[s, s, s], dx, dy, dz),
+                rtol=0.01, unit=" V*m")
+```
+
+:::{admonition} Solution — Task 5
+:class: dropdown
+
+```python
+    flux_y = (fw.area_integral(Ay[s, i1, s], dx, dz)
+              - fw.area_integral(Ay[s, i0, s], dx, dz))
+    flux_z = (fw.area_integral(Az[s, s, i1], dx, dy)
+              - fw.area_integral(Az[s, s, i0], dx, dy))
+    return flux_x + flux_y + flux_z
+```
+:::
+
+:::{admonition} Three routes, one number
+:class: important
+
+Three independent calculations. The first never examines the interior of the box, the second never examines the surface, and the third never examines the field. They agree to a fraction of a percent.
+
+The result grows with $h$ and then stops: once the cube holds nearly all the charge, enlarging it adds surface but no charge. Charge outside a closed surface contributes exactly nothing, because the field lines it sends in through one wall leave through another.
+:::
+
+### Shrinking the source to a point
+
+Run the same surface integral on the point-charge field of Lab 1's Task 7, whose divergence could not be measured at the origin because the singularity had to be masked.
+
+Rearranged, Gauss's law turns the flux into a **charge meter**: $Q_{\text{enc}} = \varepsilon_0 \oint_S \boldsymbol{E}\cdot\hat{\boldsymbol{n}}\,dS$. Weigh the charge inside each box in coulombs and compare it with the 1 nC placed there.
+
+```{code-cell} ipython3
+print("box half-width      charge it finds")
+for h in (0.6, 1.0, 1.4):
+    Q_found = epsilon_0 * closed_box_flux(Ex, Ey, Ez, h)
+    print(f"   {h:.1f} m           {Q_found * 1e12:8.2f} pC")
+print(f"\n   actually there   {Q * 1e12:8.2f} pC")
+
+# The shell between the 0.6 m and 1.4 m boxes holds no charge. Weigh it: what
+# enters the small box must leave the large one, so the difference of the two
+# fluxes is the charge in between.
+Q_shell = epsilon_0 * (closed_box_flux(Ex, Ey, Ez, 1.4)
+                       - closed_box_flux(Ex, Ey, Ez, 0.6))
+print(f"\ncharge in the shell between them: {Q_shell * 1e12:+.2f} pC "
+      f"({abs(Q_shell) / Q:.2%} of the charge at the centre)")
+```
+
+:::{admonition} Where did the charge go?
+:class: important
+
+Every box weighs the same 1 nC to a fraction of a percent, and the shell between two of them weighs nothing. All the charge lies in the only region common to every box: the origin.
+
+The whole source therefore sits at one point, where $\nabla\cdot\boldsymbol{E}$ is not a large number but undefined: $\rho_v$ has become a **Dirac delta**, zero everywhere, infinite at one point, with finite integral $Q$. The integral form survives exactly where the differential form fails.
+
+The same statement for magnetism carries no source term at all:
+
+$$ \nabla\cdot\boldsymbol{B} = 0 \qquad\Longleftrightarrow\qquad \oint_S \boldsymbol{B}\cdot\hat{\boldsymbol{n}}\,dS = 0 \ \ \text{for every closed } S $$
+
+The measurement returns zero around any closed surface anywhere: there are no magnetic monopoles, and field lines of $\boldsymbol{B}$ never begin or end.
+:::
+
+### The dipole, exactly
+
+Task 2 measured $\nabla\cdot\boldsymbol{B} = 0$ for the Earth's dipole and returned 1.8%, which is grid error rather than physics. The same claim can be tested without differentiating: put a closed surface around the dipole and weigh what crosses it.
+
+One warning before reading the numbers. A box centred on the origin is too easy a test for this dipole: with $\boldsymbol{m} = \hat{\boldsymbol{z}}$, $B_x$ and $B_y$ are odd in $z$ and $B_z$ is even, so on a $z$-symmetric box the faces cancel in pairs before any physics enters. An off-centre box is the honest test, and the cell below runs both.
+
+```{code-cell} ipython3
+r_dot_m = Z
+Bx = 3*X*r_dot_m / r_safe**5
+By = 3*Y*r_dot_m / r_safe**5
+Bz = (3*Z*r_dot_m - r_safe**2) / r_safe**5
+
+B = tuple(np.nan_to_num(q) for q in (Bx, By, Bz))
+v = tuple(np.nan_to_num(q / r_safe**3) for q in (X, Y, Z))
+
+# --- given: the same surface integral over any grid-aligned box, centred
+#     on the origin or not. Same six faces, same three pairs as Task 5.
+def box_flux(A, x0, x1, y0, y1, z0, z1):
+    i = [int(np.argmin(np.abs(axis - q))) for q in (x0, x1, y0, y1, z0, z1)]
+    sx, sy, sz = slice(i[0], i[1]+1), slice(i[2], i[3]+1), slice(i[4], i[5]+1)
+    return (fw.area_integral(A[0][i[1], sy, sz], dy, dz) - fw.area_integral(A[0][i[0], sy, sz], dy, dz)
+          + fw.area_integral(A[1][sx, i[3], sz], dx, dz) - fw.area_integral(A[1][sx, i[2], sz], dx, dz)
+          + fw.area_integral(A[2][sx, sy, i[5]], dx, dy) - fw.area_integral(A[2][sx, sy, i[4]], dx, dy))
+
+boxes = [("centred, h = 0.6", (-0.6, 0.6, -0.6, 0.6, -0.6, 0.6)),
+         ("centred, h = 1.0", (-1.0, 1.0, -1.0, 1.0, -1.0, 1.0)),
+         ("centred, h = 1.4", (-1.4, 1.4, -1.4, 1.4, -1.4, 1.4)),
+         ("lopsided in z   ", (-1.0, 1.0, -1.0, 1.0, -0.6, 1.0)),
+         ("lopsided in x, z", (-0.6, 1.0, -1.0, 1.0, -0.6, 1.0))]
+
+print("  box                 flux of B      flux of the radial flow v")
+for name, lim in boxes:
+    print(f"  {name}   {box_flux(B, *lim):+11.2e}   {box_flux(v, *lim):+12.4f}")
+print(f"\n  4*pi = {4*np.pi:.4f};  the integrator misses it by "
+      f"{4*np.pi - box_flux(v, -1, 1, -1, 1, -1, 1):.1e} on the flow")
+
+# --- self-check (leave this alone) ---
+fw.check("the dipole encloses nothing -- even in a box that is not centred on it",
+         max(abs(box_flux(B, *lim)) for _, lim in boxes) < 1e-2)
+fw.check("...and the same integrator does find the tap in the radial flow",
+         abs(box_flux(v, -1, 1, -1, 1, -1, 1) - 4*np.pi) < 0.01 * 4*np.pi)
+```
+
+:::{admonition} Two kinds of "divergence-free"
+:class: important
+
+The radial flow returns $4\pi$ through every surface, whatever its size: there is a tap at the origin, and every box finds the same one, as every box found the same 1 nC above.
+
+The dipole returns **nothing** through any of them. On the three centred boxes the result is zero to machine precision, but the warning above applies: those boxes cancel the field against itself by symmetry and could return nothing else. The off-centre boxes are the measurement that counts, and they return $-2.5\times10^{-3}$ and $-2.2\times10^{-4}$. Compare the adjacent column: the same integrator, on the same grid, misses $4\pi$ by $6.8\times10^{-3}$ on the radial flow. **The flux of the dipole is zero to better than the accuracy this method achieves on anything.**
+
+The two divergence-free fields are therefore different statements. The flow has a tap that can be located by shrinking a surface onto it; the dipole has nothing to locate, at any size or placement of the surface. This is $\nabla\cdot\boldsymbol{B} = 0$ in the form that admits no exception, and it is why the integral form is worth constructing: it settles the question at the source, where the differential form had to be masked.
+:::
+
+### Where do the 1% errors come from?
+
+Every derivative on this page is a centred difference, accurate to $O(\Delta x^{2})$: halving the spacing should reduce the error by four. Confirm it. The study is a single loop.
+
+```{code-cell} ipython3
+print(f"{'n':>4} {'dx [m]':>8} {'worst error':>12} {'ratio':>7}")
+prev = None
+for n_test in (21, 31, 41, 61):
+    ax_t = np.linspace(-L, L, n_test)
+    h_t = ax_t[1] - ax_t[0]
+    Xt, Yt, Zt = np.meshgrid(ax_t, ax_t, ax_t, indexing="ij")
+    rt = np.sqrt(Xt**2 + Yt**2 + Zt**2)
+    Rst = np.maximum(rt, 1e-12)
+    rho_t = rho_v0 * np.exp(-rt**2 / a**2)
+    E_Rt = rho_v0 / (epsilon_0 * Rst**2) * (
+        (a**3 * np.sqrt(np.pi) / 4) * erf(Rst / a)
+        - (a**2 * Rst / 2) * np.exp(-Rst**2 / a**2))
+    dv = divergence(E_Rt * Xt / Rst, E_Rt * Yt / Rst, E_Rt * Zt / Rst, h_t, h_t, h_t)
+    inner = np.zeros(Xt.shape, bool)
+    inner[2:-2, 2:-2, 2:-2] = True
+    e = np.nanmax(np.abs(dv[inner] - (rho_t / epsilon_0)[inner])) / np.nanmax(rho_t / epsilon_0)
+    ratio = "-" if prev is None else f"{prev / e:.2f}"
+    print(f"{n_test:>4} {h_t:>8.4f} {e:>11.2%} {ratio:>7}")
+    prev = e
+```
+
+:::{admonition} Second order, by measurement
+:class: important
+
+Compare each ratio with the square of the spacing ratio: $1.5^2 = 2.25$ from $n=21$ to $31$, $1.33^2 = 1.78$ from $31$ to $41$, and $1.5^2 = 2.25$ from $41$ to $61$.
+
+The 1.06% in Task 4 is therefore not noise to be tolerated but a predictable quantity that can be reduced at a known cost, and the choice of $n = 61$ in Lab 1 can now be audited rather than assumed.
+:::
+
+---
+
+## Part 3 — Curl
+
+Field **(b)**, the rotation, has zero divergence everywhere, yet it plainly circulates. The divergence cannot detect circulation. The operator that can is the **curl**, which takes a vector field and returns another vector field:
+
+$$ \nabla\times\boldsymbol{v} \;=\; \hat{\boldsymbol{x}}\left(\partial_y v_z - \partial_z v_y\right) + \hat{\boldsymbol{y}}\left(\partial_z v_x - \partial_x v_z\right) + \hat{\boldsymbol{z}}\left(\partial_x v_y - \partial_y v_x\right) $$
+
+There is no determinant to memorise. Each component pairs an even permutation of $(x,y,z)$ against an odd one, in three places at once: the direction, the differentiation, and the vector component. The $\hat{\boldsymbol{x}}$ term takes $(x,y,z)$ minus $(x,z,y)$, and the other two follow by advancing every letter one step, $x\to y\to z\to x$.
+
+The interpretation comes from a circulation integral. Take a small rectangle of side $dy$ by $dz$ around a point, walk its four edges once round, and add up the component of $\boldsymbol{v}$ along the direction of travel. Expanding each edge to first order in a Taylor series leaves
+
+$$ \oint_{\boldsymbol{r}}\boldsymbol{\tau}\cdot\boldsymbol{v}\;dl \;=\; \left(\partial_y v_z - \partial_z v_y\right)dy\,dz \;+\; \text{higher order} $$
+
+with $\boldsymbol{\tau}$ the unit tangent along the path. That is the $\hat{\boldsymbol{x}}$ component of the curl times the area of the rectangle, and rectangles perpendicular to $\hat{\boldsymbol{y}}$ and $\hat{\boldsymbol{z}}$ give the other two. So the curl is **net circulation per unit area**:
+
+$$ \hat{\boldsymbol{n}}\cdot\left(\nabla\times\boldsymbol{v}\right) \;=\; \lim_{S\to 0}\frac{\oint_{\boldsymbol{r}}\boldsymbol{\tau}\cdot\boldsymbol{v}\;dl}{A} $$
+
+where $A$ is the area of the open surface $S$ and $\hat{\boldsymbol{n}}$ is its unit normal, oriented so that a right-handed screw turned in the direction of $\boldsymbol{\tau}$ advances along $\hat{\boldsymbol{n}}$. Part 1 built the divergence from flux through a closed surface. The curl is built from circulation around a closed curve, one dimension down.
+
+### Task 6 — the operator, and the three flows again
+
+Three lines, in the pattern of the formula above. `np.gradient(Az, dy, axis=1)` is $\partial_y v_z$: the array holding the $z$-component, differentiated along the $y$-axis.
+
+Put the three fields of Task 3 through it. Field (b) rotates rigidly about the $z$-axis at $\omega = 1$ s$^{-1}$, so a paddle wheel dropped anywhere in it turns; fields (a) and (c) carry no rotation. Predict all three before running the cell.
+
+```{code-cell} ipython3
+# Task 6 -- two blanks. The x-component is given. Advance every letter one
+# step, x -> y -> z -> x, in the direction, the differentiation and the
+# component, and the other two lines write themselves.
+def curl(Ax, Ay, Az, dx, dy, dz):
+    cx = np.gradient(Az, dy, axis=1) - np.gradient(Ay, dz, axis=2)
+    cy = ___
+    cz = ___
+    return cx, cy, cz
+
+# --- given: the three fields of Task 3, through the new operator, and one
+#     more. Field (d) is the same rigid rotation turned onto the y-axis, so
+#     its curl is 2 y-hat: it is here to exercise the second line you wrote,
+#     which the three planar fields above leave at zero whatever you put in it.
+Ad = (Z, zero, -X)
+curl_a, curl_b = curl(*Aa, dx, dy, dz), curl(*Ab, dx, dy, dz)
+curl_c, curl_d = curl(*Ac, dx, dy, dz), curl(*Ad, dx, dy, dz)
+
+for name, w in [("(a) outward flow", curl_a), ("(b) rotation about z", curl_b),
+                ("(c) straining flow", curl_c), ("(d) rotation about y", curl_d)]:
+    print(f"{name:22s} curl = ({w[0].mean():+.2f}, {w[1].mean():+.2f}, {w[2].mean():+.2f})")
+
+fig, axes = plt.subplots(1, 3, figsize=(16, 4.6))
+for ax_, (name, A, w) in zip(axes, [("(a) outward flow", Aa, curl_a),
+                                    ("(b) rotation", Ab, curl_b),
+                                    ("(c) straining flow", Ac, curl_c)]):
+    fw.show_field_slice(X, Y, Z, *A[:2], background=w[2], ax=ax_, density=1.1,
+                        vmin=-3, vmax=3, colorbar=(ax_ is axes[-1]),
+                        label=r"$(\nabla\times\mathbf{A})_z$  [s$^{-1}$]", title=name)
+plt.tight_layout()
+plt.show()
+
+# --- self-check (leave this alone) ---
+# Each test reads the whole vector, not one component: a sign slip in `cy`
+# leaves every planar field untouched and would otherwise pass unnoticed.
+def _mag(w):
+    return np.sqrt(w[0]**2 + w[1]**2 + w[2]**2)
+
+fw.check_abs("(a) curl = 0 (outward flow)", _mag(curl_a), atol=1e-9)
+fw.check_abs("(c) curl = 0 (straining flow)", _mag(curl_c), atol=1e-9)
+fw.check_close("(b) curl = 2 z-hat", curl_b[2], 2.0, rtol=1e-6)
+fw.check_abs("...and nothing along x or y",
+             np.abs(curl_b[0]) + np.abs(curl_b[1]), atol=1e-9)
+fw.check_close("(d) curl = 2 y-hat", curl_d[1], 2.0, rtol=1e-6,
+               where=interior)
+fw.check_abs("...and nothing along x or z",
+             (np.abs(curl_d[0]) + np.abs(curl_d[2]))[interior], atol=1e-9)
+```
+
+:::{admonition} Solution — Task 6
+:class: dropdown
+
+```python
+    cy = np.gradient(Ax, dz, axis=2) - np.gradient(Az, dx, axis=0)
+    cz = np.gradient(Ay, dx, axis=0) - np.gradient(Ax, dy, axis=1)
+```
+:::
+
+:::{admonition} Two independent questions about one field
+:class: important
+
+The three flows of Task 3 now carry two answers each, and neither constrains the other:
+
+| | Field $\boldsymbol{A}$ | $\nabla\cdot\boldsymbol{A}$ | $\nabla\times\boldsymbol{A}$ |
+| :---: | :--- | :---: | :---: |
+| **(a)** | $x\,\hat{\boldsymbol{x}} + y\,\hat{\boldsymbol{y}} + z\,\hat{\boldsymbol{z}}$ | $3$ | $\boldsymbol{0}$ |
+| **(b)** | $-y\,\hat{\boldsymbol{x}} + x\,\hat{\boldsymbol{y}}$ | $0$ | $2\,\hat{\boldsymbol{z}}$ |
+| **(c)** | $x\,\hat{\boldsymbol{x}} - y\,\hat{\boldsymbol{y}}$ | $0$ | $\boldsymbol{0}$ |
+
+"How much is created here" and "how much does this spin here" are separate measurements. Field (c) answers zero to both and is still not the zero field: it stretches a fluid parcel along $x$ and squeezes it along $y$ at equal rates, changing its shape while conserving its volume and its orientation. Deformation is the third thing a flow can do, and neither operator reports it.
+
+Field (b) rotates at $\omega = 1$ s$^{-1}$ and its curl is $2\hat{\boldsymbol{z}}$. For rigid rotation at angular velocity $\boldsymbol{\omega}$ the curl is $2\boldsymbol{\omega}$, whatever the axis, which is why field (d) about $\hat{\boldsymbol{y}}$ returns $2\hat{\boldsymbol{y}}$. In fluid mechanics $\nabla\times\boldsymbol{v}$ is the **vorticity**, twice the local angular velocity of a fluid parcel.
+:::
+
+### Task 7 — closed streamlines are not curl
+
+Task 3 established that arrows spreading apart do not make a divergence. The same warning applies here, in the same shape, and the two errors are the same error. **The picture a streamline makes says nothing about the curl.**
+
+Two fields settle it, with the rotation of Task 6 as a control. Both are written on the distance from the $z$-axis, the cylindrical $\varrho = \sqrt{x^2+y^2}$, which is not the spherical $r$ of Parts 1 and 2.
+
+| | Field | Streamlines |
+| :---: | :--- | :--- |
+| **rotation** | $\omega\,\varrho\,\hat{\boldsymbol{\phi}} \;=\; -y\,\hat{\boldsymbol{x}} + x\,\hat{\boldsymbol{y}}$ | concentric circles |
+| **shear** | $\sigma\,y\,\hat{\boldsymbol{x}}$ | straight lines, all parallel to $x$ |
+| **line vortex** | $\dfrac{\Gamma_0}{2\pi\varrho}\hat{\boldsymbol{\phi}} \;=\; \Gamma_0\dfrac{-y\,\hat{\boldsymbol{x}} + x\,\hat{\boldsymbol{y}}}{2\pi\varrho^{2}}$ | concentric circles |
+
+with $\omega = 1$ s$^{-1}$, shear rate $\sigma = 1$ s$^{-1}$ and $\Gamma_0 = 2\pi$ m$^2$/s, so all three curls come out in s$^{-1}$ and one colour scale serves the row.
+
+The shear is the water beside a riverbank, or between two plates sliding past each other: further out it runs faster, but every parcel travels in a straight line and none of them goes round anything. The line vortex circles the axis exactly as the rotation does, and falls off as $1/\varrho$. Record your prediction of the sign of $(\nabla\times\boldsymbol{v})_z$ for each, then measure.
+
+```{code-cell} ipython3
+# +1 for anticlockwise rotation, -1 for clockwise, 0 for none.
+curl_predictions = {"rotation": ___, "shear": ___, "line vortex": ___}
+```
+
+```{code-cell} ipython3
+# --- given: distance from the z-axis, and a test region that avoids it ---
+varrho = np.sqrt(X**2 + Y**2)
+varrho_s = np.maximum(varrho, 1e-12)      # the axis kept out of the denominators
+ring = interior & (varrho > 0.4) & (varrho < 1.6)
+
+# Task 7 -- two blanks, one field each. Both lie in the z = 0 plane, so the
+# third component of each is `zero`. Take Gamma_0 / 2*pi = 1, as above.
+A_shear = ___                             # y x-hat, as a triple
+A_vortex = ___                            # (-Y, X) / varrho_s**2, as a triple
+
+# --- given ---
+curl_shear = curl(*A_shear, dx, dy, dz)
+curl_vortex = curl(*A_vortex, dx, dy, dz)
+
+# Reported as the scale-free ratio Task 2 built for the divergence, with the
+# distance from the AXIS in place of the distance from the origin: the size of
+# the curl against |v|/varrho, the size a derivative of this field would have
+# if nothing cancelled. Same reasoning, same reading: 1 means no cancellation.
+v_mag = np.sqrt(A_vortex[0]**2 + A_vortex[1]**2)
+vortex_ratio = np.median(np.abs(curl_vortex[2][ring]) / (v_mag / varrho_s)[ring])
+measured = {"rotation": curl_b[2].mean(), "shear": curl_shear[2][ring].mean(),
+            "line vortex": curl_vortex[2][ring].mean()}
+print(f"rotation    : curl_z = {measured['rotation']:+.3f} s^-1")
+print(f"shear       : curl_z = {measured['shear']:+.3f} s^-1")
+print(f"line vortex : |curl_z| / (|v|/varrho) = {vortex_ratio:.2%}  (zero, to grid error)")
+# The vortex is zero only to grid error, so the score needs a deadband:
+# anything under 5% of the largest curl in the row counts as no rotation.
+deadband = 0.05 * max(abs(v) for v in measured.values())
+for key, got in measured.items():
+    sign = 0 if abs(got) < deadband else int(np.sign(got))
+    verdict = "as predicted" if curl_predictions[key] == sign else "NOT what you predicted"
+    print(f"  {key:12s} you said {curl_predictions[key]:+d}, measured {sign:+d}  --  {verdict}")
+
+# The vortex is singular on the z-axis. What the operator returns on the few
+# cells around it is the grid's difficulty and not the field's, so those cells
+# are left blank rather than allowed to dominate the panel.
+shown = np.where(varrho < 0.3, np.nan, curl_vortex[2])
+
+fig, axes = plt.subplots(1, 3, figsize=(16, 4.6))
+panels = [("rotation", Ab, curl_b[2]), ("shear", A_shear, curl_shear[2]),
+          ("line vortex", A_vortex, shown)]
+for ax_, (name, A, w) in zip(axes, panels):
+    fw.show_field_slice(X, Y, Z, *A[:2], background=w, ax=ax_, density=1.1,
+                        vmin=-3, vmax=3, colorbar=(ax_ is axes[-1]),
+                        label=r"$(\nabla\times\mathbf{v})_z$  [s$^{-1}$]", title=name)
+plt.tight_layout()
+plt.show()
+
+# --- self-check (leave this alone) ---
+fw.check_close("shear: curl = -1 z-hat, although nothing goes round",
+               curl_shear[2][ring], -1.0, rtol=1e-6)
+fw.check(f"line vortex: curl = 0, although everything goes round ({vortex_ratio:.2%})",
+         vortex_ratio < 0.05)
+fw.check("...and it really does circle the axis: v has no radial component",
+         np.max(np.abs((A_vortex[0]*X + A_vortex[1]*Y)[ring])) < 1e-12)
+```
+
+:::{admonition} Solution — Task 7
+:class: dropdown
+
+```python
+A_shear = (Y, zero, zero)
+A_vortex = (-Y / varrho_s**2, X / varrho_s**2, zero)
+```
+:::
+
+:::{admonition} What the paddle wheel actually measures
+:class: warning
+
+Drop a small paddle wheel into each flow and watch its axle.
+
+In the **shear** it spins, at half a radian per second, clockwise. The water above the wheel runs faster than the water below, so the top blades are pushed harder than the bottom ones. Nothing in the flow travels in a circle and the curl is still $-\hat{\boldsymbol{z}}$.
+
+In the **line vortex** the wheel is carried once round the axis and comes back pointing the way it started. It orbits without spinning, like the Moon in reverse. Two separate effects cancel, and the cylindrical formula separates them. For a field $v_\phi(\varrho)\,\hat{\boldsymbol{\phi}}$,
+
+$$ (\nabla\times\boldsymbol{v})_z = \frac{1}{\varrho}\frac{\partial\left(\varrho\, v_\phi\right)}{\partial\varrho} - \frac{1}{\varrho}\frac{\partial v_\varrho}{\partial \phi} \;=\; \underbrace{\frac{dv_\phi}{d\varrho}}_{\text{shear}} + \underbrace{\frac{v_\phi}{\varrho}}_{\text{orbit}} $$
+
+the second term dropping because these fields have no radial component. The **shear** term is the blades: the inner ones sit in faster water than the outer ones, which turns the wheel backwards, at $-\Gamma_0/2\pi\varrho^{2}$. The **orbit** term is the wheel's own frame turning once per lap, forwards, at $+\Gamma_0/2\pi\varrho^{2}$. Their sum is zero at $1/\varrho$ and at no other falloff: a steeper $1/\varrho^{2}$ over-cancels and spins the wheel backwards, a shallower one spins it forwards. Their *mean* is the local angular velocity, which is where Task 6's factor of two comes from.
+
+Read the same formula the other way and the curl vanishes when $\varrho\,v_\phi$ is constant, which is $v_\phi \propto 1/\varrho$ and nothing else. Rigid rotation, $v_\phi = \omega\varrho$, gives $2\omega$ instead.
+
+That single surviving field, $\boldsymbol{H} = \dfrac{I}{2\pi\varrho}\hat{\boldsymbol{\phi}}$, is the magnetic field around a straight wire carrying a current $I$. Its curl is zero at every point outside the wire, and the current is still there. The end of this part explains how both can be true.
+:::
+
+### Task 8 — circulation per unit area, and Stokes' theorem
+
+A real vortex has a core. Stirred coffee, a tornado and the vortex trailing from a wing all rotate almost rigidly near the axis and fall off as $1/\varrho$ far from it, because viscosity spreads the vorticity over a finite radius $b$. The **Lamb–Oseen vortex** is the exact solution for that spreading, with $b^{2} = 4\nu t$ after a time $t$ in a fluid of kinematic viscosity $\nu$:
+
+$$ v_\phi(\varrho) = \frac{\Gamma}{2\pi\varrho}\left(1 - e^{-\varrho^{2}/b^{2}}\right), \qquad \Gamma = 1\ \text{m}^2\text{/s}, \qquad b = 0.5\ \text{m}. $$
+
+Inside the core this is $\Gamma\varrho/2\pi b^{2}$, the rigid rotation of Task 6. Outside it is $\Gamma/2\pi\varrho$, the irrotational vortex of Task 7. The cylindrical formula turns it into a vorticity that is a Gaussian blob:
+
+$$ (\nabla\times\boldsymbol{v})_z = \frac{1}{\varrho}\frac{d}{d\varrho}\left(\varrho\,v_\phi\right) = \frac{\Gamma}{\pi b^{2}}\,e^{-\varrho^{2}/b^{2}} $$
+
+the same shape as Task 4's blob of charge, with $\Gamma$ in the part of $Q$. The rest of this task is Part 2 run one dimension down: a closed curve instead of a closed surface, circulation instead of flux, and **Stokes' theorem** instead of the divergence theorem,
+
+$$ \oint_{\boldsymbol{r}}\boldsymbol{\tau}\cdot\boldsymbol{v}\;dl \;=\; \int_{\boldsymbol{r}\in S}\hat{\boldsymbol{n}}\cdot\left(\nabla\times\boldsymbol{v}\right)dS $$
+
+```{code-cell} ipython3
+Gamma, b_core = 1.0, 0.5          # circulation [m^2/s], core radius b [m]
+
+# --- given: the vortex, and the vorticity it should have ---
+_swirl = np.where(varrho < 1e-8, Gamma / (2*np.pi*b_core**2),
+                  Gamma / (2*np.pi*varrho_s**2) * (1 - np.exp(-varrho**2/b_core**2)))
+oseen = (-Y * _swirl, X * _swirl, zero)   # v_phi phi-hat, in Cartesian components
+w_exact = Gamma / (np.pi * b_core**2) * np.exp(-varrho**2 / b_core**2)
+curl_oseen = curl(*oseen, dx, dy, dz)
+
+# Task 8 -- two blanks, one per side of Stokes' theorem.
+#
+# LEFT SIDE. Walk the four edges of a rectangle in the z = 0 plane once
+# counter-clockwise, so the right-hand rule puts the unit normal along +z-hat,
+# and add up the component of the field along the direction of travel:
+#
+#     edge     samples along it   travelling   spacing   sign
+#     y = y0   Ax[sx, iy0, k]        +x          dx       +
+#     x = x1   Ay[ix1, sy, k]        +y          dy       +
+#     y = y1   Ax[sx, iy1, k]        -x          dx       -
+#     x = x0   Ay[ix0, sy, k]        -y          dy       -
+#
+# The two edges walked backwards enter negatively, exactly as the inward faces
+# did in Task 5. The x pair is written for you.
+
+def loop_circulation(Ax, Ay, x0, x1, y0, y1):
+    """Counter-clockwise circulation of (Ax, Ay) round a rectangle in z = 0."""
+    ix0, ix1 = [int(np.argmin(np.abs(axis - q))) for q in (x0, x1)]
+    iy0, iy1 = [int(np.argmin(np.abs(axis - q))) for q in (y0, y1)]
+    sx, sy, k = slice(ix0, ix1 + 1), slice(iy0, iy1 + 1), fw.z0_index(Z)
+    along_x = (fw.line_integral(Ax[sx, iy0, k], dx)
+               - fw.line_integral(Ax[sx, iy1, k], dx))
+    along_y = ___
+    return along_x + along_y
+
+
+# RIGHT SIDE. n-hat is +z-hat, so only the z-component of the curl crosses the
+# rectangle. Integrate it over the flat patch the same indices span: one call
+# to fw.area_integral, on the z = 0 plane, with spacings dx and dy.
+
+def curl_flux(x0, x1, y0, y1):
+    """Flux of the vorticity through the same rectangle: Stokes' right side."""
+    ix0, ix1 = [int(np.argmin(np.abs(axis - q))) for q in (x0, x1)]
+    iy0, iy1 = [int(np.argmin(np.abs(axis - q))) for q in (y0, y1)]
+    return ___
+
+
+# --- given: the limit in the definition, run as a measurement. Every side
+#     below is a whole number of grid spacings, so each loop is the one asked
+#     for rather than the nearest one the grid happens to be able to draw.
+def curl_z_area(A, x0, x1, y0, y1):
+    """Circulation per unit area round the same rectangle."""
+    return loop_circulation(A[0], A[1], x0, x1, y0, y1) / ((x1 - x0) * (y1 - y0))
+
+w_centre = curl_oseen[2][c, c, fw.z0_index(Z)]
+print(f"vorticity at the origin: measured {w_centre:.4f} s^-1, "
+      f"exact {Gamma/(np.pi*b_core**2):.4f} s^-1")
+print(f"the whole Gaussian, worst error "
+      f"{np.abs(curl_oseen[2] - w_exact)[interior].max()/w_exact.max():.2%} of peak\n")
+print(f"{'side [m]':>9} {'samples/edge':>13} {'circulation':>13} {'C/A':>9}"
+      f" {'C/A / measured':>15}")
+for m in (18, 9, 6, 3, 2, 1):
+    h = m * dx
+    print(f"{2*h:9.4f} {2*m+1:13d} {loop_circulation(*oseen[:2], -h, h, -h, h):13.5f}"
+          f" {curl_z_area(oseen, -h, h, -h, h):9.4f}"
+          f" {curl_z_area(oseen, -h, h, -h, h)/w_centre:15.4f}")
+
+# --- given: Stokes' theorem on five rectangles, two of them off-centre ---
+rects = [("centred, side 0.8", (-0.4, 0.4, -0.4, 0.4)),
+         ("centred, side 2.0", (-1.0, 1.0, -1.0, 1.0)),
+         ("centred, side 3.2", (-1.6, 1.6, -1.6, 1.6)),
+         ("off-centre, over the core", (-0.2, 1.4, -0.6, 1.0)),
+         ("off to one side", (0.4, 1.6, -0.6, 0.6))]
+print(f"\n  {'rectangle':<26} {'circulation':>12} {'flux of curl':>13} {'apart':>8}")
+for name, lim in rects:
+    C, S = loop_circulation(*oseen[:2], *lim), curl_flux(*lim)
+    print(f"  {name:<26} {C:12.5f} {S:13.5f} {abs(C - S)/abs(C):8.2%}")
+print(f"\n  all the vorticity there is: Gamma = {Gamma:.4f} m^2/s")
+
+fig, axes = plt.subplots(1, 3, figsize=(16, 4.2))
+fw.show_field_slice(X, Y, Z, *oseen[:2], background=curl_oseen[2], ax=axes[0],
+                    density=1.2, vmin=0, vmax=1.3, levels=14, cmap="inferno",
+                    symmetric=False, stream_color="w",
+                    label=r"$(\nabla\times\mathbf{v})_z$  [s$^{-1}$]",
+                    title="the vortex, over its vorticity")
+prof = np.linspace(1e-3, 2.0, 400)
+axes[1].plot(prof, Gamma/(2*np.pi*prof)*(1 - np.exp(-prof**2/b_core**2)), "k", lw=1.8,
+             label=r"$v_\phi(\varrho)$")
+axes[1].plot(prof, Gamma*prof/(2*np.pi*b_core**2), "C1--", lw=1.2,
+             label=r"core: $\Gamma\varrho/2\pi b^2$")
+axes[1].plot(prof, Gamma/(2*np.pi*prof), "C2:", lw=1.4, label=r"outside: $\Gamma/2\pi\varrho$")
+axes[1].set_ylim(0, 0.25); axes[1].set_title(r"rigid inside the core, $1/\varrho$ outside")
+axes[1].set_xlabel(r"$\varrho$ [m]"); axes[1].set_ylabel(r"$v_\phi$ [m s$^{-1}$]")
+
+# the claim under test, drawn: measured vorticity against the exact Gaussian
+k0 = fw.z0_index(Z)
+axes[2].plot(X[:, c, k0], w_exact[:, c, k0], "k", lw=2.4, alpha=0.35, label="exact")
+axes[2].plot(X[:, c, k0], curl_oseen[2][:, c, k0], "C3", lw=1.2, label="measured")
+axes[2].set_ylim(0, 1.45); axes[2].set_title("vorticity along $y = 0$")
+axes[2].set_xlabel("$x$ [m]"); axes[2].set_ylabel(r"$(\nabla\times\mathbf{v})_z$ [s$^{-1}$]")
+for ax_ in axes[1:]:
+    ax_.axvline(b_core, color="C0", lw=1, alpha=0.6)
+    ax_.grid(alpha=0.3); ax_.legend(fontsize=8)
+axes[1].annotate(r"$\varrho = b$", (b_core + 0.05, 0.02), color="C0", ha="left")
+plt.tight_layout()
+plt.show()
+
+# --- self-check (leave this alone) ---
+_h = 6 * dx
+fw.check_scalar("Stokes: circulation = flux of the curl through the loop",
+                loop_circulation(*oseen[:2], -_h, _h, -_h, _h),
+                curl_flux(-_h, _h, -_h, _h), rtol=0.01, unit=" m^2/s")
+fw.check_scalar("...and again on a rectangle that is not centred on the vortex",
+                loop_circulation(*oseen[:2], -0.2, 1.4, -0.6, 1.0),
+                curl_flux(-0.2, 1.4, -0.6, 1.0), rtol=0.01, unit=" m^2/s")
+_small = curl_z_area(oseen, -dx, dx, -dx, dx)
+fw.check(f"circulation per unit area -> the vorticity at the centre "
+         f"({_small:.4f} against {w_centre:.4f} s^-1)",
+         abs(_small - w_centre) < 0.01 * abs(w_centre))
+_wide = loop_circulation(*oseen[:2], -1.6, 1.6, -1.6, 1.6)
+fw.check(f"a loop well outside the core collects all of Gamma "
+         f"({_wide:.4f} of {Gamma:.4f} m^2/s)", abs(_wide - Gamma) < 0.01 * Gamma)
+```
+
+:::{admonition} Solution — Task 8
+:class: dropdown
+
+```python
+# in loop_circulation, the left side:
+    along_y = (fw.line_integral(Ay[ix1, sy, k], dy)
+               - fw.line_integral(Ay[ix0, sy, k], dy))
+
+# in curl_flux, the right side:
+    return fw.area_integral(curl_oseen[2][ix0:ix1+1, iy0:iy1+1, fw.z0_index(Z)],
+                            dx, dy)
+```
+:::
+
+:::{admonition} The same theorem, one dimension down
+:class: important
+
+Read the first table downwards. The loop shrinks, the circulation falls, the area falls faster, and the ratio climbs to 0.996 of the vorticity measured at the centre: 0.137 of it at side 2.4 m, 0.908 at side 0.4 m. That is the limit in the definition, evaluated rather than asserted. The comparison is against the **measured** 1.2620 s$^{-1}$ rather than the exact 1.2732, so the last 0.9% is the grid error already reported above and not a failure of the limit.
+
+Read the second table across. Five rectangles, two of them not centred on the vortex, and the two sides of Stokes' theorem agree to better than 1% on every one: to a few parts in $10^{5}$ on the largest centred loop, and worst on the smallest, which is only 13 samples across. Size, not placement, is what sets the accuracy, and it is the same second-order error Task 5 measured on the divergence theorem. The left side never looks inside the loop and the right side never looks at the boundary.
+
+| | Divergence theorem | Stokes' theorem |
+| :--- | :--- | :--- |
+| Boundary integral | flux through a closed **surface** | circulation round a closed **curve** |
+| Interior integral | $\nabla\cdot\boldsymbol{v}$ over the enclosed **volume** | $\hat{\boldsymbol{n}}\cdot(\nabla\times\boldsymbol{v})$ over the enclosed **area** |
+| Source it counts | $Q_{\text{enc}}/\varepsilon_0$ | $\Gamma$, or the enclosed current |
+
+The largest loop returns 0.9999 of $\Gamma$, exactly as the largest boxes of Task 5 weighed the whole 1 nC. Enlarging a loop that already encloses all the vorticity adds nothing, for the same reason that charge outside a closed surface contributes nothing.
+:::
+
+### The wire, and Ampère's law
+
+Shrink the vortex core to nothing, $b\to 0$, and $v_\phi$ becomes the irrotational $\Gamma/2\pi\varrho$ of Task 7 at every radius, with all the vorticity compressed onto the axis. The Gaussian collapses to a Dirac delta, as the Gaussian blob of charge did when Part 2 shrank it to a point.
+
+That field is the magnetic field of a straight wire carrying a current $I$ along $\hat{\boldsymbol{z}}$, and the statement relating them is **Ampère's law**, in the differential and integral forms Stokes' theorem connects:
+
+$$ \nabla\times\boldsymbol{H} = \boldsymbol{J} \qquad\Longleftrightarrow\qquad \oint_{\boldsymbol{r}}\boldsymbol{\tau}\cdot\boldsymbol{H}\;dl = \int_{\boldsymbol{r}\in S}\hat{\boldsymbol{n}}\cdot\boldsymbol{J}\;dS = I_{\text{enc}} $$
+
+with $\boldsymbol{H}$ in A/m, $\boldsymbol{J}$ in A/m$^2$, and $\hat{\boldsymbol{n}}$ fixed by the right-hand rule from the direction of travel. `loop_circulation` walks counter-clockwise in the $z=0$ plane, so $\hat{\boldsymbol{n}} = \hat{\boldsymbol{z}}$ and a positive answer means current flowing towards the reader.
+
+```{code-cell} ipython3
+I_wire = 1.0                              # current along +z, in amperes
+Hx, Hy = -Y * I_wire / (2*np.pi*varrho_s**2), X * I_wire / (2*np.pi*varrho_s**2)
+
+print(f"  {'loop':<28} {'oint tau.H dl [A]':>18}")
+for name, lim in [("encloses the wire, side 0.8", (-0.4, 0.4, -0.4, 0.4)),
+                  ("encloses the wire, side 2.0", (-1.0, 1.0, -1.0, 1.0)),
+                  ("encloses the wire, side 3.2", (-1.6, 1.6, -1.6, 1.6)),
+                  ("encloses it, lopsidedly    ", (-0.4, 1.6, -1.0, 0.6)),
+                  ("misses the wire            ", (0.4, 1.6, -0.6, 0.6)),
+                  ("misses it, and is large    ", (0.2, 1.8, -1.8, 1.8))]:
+    print(f"  {name:<28} {loop_circulation(Hx, Hy, *lim):18.5f}")
+print(f"\n  current actually in the wire: {I_wire:.5f} A")
+```
+
+:::{admonition} A loop that encloses nothing measurable
+:class: important
+
+Every loop enclosing the wire returns $I$ to better than two parts in a thousand, at any size and whether or not it is centred on the wire. Every loop missing it returns at most $4\times10^{-4}$ A, which against the enclosing loops' 1 A is zero. The circulation counts what passes through the loop and nothing else, exactly as the closed surface of Part 2 counted the charge inside and nothing else. Shape-independence is part of the same statement, though `loop_circulation` draws only rectangles and cannot demonstrate it; it follows from Stokes' theorem, since two loops enclosing the same current bound surfaces carrying the same flux of $\boldsymbol{J}$.
+
+Now put that beside Task 7. At every point these loops pass through, $\nabla\times\boldsymbol{H} = 0$: the field is irrotational everywhere the grid can sample it, and the loop integral is 1 A regardless. There is no contradiction. Stokes' theorem equates the circulation to the flux of $\boldsymbol{J}$ through the loop, and $\boldsymbol{J}$ is zero over the whole of the loop's interior except one line, where it is infinite. The current density is a Dirac delta on the axis, the integral form survives it, and the differential form does not, which is what happened to $\rho_v$ at the point charge in Part 2.
+
+A real wire has a finite radius and a finite $\boldsymbol{J}$ spread over its cross-section, and then both forms hold everywhere. The homework builds that wire.
+:::
+
+### Why the electric fields had a potential
+
+$\boldsymbol{E} = -\nabla V$ was written in Lab 1 without asking whether an arbitrary vector field can be written that way. The rotation, the shear and the vortex above cannot. The curl is the test:
+
+$$ \nabla\times\left(\nabla p\right) = \boldsymbol{0} \qquad \text{for every twice-differentiable } p $$
+
+because each component subtracts a pair of mixed second derivatives, $\partial_x\partial_y p - \partial_y\partial_x p$, and mixed partials commute. Run the two electric fields of this notebook through the curl.
+
+```{code-cell} ipython3
+E_point = tuple(np.nan_to_num(q) for q in (Ex, Ey, Ez))   # built as -grad V
+curl_point = curl(*E_point, dx, dy, dz)
+curl_blob = curl(Ex_b, Ey_b, Ez_b, dx, dy, dz)            # built from E_r(r) r-hat
+
+shell = interior & (r > 0.5) & (r < 1.6)
+for name, w, A in [("point charge, from -grad V", curl_point, E_point),
+                   ("blob, from the analytic E_r", curl_blob, (Ex_b, Ey_b, Ez_b))]:
+    # index first: |E| is zero inside the mask, and 0/0 there would warn
+    wm = np.sqrt(w[0]**2 + w[1]**2 + w[2]**2)[shell]
+    Am = (np.sqrt(A[0]**2 + A[1]**2 + A[2]**2) / rs)[shell]
+    print(f"  {name:28s} median |curl E| / (|E|/r) = {np.median(wm / Am):.2e}")
+
+# The circulation, on the blob field, which carries no mask for a loop to cross.
+# Reported against the natural scale for a voltage here: the strongest field on
+# the grid, carried along one metre of path.
+scale_V = float(np.max(np.sqrt(Ex_b**2 + Ey_b**2 + Ez_b**2)))
+print(f"\n  {'loop':<20} {'circulation of E':>18} {'/ (|E|max x 1 m)':>18}")
+for name, lim in [("centred, side 2.0", (-1.0, 1.0, -1.0, 1.0)),
+                  ("off-centre", (-0.2, 1.4, -0.6, 1.0)),
+                  ("off to one side", (0.4, 1.6, -0.6, 0.6))]:
+    circ = loop_circulation(Ex_b, Ey_b, *lim)
+    print(f"  {name:<20} {circ:+15.2e} V {abs(circ)/scale_V:18.1e}")
+print(f"\n  the same square, side 2.0, on the wire above: "
+      f"{loop_circulation(Hx, Hy, -1.0, 1.0, -1.0, 1.0):.3f} A")
+```
+
+:::{admonition} Why voltage is a number and not a route
+:class: important
+
+The two fields report zero curl with thirteen orders of magnitude between them, and the gap is in how each was built rather than in the physics. `Ex, Ey, Ez` came out of `-np.gradient(V)`, and the curl subtracts the same centred differences in the opposite order; the stencil obeys the identity as strictly as the algebra does, so nothing survives but the order in which floating-point numbers were added, around $10^{-16}$. The blob field was built from the analytic $E_r(r)$ and never passed through a numerical gradient, so it shows the 0.6% that a centred difference costs on this grid. Neither number measures the physics. Both are consistent with the one statement being tested.
+
+The circulations say the same thing on a closed curve: a few parts in $10^{4}$ of the natural voltage scale, dropping to round-off on the centred loop, whose symmetry cancels it exactly. Put the last line beside them. Same integrator, same grid, same size of loop, and the wire returns a full ampere.
+
+Zero circulation is what makes potential a usable idea. Carrying a charge round a circuit and back to its starting point costs no net work, so the work done between two points is independent of the route, and one number can be attached to each point. That number is $V$.
+
+Two restrictions are worth naming, and Part 3 has already demonstrated both.
+
+**Zero curl gives a potential only where the region has no holes in it.** The wire is the exception: $\nabla\times\boldsymbol{H} = \boldsymbol{0}$ at every point outside it, and $\oint\boldsymbol{\tau}\cdot\boldsymbol{H}\,dl = I \neq 0$. A loop encircling the axis cannot be shrunk to a point without crossing the current, so there is nothing for Stokes' theorem to integrate the curl over, and no single-valued potential for $\boldsymbol{H}$ exists out there. Around a point charge, by contrast, the punctured space *is* simply connected and $V$ survives.
+
+**$\nabla\times\boldsymbol{E} = \boldsymbol{0}$ holds in electrostatics.** When the magnetic field changes with time, $\nabla\times\boldsymbol{E} = -\partial\boldsymbol{B}/\partial t$, the circulation round a loop is no longer zero, and that circulation is the voltage a generator produces. At that point $V$ alone stops being enough, which is where this course is going.
+:::
+
+---
+
+## Closing
+
+The chain built across the two labs, in one line:
+
+$$ \rho_v \;\longrightarrow\; V \;\xrightarrow{\ -\nabla\ }\; \boldsymbol{E} \;\xrightarrow{\ \nabla\cdot\ }\; \rho_v/\varepsilon_0, \qquad \nabla\times\boldsymbol{E} = \boldsymbol{0} $$
+
+with the last statement the licence for the arrow labelled $-\nabla$: only a field with zero curl has a potential to be recovered from.
+
+- **Gradient.** Scalar in, vector out. Points along steepest increase, normal to the level surfaces, with length equal to the rate of increase.
+- **Divergence.** Vector in, scalar out. Net flux per unit volume, which measures what is created at a point and nothing else.
+- **Curl.** Vector in, vector out. Net circulation per unit area, about the axis its own direction gives, which measures local rotation and not the shape of a streamline.
+
+Each of the last two comes with an integral theorem, and each theorem replaces a derivative that fails at a singular source with an integral that does not. The point charge and the current-carrying wire are the same difficulty met twice.
+
+### The same three operators, elsewhere in ECT
+
+Electrostatics is a convenient place to learn these, not the only place to use them. Each row below gives a potential, its gradient, and a statement about sources. The numerical machinery written in these two labs applies unchanged to all of them:
+
+| System | Potential | Field | Source equation |
+| :--- | :--- | :--- | :--- |
+| Electrostatics | $V$ [V] | $\boldsymbol{E} = -\nabla V$ &nbsp; [V/m] | $\nabla\cdot\boldsymbol{E} = \rho_v/\varepsilon_0$ |
+| Gravitation | $\Phi$ [J/kg] | $\boldsymbol{g} = -\nabla \Phi$ &nbsp; [m/s$^2$] | $\nabla\cdot\boldsymbol{g} = -4\pi G\rho_m$ |
+| Heat conduction | $T$ [K] | $\boldsymbol{q}_T = -k\nabla T$ &nbsp; [W/m$^2$] | $\nabla\cdot\boldsymbol{q}_T = 0$ (steady, no sources) |
+| Groundwater flow | $h$ [m] | $\boldsymbol{q}_h = -K\nabla h$ &nbsp; [m/s] | $\nabla\cdot\boldsymbol{q}_h = 0$ (steady, incompressible) |
+
+with $G$ the gravitational constant, $\rho_m$ the mass density [kg m$^{-3}$], $k$ the thermal conductivity [W m$^{-1}$ K$^{-1}$] and $K$ the hydraulic conductivity [m/s].
+
+The minus signs are all the same minus sign: heat flows from hot to cold, water flows from high head to low, a positive charge falls from high potential to low. Flow runs downhill, and the gradient points uphill.
+
+The last two rows show why solenoidal fields matter in practice. $\nabla\cdot\boldsymbol{q}_h = 0$ in an aquifer states conservation of water locally, in the form a numerical model actually solves.
+
+In a **homogeneous** medium, where $k$ and $K$ are constants, every field in that table is a gradient, so every one of them has zero curl and none can circulate. Let $K$ vary from place to place, as it does in any real aquifer, and $\nabla\times\boldsymbol{q}_h = -\nabla K\times\nabla h$ need not vanish. The fields that circulate are the ones with no potential to be had, and they are the subject of the rest of the course:
+
+| Field | Circulation equation | What sets it |
+| :--- | :--- | :--- |
+| Magnetic field $\boldsymbol{H}$ [A/m] | $\nabla\times\boldsymbol{H} = \boldsymbol{J}$ | the current threading the loop |
+| Fluid velocity $\boldsymbol{v}$ [m/s] | $\nabla\times\boldsymbol{v} = \boldsymbol{\omega}_v$ | shear at a boundary, and rotation of the Earth |
+| Electric field, unsteady | $\nabla\times\boldsymbol{E} = -\partial\boldsymbol{B}/\partial t$ | a magnetic field that changes with time |
+
+where $\boldsymbol{\omega}_v \equiv \nabla\times\boldsymbol{v} = 2\boldsymbol{\omega}$ is the vorticity, twice the local angular velocity of Task 6.
+
+The first row is Ampère's law in the static limit, and the term Maxwell added to it, $\partial\boldsymbol{D}/\partial t$, is the reason light exists. The third is Faraday's law, where the potential $V$ stops being sufficient on its own. Those two together with the two source equations of Parts 1 and 2 are Maxwell's four.
+
+### Formative assessment — Chapters 1 and 2
+
+Not graded, and not handed in. It exists so you can find out what you do not yet know, while there is still time to fix it. Allow about **45 minutes**: 12 for Part A and the rest for Part B.
+
+The two chapters end here. Part A checks that you can say what the operators mean; Part B gives you a system nobody has solved for you and asks you to measure all three.
+
+#### Part A — six questions, no code
+
+:::{admonition} A1. Summing and truncating
+:class: tip
+
+The bouncing ball converged: every extra term in $\sum T_n$ brought the answer closer to $T_\infty$. The series for $(1+x)^{-1}$ did not, once $\lvert x\rvert > 1$: no number of terms helps. Both are infinite sums of shrinking-looking terms. What distinguishes them, and how would you decide which case you are in before spending an afternoon adding terms?
+:::
+
+:::{admonition} A2. Which coordinates, and which symbol
+:class: tip
+
+You are handed three fields: the temperature around a buried sphere; the magnetic field around a long straight cable; the field in a rectangular room. Which coordinate system would you compute each in, and why? Then: in cylindrical coordinates the radial distance is written $\varrho$ and in spherical it is written $r$. Give one calculation that goes wrong if you conflate them.
+:::
+
+:::{admonition} A3. The gradient of a distance
+:class: tip
+
+Without computing anything: what is $\lvert\nabla r\rvert$, and why must it be that number for every $r > 0$? What direction does $\nabla r$ point, and what does that say about the surfaces $r = \text{constant}$?
+:::
+
+:::{admonition} A4. Two fields that look like sources
+:class: tip
+
+$\boldsymbol{A} = x\,\hat{\boldsymbol{x}} - y\,\hat{\boldsymbol{y}}$ has arrows that fly apart along the $x$-axis, and $\nabla\cdot\boldsymbol{A} = 0$. $\boldsymbol{E}$ outside a charged blob has arrows that fly apart in every direction, and $\nabla\cdot\boldsymbol{E} = 0$ as well. Are these the same statement twice? Explain each with a box, not with algebra.
+:::
+
+:::{admonition} A5. The wire that circulates without curling
+:class: tip
+
+Outside a straight current-carrying wire, $\nabla\times\boldsymbol{H} = \boldsymbol{0}$ at every point you can measure, and yet $\oint\boldsymbol{\tau}\cdot\boldsymbol{H}\,dl = I \neq 0$ around any loop enclosing it. Stokes' theorem says these two are equal. Resolve it. Then say what goes wrong if you try to define a potential for $\boldsymbol{H}$ outside the wire.
+:::
+
+:::{admonition} A6. Why keep both forms
+:class: tip
+
+Each operator came with an integral theorem. Name the one situation, met twice in these labs, in which the differential form fails and the integral form still works, and say what the integral form is doing that the derivative cannot.
+:::
+
+#### Part B — a buried heating panel
+
+A rectangular electrical heating element, $1.0 \times 0.6$ m, is buried in soil and dissipates $P = 100$ W. Nothing here has been solved for you; the tools are the ones you built.
+
+A steady point source of power $P$ in a medium of thermal conductivity $k$ raises the temperature above ambient by $P/4\pi k r$, the same $1/r$ used throughout these labs. Split the panel into $N = 20\times12$ sub-sources, give each an equal share of the power, and superpose:
+
+$$ T(\boldsymbol{r}) = \frac{P}{4\pi k N}\sum_{i=1}^{N}\frac{1}{\lvert\boldsymbol{r}-\boldsymbol{r}_i\rvert}, \qquad P = 100\ \text{W}, \qquad k_{\text{soil}} = 1.5\ \text{W m}^{-1}\text{K}^{-1} $$
+
+Check the dimensions before you code. $[P]/[k] = \text{W}/(\text{W m}^{-1}\text{K}^{-1}) = \text{m}\cdot\text{K}$, divided by a distance, so $T$ comes out in kelvin. A temperature formula that does not reduce to kelvin has an error in it.
+
+```{code-cell} ipython3
+# --- given: the panel, and the grid it sits on (the cube from Part 0) ---
+P_heat, k_soil = 100.0, 1.5               # W, and W/m/K for soil
+panel_x = np.linspace(-0.5, 0.5, 20)      # sub-source positions, in the z = 0 plane
+panel_y = np.linspace(-0.3, 0.3, 12)
+N_sub = panel_x.size * panel_y.size
+
+# B1 -- two blanks, inside the loop. `d_min` records how close each grid point
+# comes to the nearest sub-source; the mask below uses it.
+T_sum = np.zeros_like(X)
+d_min = np.full(X.shape, np.inf)
+for x0 in panel_x:
+    for y0 in panel_y:
+        d = ___                           # distance from (x0, y0, 0) to every grid point
+        d_min = np.minimum(d_min, d)
+        T_sum += 1.0 / np.maximum(d, 1e-12)
+T_panel = ___                             # the prefactor, applied once at the end
+
+# --- given: the panel is a set of singularities, so keep a shell around it ---
+T_panel = np.where(d_min < 0.15, np.nan, T_panel)
+print(f"{N_sub} sub-sources, each {P_heat/N_sub:.3f} W")
+print(f"T ranges {np.nanmin(T_panel):.2f} to {np.nanmax(T_panel):.2f} K above ambient")
+
+# --- self-check (leave this alone) ---
+fw.check_shape("T has the shape of the grid", T_panel, X.shape)
+_i1 = int(np.argmin(np.abs(axis - 1.0)))
+fw.check_scalar("1 m directly above the centre of the panel", T_panel[c, c, _i1],
+                5.007, rtol=0.01, unit=" K")
+fw.check("...and the prefactor was applied, not left out",
+         np.nanmax(T_panel) < 100.0)
+```
+
+Heat flows down the temperature gradient, with the same minus sign and the same reason as $\boldsymbol{E} = -\nabla V$. Fourier's law is
+
+$$ \boldsymbol{q}_T = -k\nabla T \qquad [\text{W m}^{-2}] $$
+
+and at steady state, away from the panel, no heat is created or destroyed, so $\boldsymbol{q}_T$ should be solenoidal there. It is also a gradient field, so its curl should vanish.
+
+```{code-cell} ipython3
+# B2 -- three blanks. Reuse the operators you wrote: `divergence` from Task 1
+# and `curl` from Task 6. Both need arrays without NaN, so pass them through
+# np.nan_to_num first, as Task 2 did for the dipole.
+qx, qy, qz = ___                          # Fourier's law, as three arrays
+q_T = tuple(np.nan_to_num(v) for v in (qx, qy, qz))
+div_q = ___
+curl_q = ___
+
+# --- given: both reported scale-free, against |q|/d, exactly as Task 2 and
+#     Task 7 did. `far` is the region well clear of the panel.
+far = interior & (d_min > 0.6)
+q_mag = np.sqrt(q_T[0]**2 + q_T[1]**2 + q_T[2]**2)
+yardstick = (q_mag / np.maximum(d_min, 1e-12))[far]
+curl_mag = np.sqrt(curl_q[0]**2 + curl_q[1]**2 + curl_q[2]**2)
+print(f"  |div q| / (|q|/d), away from the panel : "
+      f"{np.median(np.abs(div_q[far]) / yardstick):.3%}")
+print(f"  |curl q| / (|q|/d)                     : "
+      f"{np.median(curl_mag[far] / yardstick):.2e}")
+
+# --- self-check (leave this alone) ---
+fw.check(f"q points away from the panel, so heat flows outward "
+         f"({np.median((q_T[0]*X + q_T[1]*Y + q_T[2]*Z)[far]):+.3f})",
+         np.median((q_T[0]*X + q_T[1]*Y + q_T[2]*Z)[far]) > 0)
+fw.check(f"no heat is created away from the panel "
+         f"({np.median(np.abs(div_q[far]) / yardstick):.2%})",
+         np.median(np.abs(div_q[far]) / yardstick) < 0.05)
+fw.check(f"and the flux of a gradient cannot circulate "
+         f"({np.median(curl_mag[far] / yardstick):.1e})",
+         np.median(curl_mag[far] / yardstick) < 1e-10)
+```
+
+The divergence is zero away from the panel and the panel is certainly a source, so the differential form has nothing to say about how strong it is. Put a closed surface around it instead. `closed_box_flux` from Task 5 works unchanged.
+
+```{code-cell} ipython3
+# --- given: the divergence theorem as an instrument, reading in watts ---
+print("  box half-width      power it finds")
+for h in (1.0, 1.4):
+    print(f"     {h:.1f} m            {closed_box_flux(*q_T, h):8.3f} W")
+print(f"\n     actually buried  {P_heat:8.3f} W")
+
+# The same box at h = 0.6 m returns 84.3 W. Its faces pass 0.1 m from the edge
+# of the panel, inside the shell that was masked out above; np.nan_to_num then
+# integrated the deleted samples as zeros. With no mask it returns 100.2 W.
+# Question B3 asks what the general rule is.
+
+# --- self-check (leave this alone) ---
+fw.check_scalar("closed-surface flux of q = the power buried inside",
+                closed_box_flux(*q_T, 1.0), P_heat, rtol=0.01, unit=" W")
+fw.check("...and a larger box finds the same power, not more",
+         abs(closed_box_flux(*q_T, 1.4) - closed_box_flux(*q_T, 1.0)) < 0.01 * P_heat)
+```
+
+Far from the panel its shape should stop mattering. Test that against the single term a point source would give.
+
+```{code-cell} ipython3
+# --- given: the panel against one point source of the same total power ---
+print(f"  {'distance':>10} {'along x':>10} {'along z':>10} {'point source':>14} {'spread':>8}")
+for d_ in (0.8, 1.2, 1.6):
+    i = int(np.argmin(np.abs(axis - d_)))
+    T_x, T_z = T_panel[i, c, c], T_panel[c, c, i]
+    print(f"  {d_:8.1f} m {T_x:9.3f} K {T_z:9.3f} K "
+          f"{P_heat/(4*np.pi*k_soil*d_):13.3f} K {abs(T_x-T_z)/T_x:8.1%}")
+
+# The blank shell in both panels is the masked region, 0.15 m around the
+# panel. Its outline in the plan view is the panel's own shape.
+fig, axes = plt.subplots(1, 2, figsize=(12.5, 4.6))
+fw.show_field_slice(X, Y, Z, q_T[0], q_T[2], background=T_panel, ax=axes[0],
+                    plane="y", density=1.2, vmin=0, vmax=15, levels=16,
+                    cmap="inferno", symmetric=False, stream_color="w",
+                    label="$T$ above ambient  [K]", title="vertical section, $y = 0$")
+fw.show_field_slice(X, Y, Z, q_T[0], q_T[1], background=T_panel, ax=axes[1],
+                    plane="z", density=1.2, vmin=0, vmax=15, levels=16,
+                    cmap="inferno", symmetric=False, stream_color="w",
+                    label="$T$ above ambient  [K]", title="plan view, $z = 0$")
+plt.tight_layout()
+plt.show()
+```
+
+:::{admonition} B3. Four questions on what you just measured
+:class: tip
+
+Answer these in writing.
+
+1. The closed surface returned 100.05 W and the divergence returned zero everywhere you could measure it. Both are correct. What does each one tell you that the other cannot?
+2. The $h = 0.6$ m box returns 84.3 W with the mask in place and 100.2 W without it. State the general rule this illustrates about masked samples and surface integrals.
+3. `curl_q` came back at $10^{-15}$ rather than at the fraction of a percent `div_q` shows. Why is it so much smaller, and is that a better measurement or a different kind of statement?
+4. In the plan view the isotherms near the panel are rounded rectangles and far away they are circles, and the table shows the difference between the two directions falling from 19% to 5%. What has been lost, and what does that have to do with truncating a series?
+:::
+
+:::{admonition} B4. The number that is wrong
+:class: tip
+
+Everything above was computed in soil, $k = 1.5$ W m⁻¹K⁻¹, and one metre above the panel it predicts $+5.0$ K. Re-run it for the same panel hanging in **air**, $k_{\text{air}} = 0.026$ W m⁻¹K⁻¹. You do not need to recompute anything: $T \propto 1/k$, so the answer is $5.0 \times 1.5/0.026$.
+
+The arithmetic is right and the answer is absurd. Identify the assumption that failed. Two are worth naming.
+:::
diff --git a/book/1_gradient_divergence_curl/sums_series_approx.md b/book/1_gradient_divergence_curl/sums_series_approx.md
new file mode 100644
index 0000000..3306266
--- /dev/null
+++ b/book/1_gradient_divergence_curl/sums_series_approx.md
@@ -0,0 +1,108 @@
+# Sums, series, approximations
+
+In all classical physics domains, it is useful to write a quantity as a sum of a large number of terms. In many cases, each such term has a physical interpretation. Often, the large sum of terms can be approximated by neglecting most terms based on physical arguments, and keeping only a few terms as the approximate solution for the quantity we investigate.
+
+A simple example is position as a function of time. A particle that is not moving has a fixed position and we can denote it as
+
+$$
+x(t) = x_0 .
+$$
+
+This expression states that for any time value the particle is located at position $x_0$. If at $t=0$ the particle starts to move with a constant velocity $v_0$, the position becomes a linear function of time. We can express it as
+
+$$
+x(t) = x_0 + v_0 t .
+$$
+
+If we take $t=0$ in this expression, we find the starting position $x(t=0) = x_0$. To find the velocity, we must differentiate the position with respect to time and find
+
+$$
+v_0 = \frac{\mathrm{d}x(t)}{\mathrm{d}t} .
+$$
+
+As long as the velocity is constant, we can evaluate this expression at any time instant, but if the velocity is variable we must evaluate the expression at $t=0$. Hence, we express velocity as
+
+$$
+v_0 = \lim_{t\downarrow 0}\frac{\mathrm{d}x(t)}{\mathrm{d}t}
+    = \left.\frac{\mathrm{d}x(t)}{\mathrm{d}t}\right|_{t\downarrow 0} .
+$$
+
+:::{admonition} The unit-step function
+:class: note
+Note that in this expression we take the limit from positive values of $t$ to zero, because the derivative of the position of the particle is not continuous at $t=0$. This can be seen because $x(t) = x_0$ for $t<0$, and the velocity of the particle is $v=0$ for $t<0$ and $v=v_0$ for $t>0$. We express this as
+
+$$
+v(t) = v_0\, u(t) ,
+$$
+
+where $u(t)$ is known as the unit-step function, given by
+
+$$
+u(t) = \left\{
+\begin{array}{ll}
+0    & t < 0 \\
+1/2  & t = 0 \\
+1    & t > 0
+\end{array}\right. .
+$$
+
+This function is also known as the Heaviside function, after Oliver Heaviside. The value at $t=0$ for $u(t)$ is obtained by taking the limit on both sides to $t=0$ and keeping the average. This is called the principal value. To avoid differentiating a function with a step discontinuity at this moment, we have used the differentiation in the time window where the position of the particle is continuous and continuously differentiable. We deal with differentiating across discontinuities later.
+:::
+
+Suppose the particle also has a constant acceleration $a_0$. From classical mechanics we know that the position of the particle as a function of time can then be expressed as
+
+$$
+x(t) = x_0 + v_0 t + \tfrac{1}{2}a_0 t^2 .
+$$
+
+To find $x_0$ and $v_0$ we can use the recipes above, while to obtain $a_0$ we must differentiate $x(t)$ twice with respect to $t$. Now we assume the position of the particle is twice continuously differentiable with respect to time, and this is true for both negative and positive times but not for $t=0$. Hence, to find the acceleration, we should evaluate
+
+$$
+a_0 = \lim_{t\downarrow 0}\frac{\mathrm{d}^2 x(t)}{\mathrm{d}t^2}
+    = \left.\frac{\mathrm{d}^2 x(t)}{\mathrm{d}t^2}\right|_{t\downarrow 0} .
+$$
+
+It shows the intuitive knowledge that acceleration is in the second derivative of position with respect to time.
+
+Now, if we generalise this notion, we can think of a position that changes location in an arbitrary way and has many more non-zero derivatives that are all smooth functions of time except across the start of the motion. This is one of the aspects of the principle of causality: a response cannot be present before an action happens. In classical mechanics it means the position cannot change unless the particle is already in motion or a force acts on it, in which case it has a non-zero acceleration. This notion belongs to the concept of the generation of fields and waves and we discuss it in more detail later. Here we state that the position of the particle can be generally expressed as
+
+$$
+x(t) = c_0 + c_1 t + c_2 t^2 + \cdots = \sum_{m=0}^{\infty} c_m t^m ,
+$$ (eq:pm)
+
+and
+
+$$
+c_m = \left.\frac{1}{m!}\frac{\mathrm{d}^m x(t)}{\mathrm{d}t^m}\right|_{t\downarrow 0} .
+$$
+
+## The Taylor series
+
+This series, where a function is expressed as a sum of terms with increasing powers of the independent variable, is called a Taylor series. It demonstrates that we can know a function if we know its value at every point in time, $x(t)$, **or** when we know all of its derivatives at one single time instant. This is under the assumptions that
+
+1. the function is continuously differentiable infinitely many times, and
+2. the series sums up to a finite result that represents the function.
+
+The Taylor series can of course be used for any function, for any variable, and in more than one dimension. We can therefore write an arbitrary function of position $x$ as $f(x)$ and express it as
+
+$$
+f(x) = \sum_{m=0}^{\infty}\left.\frac{1}{m!}\frac{\mathrm{d}^m f(x)}{\mathrm{d}x^m}\right|_{x=0} x^m
+     = f(0) + x f^{(1)}(0) + \frac{f^{(2)}(0)}{2}x^2 + \cdots ,
+$$
+
+where $f^{(m)}(0)$ is short-hand notation for $\left.\dfrac{\mathrm{d}^m f(x)}{\mathrm{d}x^m}\right|_{x=0}$.
+
+It is not necessary to expand a function around zero, and we can expand it around any point $x=a$. It is given by
+
+$$
+f(x) = \sum_{m=0}^{\infty}\frac{f^{(m)}(x=a)}{m!}(x-a)^m
+     = f(a) + f^{(1)}(a)(x-a) + \frac{f^{(2)}(a)}{2}(x-a)^2 + \cdots .
+$$ (eq:Taylor)
+
+We saw for the particle motion that for negative $t$ the particle is at rest and has position $x_0$. The reason is that at $t=0$ something happens and the information is not present in $x(t)$ for negative times. We saw that $x(t)$ changes continuously, but the slope does not change continuously. Functions that have a discontinuity in one or more derivatives are called non-analytic. Taylor series cannot be used for such functions. Until you reach a derivative that is not continuous, a truncated Taylor series expansion can still be useful.
+
+## Exercises
+
+1. Expand the particle motion of {eq}`eq:pm` around $t=-1$, using $f(x) = x(t)$ in {eq}`eq:Taylor` with $a=-1$, and explain why it is not giving you more than $x(t) = x_0$.
+2. Find the expansions for $\sin(x)$, $\cos(x)$, $\exp(-x)$, $(1+x)^{-1}$ around $x=0$ and determine whether the Taylor series converges.
+3. What happens if you try a Taylor series expansion for $\sqrt{t}$ and $1/t$?
diff --git a/book/_config.yml b/book/_config.yml
index 2044f6d..a60db27 100644
--- a/book/_config.yml
+++ b/book/_config.yml
@@ -2,7 +2,7 @@
 
 # This config includes an opinionated list of configuration options of TeachBooks, Jupyterbook v1, Sphinx and other extensions part of TeachBooks-Favourites.
 
-author: Author from Delft University of Technology, built with <a href="https://teachbooks.io">TeachBooks</a>, <a href="https://creativecommons.org/licenses/by/4.0/" target="_blank" rel="license noopener noreferrer"">CC BY 4.0<img style="height:22px!important;margin-left:3px;vertical-align:text-bottom;" src="https://mirrors.creativecommons.org/presskit/icons/cc.svg" alt=""><img style="height:22px!important;margin-left:3px;vertical-align:text-bottom;" src="https://mirrors.creativecommons.org/presskit/icons/by.svg" alt=""></a>
+author: Paco Lopez Dekker, Evert Slob, Guy Drijkoningen and Jinqiang Chen from Delft University of Technology, built with <a href="https://teachbooks.io">TeachBooks</a>, <a href="https://creativecommons.org/licenses/by/4.0/" target="_blank" rel="license noopener noreferrer"">CC BY 4.0<img style="height:22px!important;margin-left:3px;vertical-align:text-bottom;" src="https://mirrors.creativecommons.org/presskit/icons/cc.svg" alt=""><img style="height:22px!important;margin-left:3px;vertical-align:text-bottom;" src="https://mirrors.creativecommons.org/presskit/icons/by.svg" alt=""></a>
                                                      # Replace TeachBooks Teams with your own name in the line above, configuratin part of JupyterBook: https://jupyterbook.org/en/stable/customize/config.html
 execute:
   execute_notebooks: "auto"                          # Execute notebooks during build when outputs are missing or stale, configuration part of Jupyterbook v1 (https://jupyterbook.org/en/stable/customize/config.html)
@@ -27,13 +27,13 @@ sphinx:                                              # Options passed on to the
       use_thebe_lite: true                           # Required for live code, part of TeachBooks-Sphinx-Thebe (https://teachbooks.io/manual/features/live_code.html)
       exclude_patterns: ["**/_*.yml", "**/*.md", "**/*.ipynb"] #exclude files which should not be accessible for live code, configuration part of TeachBooks-Sphinx-Thebe (https://teachbooks.io/manual/features/live_code.html)
     #html_favicon:                                   # Default value not shown (line can be removed), allows to add your own favicon (disabling of TU Delft favicon is required), configuration part of Sphinx (https://www.sphinx-doc.org/en/master/index.html)
-    html_baseurl: "https://oit.tudelft.nl>/<repo>"   # Replace this with your own URL, configuration part of Sphinx (https://www.sphinx-doc.org/en/master/index.html)
+    html_baseurl: "https://bscect.github.io/FieldWaves/main/" # Replace this with your own URL, configuration part of Sphinx (https://www.sphinx-doc.org/en/master/index.html)
     html_theme_options:
       logo:
-        text: TU Delft GitHub OILM Template          # Replace this with your own open interactive learning material title, configuration part of Sphinx (https://www.sphinx-doc.org/en/master/index.html)
+        text: Fields and Waves                        # Replace this with your own open interactive learning material title, configuration part of Sphinx (https://www.sphinx-doc.org/en/master/index.html)
         # image_light:                               # Default value not shown (line can be removed), adds your logo for the light mode here (can be the same as image_dark) (disabling of TU Delft logo is required), configuration part of Sphinx (https://www.sphinx-doc.org/en/master/index.html)
         # image_dark:                                # Default value not shown (line can be removed), adds your logo for the dark mode here (can be the same as image_light) (disabling of TU Delft logo is required), configuration part of Sphinx (https://www.sphinx-doc.org/en/master/index.html)
-      repository_url: "https://github.com/TUDelft-books/repo" # Add your own repo URL here, configuration part of JupyterBook v1 (https://jupyterbook.org/en/stable/customize/config.html)
+      repository_url: "https://github.com/BScECT/FieldWaves" # Add your own repo URL here, configuration part of JupyterBook v1 (https://jupyterbook.org/en/stable/customize/config.html)
       path_to_docs:          "book"                  # Required for edit_page_button, should be book if you're using TeachBooks package (https://github.com/TeachBooks/TeachBooks) or the TeachBooks deploy book workflow (https://teachbooks.io/manual/external/deploy-book-workflow/README.html), configuration part of Jupyterbook v1 (https://jupyterbook.org/en/stable/customize/config.html)
       repository_branch:     "main"                  # Replace when you change the name of your published branch (required for edit_page_button), configuration part of Jupyterbook v1 (https://jupyterbook.org/en/stable/customize/config.html)
       use_edit_page_button:  true                    # Replace with false if you don't want the edit page button (i.e. if you don't user to propose changes themselves or have a private repo), configuration part of Jupyterbook v1 (https://jupyterbook.org/en/stable/customize/config.html)
diff --git a/book/_toc.yml b/book/_toc.yml
index e807a50..ae6c604 100644
--- a/book/_toc.yml
+++ b/book/_toc.yml
@@ -11,9 +11,13 @@ parts:
   - caption: Gradient, divergence, and curl
     chapters:
     - file: 1_gradient_divergence_curl/intro.md
-    # - file: 1_gradient_divergence_curl/gradient.md
-    # - file: 1_gradient_divergence_curl/divergence.md
-    # - file: 1_gradient_divergence_curl/curl.md
+    - file: 1_gradient_divergence_curl/sums_series_approx.md
+    - file: 1_gradient_divergence_curl/coord_sys.md
+    - file: 1_gradient_divergence_curl/gradient.md
+    - file: 1_gradient_divergence_curl/divergence.md
+    - file: 1_gradient_divergence_curl/curl.md
+    - file: 1_gradient_divergence_curl/labs/week01_series_grad.md
+    - file: 1_gradient_divergence_curl/labs/week02_div_curl.md
   - caption: Potential Fields
     chapters:
     - file: 2_potential_fields/introduction/intro.md
diff --git a/book/changelog.md b/book/changelog.md
index b9de8fa..8911db5 100644
--- a/book/changelog.md
+++ b/book/changelog.md
@@ -1,11 +1,6 @@
 # Changelog
 
-## `<latest version>`: `<date>`
-- `<Added/modified/deleted>` [](`<relative link to changed file>`)
-- ...
-- Full Changelog: `[<previous version>...<current version>](<link to diff as provided by GitHub>)`
-
-## `<previous version>`: <...>
-- <...>
-
-<...>
+## 2026-09-03
+- Added the Chapter 1 lecture notes: introduction, sums and series, coordinate systems, gradient, divergence and curl.
+- Added the two Chapter 1 computer labs, on series and the gradient, and on divergence and curl, with the shared `fwtools` helper module.
+- Set the book title, the repository link and the published URL.
diff --git a/book/credits.md b/book/credits.md
index 2168738..924bdc7 100644
--- a/book/credits.md
+++ b/book/credits.md
@@ -3,16 +3,16 @@
 
 You can refer to this book as:
 
-> `<editors>` from Delft University of Technology (`<year>`) _`<title>`_. `<url to book website>`. Source files at `<link to github repo`. CC BY 4.0.
+> Lopez Dekker, P., Slob, E., Drijkoningen, G. and Chen, J. from Delft University of Technology (2026) _Fields and Waves_. <https://bscect.github.io/FieldWaves/main/>. Source files at <https://github.com/BScECT/FieldWaves>. CC BY 4.0.
 
 You can refer to individual chapters or pages within this book as:
 
-> `<Title of Chapter or Page>`. In `<editors>` from Delft University of Technology (`<year>`) _`<title>`_. `<url to specific page on book website>`. Source files at `<link to specific commit / file in github repo`. CC BY 4.0.
+> `<Title of Chapter or Page>`. In Lopez Dekker, P., Slob, E., Drijkoningen, G. and Chen, J. from Delft University of Technology (2026) _Fields and Waves_. `<url to specific page on book website>`. Source files at `<link to specific commit / file in github repo>`. CC BY 4.0.
 
 We anticipate that the content of this book will change significantly. Therefore, we recommend using the source code directly with the citation above that refers to the GitHub repository and lists the date and name of the file. Although content will be added over time, chapter titles and URL's in this book are expected to remain relatively static. However, we make no guarantee, so if it is important for you to reference a specific location/commit within the book.
 
 ## How the book is made
-This website is written in markdown and jupyter notebooks files, which are converted to html using tools from [TeachBooks](https://teachbooks.io/). The files are stored on a [public GitHub repository](`<link to GitHub repo>`). The website can be viewed at `<link to book website url>`.
+This website is written in markdown and jupyter notebooks files, which are converted to html using tools from [TeachBooks](https://teachbooks.io/). The files are stored on a [public GitHub repository](https://github.com/BScECT/FieldWaves). The website can be viewed at <https://bscect.github.io/FieldWaves/main/>.
 
 To recreate the website you have two options (more information in the [TeachBooks manual](https://teachbooks.io/manual/):
 - In the GitHub interface: fork this repository, enable Github Pages from the source GitHub actions (Settings - Code and automation - Pages - Build and deployment - Source - GitHub Actions), enable workflows (Actions - I understand my workflows, go ahead and enable them) and run the call-deploy-book workflow (Actions - call-deploy-book - Run workflow - Run workflow). The website is released on the URL as shown on the workflow summary when the workflow has finished (Actions - call-deploy-book - call-deploy-book - Summary).
@@ -26,7 +26,7 @@ This book is [CC BY 4.0 licensed](https://creativecommons.org/licenses/by/4.0/)
 
 Parts of this book are taken from other external resources and reused in various ways. If an author is not listed on a particular page, it is by the Authors, except as follows:
 
-The following pages are included directly from an external resource and is not edited by `<Editor>`:
+The following pages are included directly from an external resource and are not edited by the course team:
 - The following pages includes text from {cite:t}`template`. Original content licensed under CC BY 4.0 License:
   - [](./exercises.md)
   - [](./exercises/002.md)
@@ -42,7 +42,7 @@ The following pages are included directly from an external resource and is not e
   - [](./syntax_exercises/011.md)
   - [](./exercises/summary.md) 
 
-The following pages contain content written by others, part of has been reused and/or modified by `<Editor>`
+The following pages contain content written by others, part of which has been reused and/or modified by the course team:
 - Page [](./exercises/001.md) includes text from {cite:t}`template` and is edited to be made TU Delft-specific. Original content licensed under CC BY 4.0 License. 
 - Page [](./syntax_exercises/012.md) includes text from {cite:t}`template` and is edited to be made TU Delft-specific. Original content licensed under CC BY 4.0 License. 
 
@@ -50,4 +50,11 @@ The following pages contain content written by others, part of has been reused a
 (editor)=
 ## About the Editors
 
+This book is written and maintained by the ECTB2140 teaching team at Delft University of Technology:
+
+- **Paco Lopez Dekker**, Department of Space Engineering, Faculty of Aerospace Engineering
+- **Evert Slob**, Department of Geoscience and Engineering, Faculty of Civil Engineering and Geosciences
+- **Guy Drijkoningen**, Department of Geoscience and Engineering, Faculty of Civil Engineering and Geosciences
+- **Jinqiang Chen**, Department of Geoscience and Engineering, Faculty of Civil Engineering and Geosciences
+
 ### Acknowledgements
diff --git a/book/intro.md b/book/intro.md
index 1b99fce..400ed64 100644
--- a/book/intro.md
+++ b/book/intro.md
@@ -4,12 +4,28 @@ title: ECTB2140 Fields and Waves 2026
 
 # Fields and Waves
 
+:::{admonition} This book is under development
+:class: warning
+
+The book is being written during the 2026-2027 run of ECTB2140, and the material is incomplete beyond the first chapter.
+
+- **Chapter 1, Gradient, divergence and curl**, is complete: the lecture notes and both computer labs.
+- **Chapter 2, Potential fields**, is partly written.
+- Later chapters are not yet available here.
+- The **Exercises** section in the sidebar is left over from the book template. It explains how TeachBooks works and is not course material.
+
+Brightspace remains the authoritative source for the schedule, assessment and announcements. Pages here change between weeks, so download a notebook again rather than relying on a copy saved earlier.
+:::
+
 ## Course description
 
 Physical objects can be described in a quantitative way only with the aid of mathematics. In classical physics, all physical objects are geometric objects. The course “Fields and Waves” combines the physics of fields and waves with the mathematical tools required to describe these phenomena. The course enables students to develop the essential understanding of the physical interpretation of mathematical formulations. Examples are radar waves that are used for earth surface and subsurface observations from antennas placed on satellites, airplanes, and close to or on the ground surface; sound waves, electromagnetic diffusion fields, electric and magnetic potential fields, and the gravity field to probe the earth’s interior. The approach is to start with relevant physical problems (1D wave equation, stretched membrane, electrostatic fields, heat diffusion, 2D and 3D waves) introducing mathematical tools and concepts (Partial Differential Equations (elliptic, parabolic, hyperbolic), separation of variables, Fourier series and Fourier transformations) during the course as needed to solve specific problems. While being able to find solutions for some specific cases is an important part of the course, the emphasis of the course is on learning how to describe the relevant physics with the corresponding math, to understand the solution space, and to be able to interpret solutions.
 
-### Expected prior knowledge:
-...
+### Expected prior knowledge
+
+- Calculus
+- Mechanics and Thermodynamics
+- Linear Algebra
 
 ## What to expect?
 <!-- This Jupyter book provides the road map for the course. All activities are introduced in the following chapters and there are a number of assignments you need to do. The didactic approach is based on guided tutorials where you study and practice the material yourself. The amount of lectures is kept to a minimum. In order to benefit most from this style of learning we urge you to prepare each lecture. We have scheduled time for this preparation in the time planner.
@@ -20,4 +36,4 @@ Brightspace will be used for providing day to day information relevant to this c
 
 
 
-<a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-sa/4.0/88x31.png" /></a>
\ No newline at end of file
+<a rel="license" href="https://creativecommons.org/licenses/by/4.0/"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by/4.0/88x31.png" /></a>
\ No newline at end of file