diff --git a/docs/conf.py b/docs/conf.py
index 366b53f..d7180f2 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -58,8 +58,8 @@
 
 # General information about the project.
 project = u'fmm2d'
-copyright = u'2018-2019 The Simons Foundation, Inc. - All Rights Reserved'
-author = u"Zydrunas Gimbutas, Leslie Greengard, Mike O'Neil, Manas Rachh, and Vladimir Rokhlin"
+copyright = u'2018-2023 The Simons Foundation, Inc. - All Rights Reserved'
+author = u"Travis Askham, Zydrunas Gimbutas, Leslie Greengard, Libin Lu, Mike O'Neil, Manas Rachh, and Vladimir Rokhlin"
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
diff --git a/docs/math.rst b/docs/math.rst
index d339bda..944c8fb 100644
--- a/docs/math.rst
+++ b/docs/math.rst
@@ -1,22 +1,23 @@
 Definitions 
 ===========
 Let $x_{j} = (x_{j,1}, x_{j,2}) \in \mathbb{R}^{2}$, $i=1,2,\ldots N$, denote a collection
-of source locations and let $t_{i} = (t_{i,1}, t_{i,2}) \in \mathbb{R}^{2}$ denote a collection
-of target locations. 
-
+of source locations and let $x = (x_{1}, x_{2}) \in \mathbb{R}^{2}$ denote a
+target location. Unless stated otherwise, the gradients are with
+respect to the target variable $x$.
 
 Laplace FMM
 *************
 The Laplace FMM comes in three varieties
 
 * rfmm2d: Charges, dipole strengths, potentials, their
-  gradients, and hessians are all double precision 
+  gradients, and Hessians are all double precision 
 * lfmm2d: All of the above quantities are double complex
 * cfmm2d: Same as lfmm2d except that there are no dipole orientation vectors,
-  and the dipole kernel is given by $1/(z_{i} - \xi_{j})$ where $z_{i}, \xi_{j}$
-  are the complex versions of the source and target locations given by
-  $z_{i} = t_{i,1} + i\cdot t_{i,2}$, and $\xi_{j} = x_{j,1} + i \cdot x_{j,2}$. 
-  Moreover, the gradients, and hessians
+  and the dipole kernel is given by $1/(z - \xi_{j})$ where $z, \xi_{j}$
+  are the target, and source locations viewed as points in the complex
+  plane given by
+  $z = x_{1} + i\cdot x_{2}$, and $\xi_{j} = x_{j,1} + i \cdot x_{j,2}$. 
+  Moreover, the gradients, and Hessians
   in this case are derivatives with respect to $z$, and in a slight abuse of
   notation, we set $\frac{\mathrm{d}}{\mathrm{d} z} \log(\|z\|) = 1/z$
  
@@ -31,8 +32,7 @@ denote a collection of dipole strengths, and $d_{j} \in \mathbb{R}^{2}$,
 $j=1,2,\ldots N$, denote the corresponding dipole orientation vectors.
 
 The Laplace FMM (rfmm2d) computes 
-the potential $u(x) \in \mathbb{R}$ and its gradient $\nabla u(x) \in
-\mathbb{R}^{2}$
+the potential $u(x) \in \mathbb{R}$
 given by
 
 .. math::
@@ -40,6 +40,7 @@ given by
 
     u(x) = \sum_{j=1}^{N} c_{j} \log{(\|x-x_{j}\|)} - v_{j} d_{j}\cdot \nabla \log{(\|x-x_{j}\|)}  \, , 
 
+and its gradient $\nabla u(x) \in \mathbb{R}^{2}$
 at the source and target locations. When $x=x_{j}$, the term
 corresponding to $x_{j}$ is dropped from the sum.
 
@@ -54,8 +55,7 @@ denote a collection of dipole strengths, and $d_{j} \in \mathbb{R}^{2}$,
 $j=1,2,\ldots N$, denote the corresponding dipole orientation vectors.
 
 The Laplace FMM (rfmm2d) computes 
-the potential $u(x) \in \mathbb{C}$ and its gradient $\nabla u(x) \in
-\mathbb{C}^{2}$
+the potential $u(x) \in \mathbb{C}$ 
 given by
 
 .. math::
@@ -63,6 +63,7 @@ given by
 
     u(x) = \sum_{j=1}^{N} c_{j} \log{(\|x-x_{j}\|)} - v_{j} d_{j}\cdot \nabla \log{(\|x-x_{j}\|)}  \, , 
 
+and its gradient $\nabla u(x) \in \mathbb{C}^{2}$
 at the source and target locations. When $x=x_{j}$, the term
 corresponding to $x_{j}$ is dropped from the sum.
 
@@ -76,8 +77,7 @@ $j=1,2,\ldots N$,
 denote a collection of dipole strengths.
 
 The Laplace FMM (cfmm2d) computes 
-the potential $u(z) \in \mathbb{C}$ and its gradient $\frac{\mathrm{d}}{\mathrm{d}z} u(z) \in
-\mathbb{C}$
+the potential $u(z) \in \mathbb{C}$ 
 given by
 
 .. math::
