diff --git a/examples/fdtd.py b/examples/fdtd.py
index be3942b..3ffa077 100644
--- a/examples/fdtd.py
+++ b/examples/fdtd.py
@@ -20,9 +20,10 @@ def perturbed_l3(a: float, radius: float, **kwargs) -> Pattern:
     """
     Generate a masque.Pattern object containing a perturbed L3 cavity.
 
-    :param a: Lattice constant.
-    :param radius: Hole radius, in units of a (lattice constant).
-    :param kwargs: Keyword arguments:
+    Args:
+        a: Lattice constant.
+        radius: Hole radius, in units of a (lattice constant).
+        **kwargs: Keyword arguments:
         hole_dose, trench_dose, hole_layer, trench_layer: Shape properties for Pattern.
                 Defaults *_dose=1, hole_layer=0, trench_layer=1.
         shifts_a, shifts_r: passed to pcgen.l3_shift; specifies lattice constant (1 -
@@ -30,11 +31,13 @@ def perturbed_l3(a: float, radius: float, **kwargs) -> Pattern:
                 holes adjacent to the defect (same row). Defaults are 0.15 shift for
                 first hole, 0.075 shift for third hole, and no radius change.
         xy_size: [x, y] number of mirror periods in each direction; total size is
-                2 * n + 1 holes in each direction. Default [10, 10].
+                `2 * n + 1` holes in each direction. Default `[10, 10]`.
         perturbed_radius: radius of holes perturbed to form an upwards-driected beam
                 (multiplicative factor). Default 1.1.
         trench width: Width of the undercut trenches. Default 1.2e3.
-    :return: masque.Pattern object containing the L3 design
+
+    Return:
+        `masque.Pattern` object containing the L3 design
     """
 
     default_args = {'hole_dose':    1,
diff --git a/meanas/fdfd/bloch.py b/meanas/fdfd/bloch.py
index 28b82df..db055f6 100644
--- a/meanas/fdfd/bloch.py
+++ b/meanas/fdfd/bloch.py
@@ -4,47 +4,54 @@ Bloch eigenmode solver/operators
 This module contains functions for generating and solving the
  3D Bloch eigenproblem. The approach is to transform the problem
  into the (spatial) fourier domain, transforming the equation
-   1/mu * curl(1/eps * curl(H)) = (w/c)^2 H
+
+    1/mu * curl(1/eps * curl(H)) = (w/c)^2 H
+
  into
-   conv(1/mu_k, ik x conv(1/eps_k, ik x H_k)) = (w/c)^2 H_k
+
+    conv(1/mu_k, ik x conv(1/eps_k, ik x H_k)) = (w/c)^2 H_k
+
  where:
-  - the _k subscript denotes a 3D fourier transformed field
-  - each component of H_k corresponds to a plane wave with wavevector k
-  - x is the cross product
-  - conv denotes convolution
 
- Since k and H are orthogonal for each plane wave, we can use each
- k to create an orthogonal basis (k, m, n), with k x m = n, and
- |m| = |n| = 1. The cross products are then simplified with
+  - the `_k` subscript denotes a 3D fourier transformed field
+  - each component of `H_k` corresponds to a plane wave with wavevector `k`
+  - `x` is the cross product
+  - `conv()` denotes convolution
 
- k @ h = kx hx + ky hy + kz hz = 0 = hk
- h = hk + hm + hn = hm + hn
- k = kk + km + kn = kk = |k|
+ Since `k` and `H` are orthogonal for each plane wave, we can use each
+ `k` to create an orthogonal basis (k, m, n), with `k x m = n`, and
+ `|m| = |n| = 1`. The cross products are then simplified with
 
- k x h = (ky hz - kz hy,
-          kz hx - kx hz,
-          kx hy - ky hx)
-       = ((k x h) @ k, (k x h) @ m, (k x h) @ n)_kmn
-       = (0, (m x k) @ h, (n x k) @ h)_kmn         # triple product ordering
-       = (0, kk (-n @ h), kk (m @ h))_kmn          # (m x k) = -|k| n, etc.
-       = |k| (0, -h @ n, h @ m)_kmn
+    k @ h = kx hx + ky hy + kz hz = 0 = hk
+    h = hk + hm + hn = hm + hn
+    k = kk + km + kn = kk = |k|
 
- k x h = (km hn - kn hm,
-          kn hk - kk hn,
-          kk hm - km hk)_kmn
-       = (0, -kk hn, kk hm)_kmn
-       = (-kk hn)(mx, my, mz) + (kk hm)(nx, ny, nz)
-       = |k| (hm * (nx, ny, nz) - hn * (mx, my, mz))
+    k x h = (ky hz - kz hy,
+             kz hx - kx hz,
+             kx hy - ky hx)
+          = ((k x h) @ k, (k x h) @ m, (k x h) @ n)_kmn
+          = (0, (m x k) @ h, (n x k) @ h)_kmn         # triple product ordering
+          = (0, kk (-n @ h), kk (m @ h))_kmn          # (m x k) = -|k| n, etc.
+          = |k| (0, -h @ n, h @ m)_kmn
 
- where h is shorthand for H_k, (...)_kmn deontes the (k, m, n) basis,
- and e.g. hm is the component of h in the m direction.
+    k x h = (km hn - kn hm,
+             kn hk - kk hn,
+             kk hm - km hk)_kmn
+          = (0, -kk hn, kk hm)_kmn
+          = (-kk hn)(mx, my, mz) + (kk hm)(nx, ny, nz)
+          = |k| (hm * (nx, ny, nz) - hn * (mx, my, mz))
 
- We can also simplify conv(X_k, Y_k) as fftn(X * ifftn(Y_k)).
+ where `h` is shorthand for `H_k`, `(...)_kmn` deontes the `(k, m, n)` basis,
+ and e.g. `hm` is the component of `h` in the `m` direction.
+
+ We can also simplify `conv(X_k, Y_k)` as `fftn(X * ifftn(Y_k))`.
+
+ Using these results and storing `H_k` as `h = (hm, hn)`, we have
+
+    e_xyz = fftn(1/eps * ifftn(|k| (hm * n - hn * m)))
+    b_mn = |k| (-e_xyz @ n, e_xyz @ m)
+    h_mn = fftn(1/mu * ifftn(b_m * m + b_n * n))
 
- Using these results and storing H_k as h = (hm, hn), we have
-  e_xyz = fftn(1/eps * ifftn(|k| (hm * n - hn * m)))
-  b_mn = |k| (-e_xyz @ n, e_xyz @ m)
-  h_mn = fftn(1/mu * ifftn(b_m * m + b_n * n))
  which forms the operator from the left side of the equation.
 
  We can then use a preconditioned block Rayleigh iteration algorithm, as in
@@ -120,12 +127,15 @@ def generate_kmn(k0: numpy.ndarray,
     """
     Generate a (k, m, n) orthogonal basis for each k-vector in the simulation grid.
 
-    :param k0: [k0x, k0y, k0z], Bloch wavevector, in G basis.
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param shape: [nx, ny, nz] shape of the simulation grid.
-    :return: (|k|, m, n) where |k| has shape tuple(shape) + (1,)
-            and m, n have shape tuple(shape) + (3,).
-            All are given in the xyz basis (e.g. |k|[0,0,0] = norm(G_matrix @ k0)).
+    Args:
+        k0: [k0x, k0y, k0z], Bloch wavevector, in G basis.
+        G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
+        shape: [nx, ny, nz] shape of the simulation grid.
+
+    Returns:
+        `(|k|, m, n)` where `|k|` has shape `tuple(shape) + (1,)`
+            and `m`, `n` have shape `tuple(shape) + (3,)`.
+            All are given in the xyz basis (e.g. `|k|[0,0,0] = norm(G_matrix @ k0)`).
     """
     k0 = numpy.array(k0)
 
@@ -159,21 +169,27 @@ def maxwell_operator(k0: numpy.ndarray,
                      ) -> Callable[[numpy.ndarray], numpy.ndarray]:
     """
     Generate the Maxwell operator
+
         conv(1/mu_k, ik x conv(1/eps_k, ik x ___))
+
     which is the spatial-frequency-space representation of
+
         1/mu * curl(1/eps * curl(___))
 
-    The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size)
+    The operator is a function that acts on a vector h_mn of size `2 * epsilon[0].size`
 
-    See the module-level docstring for more information.
+    See the `meanas.fdfd.bloch` docstring for more information.
 
-    :param k0: Bloch wavevector, [k0x, k0y, k0z].
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param epsilon: Dielectric constant distribution for the simulation.
-        All fields are sampled at cell centers (i.e., NOT Yee-gridded)
-    :param mu: Magnetic permability distribution for the simulation.
-        Default None (1 everywhere).
-    :return: Function which applies the maxwell operator to h_mn.
+    Args:
+        k0: Bloch wavevector, `[k0x, k0y, k0z]`.
+        G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
+        epsilon: Dielectric constant distribution for the simulation.
+                 All fields are sampled at cell centers (i.e., NOT Yee-gridded)
+        mu: Magnetic permability distribution for the simulation.
+            Default None (1 everywhere).
+
+    Returns:
+        Function which applies the maxwell operator to h_mn.
     """
 
     shape = epsilon[0].shape + (1,)
@@ -189,8 +205,11 @@ def maxwell_operator(k0: numpy.ndarray,
 
         h is complex 2-field in (m, n) basis, vectorized
 
-        :param h: Raveled h_mn; size (2 * epsilon[0].size).
-        :return: Raveled conv(1/mu_k, ik x conv(1/eps_k, ik x h_mn)).
+        Args:
+            h: Raveled h_mn; size `2 * epsilon[0].size`.
+
+        Returns:
+            Raveled conv(1/mu_k, ik x conv(1/eps_k, ik x h_mn)).
         """
         hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
 
@@ -231,18 +250,22 @@ def hmn_2_exyz(k0: numpy.ndarray,
                ) -> Callable[[numpy.ndarray], fdfield_t]:
     """
     Generate an operator which converts a vectorized spatial-frequency-space
-     h_mn into an E-field distribution, i.e.
+     `h_mn` into an E-field distribution, i.e.
+
         ifft(conv(1/eps_k, ik x h_mn))
 
-    The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size)
+    The operator is a function that acts on a vector `h_mn` of size `2 * epsilon[0].size`.
 
-    See the module-level docstring for more information.
+    See the `meanas.fdfd.bloch` docstring for more information.
 
-    :param k0: Bloch wavevector, [k0x, k0y, k0z].
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param epsilon: Dielectric constant distribution for the simulation.
-        All fields are sampled at cell centers (i.e., NOT Yee-gridded)
-    :return: Function for converting h_mn into E_xyz
+    Args:
+        k0: Bloch wavevector, `[k0x, k0y, k0z]`.
+        G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
+        epsilon: Dielectric constant distribution for the simulation.
+                 All fields are sampled at cell centers (i.e., NOT Yee-gridded)
+
+    Returns:
+        Function for converting `h_mn` into `E_xyz`
     """
     shape = epsilon[0].shape + (1,)
     epsilon = numpy.stack(epsilon, 3)
@@ -266,18 +289,22 @@ def hmn_2_hxyz(k0: numpy.ndarray,
                ) -> Callable[[numpy.ndarray], fdfield_t]:
     """
     Generate an operator which converts a vectorized spatial-frequency-space
-     h_mn into an H-field distribution, i.e.
+     `h_mn` into an H-field distribution, i.e.
+
         ifft(h_mn)
 
-    The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size)
+    The operator is a function that acts on a vector `h_mn` of size `2 * epsilon[0].size`.
 
-    See the module-level docstring for more information.
+    See the `meanas.fdfd.bloch` docstring for more information.
 
-    :param k0: Bloch wavevector, [k0x, k0y, k0z].
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param epsilon: Dielectric constant distribution for the simulation.
-        Only epsilon[0].shape is used.
-    :return: Function for converting h_mn into H_xyz
+    Args:
+        k0: Bloch wavevector, `[k0x, k0y, k0z]`.
+        G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
+        epsilon: Dielectric constant distribution for the simulation.
+                 Only `epsilon[0].shape` is used.
+
+    Returns:
+        Function for converting `h_mn` into `H_xyz`
     """
     shape = epsilon[0].shape + (1,)
     _k_mag, m, n = generate_kmn(k0, G_matrix, shape)
@@ -298,18 +325,23 @@ def inverse_maxwell_operator_approx(k0: numpy.ndarray,
                                     ) -> Callable[[numpy.ndarray], numpy.ndarray]:
     """
     Generate an approximate inverse of the Maxwell operator,
+
         ik x conv(eps_k, ik x conv(mu_k, ___))
+
      which can be used to improve the speed of ARPACK in shift-invert mode.
 
-    See the module-level docstring for more information.
+    See the `meanas.fdfd.bloch` docstring for more information.
 
-    :param k0: Bloch wavevector, [k0x, k0y, k0z].
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param epsilon: Dielectric constant distribution for the simulation.
-        All fields are sampled at cell centers (i.e., NOT Yee-gridded)
-    :param mu: Magnetic permability distribution for the simulation.
-        Default None (1 everywhere).
-    :return: Function which applies the approximate inverse of the maxwell operator to h_mn.
+    Args:
+        k0: Bloch wavevector, `[k0x, k0y, k0z]`.
+        G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
+        epsilon: Dielectric constant distribution for the simulation.
+                 All fields are sampled at cell centers (i.e., NOT Yee-gridded)
+        mu: Magnetic permability distribution for the simulation.
+            Default None (1 everywhere).
+
+    Returns:
+        Function which applies the approximate inverse of the maxwell operator to `h_mn`.
     """
     shape = epsilon[0].shape + (1,)
     epsilon = numpy.stack(epsilon, 3)
@@ -325,8 +357,11 @@ def inverse_maxwell_operator_approx(k0: numpy.ndarray,
 
         h is complex 2-field in (m, n) basis, vectorized
 
-        :param h: Raveled h_mn; size (2 * epsilon[0].size).
-        :return: Raveled ik x conv(eps_k, ik x conv(mu_k, h_mn))
+        Args:
+            h: Raveled h_mn; size `2 * epsilon[0].size`.
+
+        Returns:
+            Raveled ik x conv(eps_k, ik x conv(mu_k, h_mn))
         """
         hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
 
@@ -376,16 +411,20 @@ def find_k(frequency: float,
     """
     Search for a bloch vector that has a given frequency.
 
-    :param frequency: Target frequency.
-    :param tolerance: Target frequency tolerance.
-    :param direction: k-vector direction to search along.
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param epsilon: Dielectric constant distribution for the simulation.
-        All fields are sampled at cell centers (i.e., NOT Yee-gridded)
-    :param mu: Magnetic permability distribution for the simulation.
-        Default None (1 everywhere).
-    :param band: Which band to search in. Default 0 (lowest frequency).
-    return: (k, actual_frequency) The found k-vector and its frequency
+    Args:
+        frequency: Target frequency.
+        tolerance: Target frequency tolerance.
+        direction: k-vector direction to search along.
+        G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
+        epsilon: Dielectric constant distribution for the simulation.
+                 All fields are sampled at cell centers (i.e., NOT Yee-gridded)
+        mu: Magnetic permability distribution for the simulation.
+            Default None (1 everywhere).
+        band: Which band to search in. Default 0 (lowest frequency).
+
+    Returns:
+        `(k, actual_frequency)`
+        The found k-vector and its frequency.
     """
 
     direction = numpy.array(direction) / norm(direction)
@@ -419,16 +458,19 @@ def eigsolve(num_modes: int,
     Find the first (lowest-frequency) num_modes eigenmodes with Bloch wavevector
      k0 of the specified structure.
 
-    :param k0: Bloch wavevector, [k0x, k0y, k0z].
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param epsilon: Dielectric constant distribution for the simulation.
-        All fields are sampled at cell centers (i.e., NOT Yee-gridded)
-    :param mu: Magnetic permability distribution for the simulation.
-        Default None (1 everywhere).
-    :param tolerance: Solver stops when fractional change in the objective
-        trace(Z.H @ A @ Z @ inv(Z Z.H)) is smaller than the tolerance
-    :return: (eigenvalues, eigenvectors) where eigenvalues[i] corresponds to the
-        vector eigenvectors[i, :]
+    Args:
+        k0: Bloch wavevector, `[k0x, k0y, k0z]`.
+        G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
+        epsilon: Dielectric constant distribution for the simulation.
+                 All fields are sampled at cell centers (i.e., NOT Yee-gridded)
+        mu: Magnetic permability distribution for the simulation.
+            Default `None` (1 everywhere).
+        tolerance: Solver stops when fractional change in the objective
+                   `trace(Z.H @ A @ Z @ inv(Z Z.H))` is smaller than the tolerance
+
+    Returns:
+        `(eigenvalues, eigenvectors)` where `eigenvalues[i]` corresponds to the
+        vector `eigenvectors[i, :]`
     """
     h_size = 2 * epsilon[0].size
 
diff --git a/meanas/fdfd/farfield.py b/meanas/fdfd/farfield.py
index d0ec008..ea11224 100644
--- a/meanas/fdfd/farfield.py
+++ b/meanas/fdfd/farfield.py
@@ -21,27 +21,31 @@ def near_to_farfield(E_near: fdfield_t,
 
     The input fields should be complex phasors.
 
-    :param E_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse
-        E fields (e.g. [Ex, Ey] for calculating the farfield toward the z-direction).
-    :param H_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse
-        H fields (e.g. [Hx, hy] for calculating the farfield towrad the z-direction).
-    :param dx: Cell size along x-dimension, in units of wavelength.
-    :param dy: Cell size along y-dimension, in units of wavelength.
-    :param padded_size: Shape of the output. A single integer `n` will be expanded to `(n, n)`.
-        Powers of 2 are most efficient for FFT computation.
-        Default is the smallest power of 2 larger than the input, for each axis.
-    :returns: Dict with keys
-        'E_far': Normalized E-field farfield; multiply by
-            (i k exp(-i k r) / (4 pi r)) to get the actual field value.
-        'H_far': Normalized H-field farfield; multiply by
-            (i k exp(-i k r) / (4 pi r)) to get the actual field value.
-        'kx', 'ky': Wavevector values corresponding to the x- and y- axes in E_far and H_far,
-            normalized to wavelength (dimensionless).
-        'dkx', 'dky': step size for kx and ky, normalized to wavelength.
-        'theta': arctan2(ky, kx) corresponding to each (kx, ky).
-            This is the angle in the x-y plane, counterclockwise from above, starting from +x.
-        'phi': arccos(kz / k) corresponding to each (kx, ky).
-            This is the angle away from +z.
+    Args:
+        E_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse
+                E fields (e.g. [Ex, Ey] for calculating the farfield toward the z-direction).
+        H_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse
+                H fields (e.g. [Hx, hy] for calculating the farfield towrad the z-direction).
+        dx: Cell size along x-dimension, in units of wavelength.
+        dy: Cell size along y-dimension, in units of wavelength.
+        padded_size: Shape of the output. A single integer `n` will be expanded to `(n, n)`.
+                     Powers of 2 are most efficient for FFT computation.
+                     Default is the smallest power of 2 larger than the input, for each axis.
+
+    Returns:
+        Dict with keys
+
+        -   `E_far`: Normalized E-field farfield; multiply by
+                (i k exp(-i k r) / (4 pi r)) to get the actual field value.
+        -   `H_far`: Normalized H-field farfield; multiply by
+                (i k exp(-i k r) / (4 pi r)) to get the actual field value.
+        -   `kx`, `ky`: Wavevector values corresponding to the x- and y- axes in E_far and H_far,
+                normalized to wavelength (dimensionless).
+        -   `dkx`, `dky`: step size for kx and ky, normalized to wavelength.
+        -   `theta`: arctan2(ky, kx) corresponding to each (kx, ky).
+                This is the angle in the x-y plane, counterclockwise from above, starting from +x.
+        -   `phi`: arccos(kz / k) corresponding to each (kx, ky).
+                This is the angle away from +z.
     """
 
     if not len(E_near) == 2:
@@ -129,23 +133,27 @@ def far_to_nearfield(E_far: fdfield_t,
 
     The input fields should be complex phasors.
 
-    :param E_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse
-        E fields (e.g. [Ex, Ey] for calculating the nearfield toward the z-direction).
-        Fields should be normalized so that
-            E_far = E_far_actual / (i k exp(-i k r) / (4 pi r))
-    :param H_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse
-        H fields (e.g. [Hx, hy] for calculating the nearfield toward the z-direction).
-        Fields should be normalized so that
-            H_far = H_far_actual / (i k exp(-i k r) / (4 pi r))
-    :param dkx: kx discretization, in units of wavelength.
-    :param dky: ky discretization, in units of wavelength.
-    :param padded_size: Shape of the output. A single integer `n` will be expanded to `(n, n)`.
-        Powers of 2 are most efficient for FFT computation.
-        Default is the smallest power of 2 larger than the input, for each axis.
-    :returns: Dict with keys
-        'E': E-field nearfield
-        'H': H-field nearfield
-        'dx', 'dy': spatial discretization, normalized to wavelength (dimensionless)
+    Args:
+        E_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse
+                E fields (e.g. [Ex, Ey] for calculating the nearfield toward the z-direction).
+                Fields should be normalized so that
+                E_far = E_far_actual / (i k exp(-i k r) / (4 pi r))
+        H_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse
+                H fields (e.g. [Hx, hy] for calculating the nearfield toward the z-direction).
+                Fields should be normalized so that
+                H_far = H_far_actual / (i k exp(-i k r) / (4 pi r))
+        dkx: kx discretization, in units of wavelength.
+        dky: ky discretization, in units of wavelength.
+        padded_size: Shape of the output. A single integer `n` will be expanded to `(n, n)`.
+                     Powers of 2 are most efficient for FFT computation.
+                     Default is the smallest power of 2 larger than the input, for each axis.
+
+    Returns:
+        Dict with keys
+
+        -   `E`: E-field nearfield
+        -   `H`: H-field nearfield
+        -   `dx`, `dy`: spatial discretization, normalized to wavelength (dimensionless)
     """
 
     if not len(E_far) == 2:
diff --git a/meanas/fdmath/vectorization.py b/meanas/fdmath/vectorization.py
index 8e8099b..63d78ef 100644
--- a/meanas/fdmath/vectorization.py
+++ b/meanas/fdmath/vectorization.py
@@ -1,6 +1,6 @@
 """
-Functions for moving between a vector field (list of 3 ndarrays, [f_x, f_y, f_z])
-and a 1D array representation of that field [f_x0, f_x1, f_x2,... f_y0,... f_z0,...].
+Functions for moving between a vector field (list of 3 ndarrays, `[f_x, f_y, f_z]`)
+and a 1D array representation of that field `[f_x0, f_x1, f_x2,... f_y0,... f_z0,...]`.
 Vectorized versions of the field use row-major (ie., C-style) ordering.
 """
 
@@ -17,11 +17,14 @@ def vec(f: fdfield_t) -> vfdfield_t:
     """
     Create a 1D ndarray from a 3D vector field which spans a 1-3D region.
 
-    Returns None if called with f=None.
+    Returns `None` if called with `f=None`.
 
-    :param f: A vector field, [f_x, f_y, f_z] where each f_ component is a 1 to
-        3D ndarray (f_* should all be the same size). Doesn't fail with f=None.
-    :return: A 1D ndarray containing the linearized field (or None)
+    Args:
+        f: A vector field, `[f_x, f_y, f_z]` where each `f_` component is a 1- to
+           3-D ndarray (`f_*` should all be the same size). Doesn't fail with `f=None`.
+
+    Returns:
+        1D ndarray containing the linearized field (or `None`)
     """
     if numpy.any(numpy.equal(f, None)):
         return None
@@ -31,15 +34,17 @@ def vec(f: fdfield_t) -> vfdfield_t:
 def unvec(v: vfdfield_t, shape: numpy.ndarray) -> fdfield_t:
     """
     Perform the inverse of vec(): take a 1D ndarray and output a 3D field
-     of form [f_x, f_y, f_z] where each of f_* is a len(shape)-dimensional
+     of form `[f_x, f_y, f_z]` where each of `f_*` is a len(shape)-dimensional
      ndarray.
 
-    Returns None if called with v=None.
+    Returns `None` if called with `v=None`.
 
-    :param v: 1D ndarray representing a 3D vector field of shape shape (or None)
-    :param shape: shape of the vector field
-    :return: [f_x, f_y, f_z] where each f_ is a len(shape) dimensional ndarray
-     (or None)
+    Args:
+        v: 1D ndarray representing a 3D vector field of shape shape (or None)
+        shape: shape of the vector field
+
+    Returns:
+        `[f_x, f_y, f_z]` where each `f_` is a `len(shape)` dimensional ndarray (or `None`)
     """
     if numpy.any(numpy.equal(v, None)):
         return None