@@ -85,9 +85,11 @@ given by
 
     u(z) = \sum_{j=1}^{N} c_{j} \log{(\|z-\xi_{j}\|)} - \frac{v_{j}}{z-\xi_{j}}  \, , 
 
+and its gradient $\frac{\mathrm{d}}{\mathrm{d}z} u(z) \in
+\mathbb{C}$
 at the source and target locations. When $z=\xi_{j}$, the term
 corresponding to $\xi_{j}$ is dropped from the sum. Here
-$z = x(1) + i x(2)$, and $\xi_{j} = x_{j}(1) + ix_{j}(2)$, are the complex
+$z = x(1) + i x(2)$, and $\xi_{j} = x_{j,1} + ix_{j,2}$, are the complex
 numbers corresponding to $x$ and $x_{j}$ respectively.
 
 
@@ -103,7 +105,7 @@ Let $k\in\mathbb{C}$ denote the wave number or the Helmholtz
 parameter. 
 
 The Helmholtz FMM computes 
-the potential $u(x)$ and its gradient $\nabla u(x)$
+the potential $u(x) \in \mathbb{C}$ 
 given by
 
 .. math::
@@ -111,6 +113,7 @@ given by
 
     u(x) = \sum_{j=1}^{N} c_{j} H_{0}^{(1)}(k\|x-x_{j}\|) - v_{j} d_{j}\cdot \nabla H_{0}^{(1)}(k\|x-x_{j}\|)  \, , 
 
+and its gradient $\nabla u(x) \in \mathbb{C}^{2}$
 at the source and target locations, where $H_{0}^{(1)}$ is the Hankel function
 of the first kind of order $0$. When $x=x_{j}$, the term
 corresponding to $x_{j}$ is dropped from the sum.
@@ -126,7 +129,8 @@ $j=1,2,\ldots N$,
 denote a collection of dipole strengths.
 
 The Biharmonic FMM computes 
-the potential $u(x)$ and its `gradient` = $(P_{z} \frac{\mathrm{d}}{\mathrm{d}z}, P_{\overline{z}} \frac{\mathrm{d}}{\mathrm{d}z}, \frac{\mathrm{d}}{\mathrm{d}\overline{z}})$
+the potential $u(x)$ and its `gradient` = 
+$(P_{z} \frac{\mathrm{d}}{\mathrm{d}z}, P_{\overline{z}} \frac{\mathrm{d}}{\mathrm{d}z}, \frac{\mathrm{d}}{\mathrm{d}\overline{z}})$
 given by
 
 .. math::
@@ -148,12 +152,19 @@ given by
     \frac{v_{j,3}}{(\overline{z-\xi_{j}})^2} - 
     2v_{j,2} \frac{z - \xi_{j}}{(\overline{z-\xi_{j}})^3} \, , \\
 
-at the source and target locations. When $z=\xi_{j}$, the term
+at the source and target locations. 
+When $z=\xi_{j}$, the term
 corresponding to $\xi_{j}$ is dropped from the sum. Here
 $z = x(1) + i x(2)$, and $\xi_{j} = x_{j}(1) + ix_{j}(2)$, are the complex
 numbers corresponding to $x$ and $x_{j}$ respectively.
 When $x=x_{j}$, the term
 corresponding to $x_{j}$ is dropped from the sum.
+In the above, $(P_{z} \frac{\mathrm{d}}{\mathrm{d}z}$
+is the part of the gradient $\frac{\mathrm{d}}{\mathrm{d} z}$
+which is purely a function of $z$, and similarly
+$(P_{\overline{z}} \frac{\mathrm{d}}{\mathrm{d}z}$
+is the part of the gradient $\frac{\mathrm{d}}{\mathrm{d} z}$
+which is purely a function of $\overline{z}$.
 
 
 Modified Biharmonic FMM
@@ -169,23 +180,23 @@ where $K_{0}$ is the modified Bessel function of order $0$, and $\beta$ is the
 Modified Biharmonic wavenumber.
 
 
-Let $c_{j} \in \mathbb{C}$, 
+Let $c_{j} \in \mathbb{R}$, 
 denote a collection of charge strengths, 
-$v_{j} \in \mathbb{C}$, 
+$v_{j} \in \mathbb{R}$, 
 denote a collection of dipole strengths,
 $d_{j} = (d_{j,1}, d_{j,2})$ denote a collection
 of dipole vectors, 
-$q_{j} \in \mathbb{C}$ denote a collection of 
+$q_{j} \in \mathbb{R}$ denote a collection of 
 quadrupole strengths, 
-$w_{j} = (w_{j,1}, w_{j,2}, w_{j,3}) \in \mathbb{C}^{3}$, 
+$w_{j} = (w_{j,1}, w_{j,2}, w_{j,3}) \in \mathbb{R}^{3}$, 
 denote a collection of quadrupole orientation vectors, 
 $o_{j}$ denote a collection of octopole strengths, and
-$p_{j} = (p_{j,1}, p_{j,2}, p_{j,3}, p_{j,4}) \in \mathbb{C}^{4}$, 
-denote a collceetion of octopole strengths. For all the vectors,
+$p_{j} = (p_{j,1}, p_{j,2}, p_{j,3}, p_{j,4}) \in \mathbb{R}^{4}$, 
+denote a collection of octopole strengths. For all the vectors,
 $j=1,2,\ldots N$.
 
 The Modified Biharmonic FMM computes 
-the potential $u(x)$ and its gradients, and hessians, 
+the potential $u(x)\in \mathbb{R}$ 
 given by
 
 .. math::
@@ -201,6 +212,7 @@ given by
     &p_{j,3} \partial_{y_{1},y_{2},y_{2}} G^{\textrm{mbh}}(x,x_{j}) +
     p_{j,4} \partial_{y_{2},y_{2},y_{2}} G^{\textrm{mbh}}(x,x_{j}) \big) \, ,
 
+and its gradients $\nabla u(x) \in \mathbb{R}^{2}$
 at the source and target locations. When $x=x_{j}$, the term
 corresponding to $x_{j}$ is dropped from the sum.
 
@@ -236,10 +248,10 @@ is given by
     (x_{1}-y_{1})(x_{2}-y_{2}) & (x_{2} - y_{2})^2
     \end{bmatrix} \, ,
 
-and finally let $\Pi^{\textrm{stok}} (x,y)$ denote its assocaited pressure tensor given by
+and finally let $\Pi^{\textrm{stok}} (x,y)$ denote its associated pressure tensor given by
 
 .. math::
-    v\cdot \Pi(x,y)^{textrm{stok}} = -\frac{v}{\|x-y \|^2} + \frac{2 v \cdot(x-y)}{\|x-y \|^4}
+    v\cdot \Pi(x,y)^{\textrm{stok}} = -\frac{v}{\|x-y \|^2} + \frac{2 v \cdot(x-y)}{\|x-y \|^4}
     \begin{bmatrix}
     (x_{1}-y_{1}) \\
     (x_{2}-y_{2})
@@ -247,10 +259,10 @@ and finally let $\Pi^{\textrm{stok}} (x,y)$ denote its assocaited pressure tenso
 
 Let $c_{j} \in \mathbb{R}^2$, 
 $j=1,2,\ldots N$, 
-denote a collection of stokeslet strengths, $v_{j} \in \mathbb{R}^2$,
+denote a collection of Stokeslet strengths, $v_{j} \in \mathbb{R}^2$,
 $j=1,2,\ldots N$, 
-denote a collection of stresslet strengths, and $d_{j} \in \mathbb{R}^{2}$,
-$j=1,2,\ldots N$, denote the corresponding stresslet orientation vectors.
+denote a collection of Stresslet strengths, and $d_{j} \in \mathbb{R}^{2}$,
+$j=1,2,\ldots N$, denote the corresponding Stresslet orientation vectors.
 
 The Stokes FMM computes 
 the potential $u(x) \in \mathbb{R}^2$, its gradient $\nabla u(x) \in
@@ -269,24 +281,28 @@ corresponding to $x_{j}$ is dropped from the sum.
 
 Vectorized versions   
 *******************
-The vectorized versions of the Helmholtz FMM, 
-computes repeated FMMs for new charge and dipole strengths
-located at the same source locations, where the potential and its
-gradient are evaluated at the same set of target locations.
+The vectorized versions of the FMMs compute the same sums
+as above for a set of problems in which the source and target
+locations are constant but multiple values of the charges, dipoles, etc 
+are specified at each source. Given a set of problems with this structure,
+the vectorized versions are faster than calling the standard
+FMM multiple times in sequence. 
 
 For example, let $c_{\ell,j}\in\mathbb{C}$, 
 $j=1,2,\ldots N$, $\ell=1,2,\ldots n_{d}$
 denote a collection of $n_{d}$ charge strengths, and
 let $v_{\ell,j} \in \mathbb{C}$, $d_{\ell,j} \in \mathbb{R}^2$ 
 denote a collection of $n_{d}$ dipole strengths and orientation vectors. 
-Then the vectorized Helmholtz FMM computes the potentials $u_{\ell}(x)$ 
-and its gradients $\nabla u_{\ell}(x)$ defined by the formula
+Then the vectorized Helmholtz FMM computes the potentials $u_{\ell}(x) \in
+\mathbb{C}$
+given by
 
 .. math::
     :label: helm_nbody_vec
 
-    u_{\ell}(x) = \sum_{j=1}^{N} c_{\ell,j} H_{0}^{(1)}(k\|x-x_{j}\|) - v_{\ell,j} d_{j}\cdot \nabla H_{0}^{(1)}(k\|x-x_{j}\|)  \, , 
+    u_{\ell}(x) = \sum_{j=1}^{N} c_{\ell,j} H_{0}^{(1)}(k\|x-x_{j}\|) - v_{\ell,j} d_{\ell,j}\cdot \nabla H_{0}^{(1)}(k\|x-x_{j}\|)  \, , 
 
+and its gradients $\nabla u_{\ell}(x) \in \mathbb{C}^{2}$
 at the source and target locations. 
 
 .. note::