updating documentation and adding python file

docs/math.rst

@@ -1,7 +1,7 @@
 Definitions 
 ===========
-Let $x_{j} \in \mathbb{R}^{2}$, $i=1,2,\ldots N$, denote a collection
-of source locations and let $t_{i} \in \mathbb{R}^{2}$ denote a collection
+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. 
 
 
@@ -13,12 +13,12 @@ The Laplace FMM comes in three varieties
   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/(zt_{i} - zx_{j})$ where $zt_{i}, zx_{j}$
+  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
-  $zt_{i} = t_{i}(1) + i*t_{i}(2)$, and $zx_{j} = x_{j}(1) + i*y_{j}(2)$. 
+  $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
   in this case are derivatives with respect to $z$, and in a slight abuse of
-  notation, we set $d/dz log(|z|) = 1/z$
+  notation, we set $\frac{\mathrm{d}}{\mathrm{d} z} \log(\|z\|) = 1/z$
 
 
 Laplace FMM (rfmm2d)
@@ -31,14 +31,14 @@ 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 the its gradient $\nabla u(x) \in
+the potential $u(x) \in \mathbb{R}$ and its gradient $\nabla u(x) \in
 \mathbb{R}^{2}$
 given by
 
 .. math::
    :label: rlap_nbody
 
-    u(x) = \sum_{j=1}^{N} c_{j} log(\|x-x_{j}\|) - v_{j} d_{j}\cdot \nabla log(\|x-x_{j}\|)  \, , 
+    u(x) = \sum_{j=1}^{N} c_{j} \log{(\|x-x_{j}\|)} - v_{j} d_{j}\cdot \nabla \log{(\|x-x_{j}\|)}  \, , 
 
 at the source and target locations. When $x=x_{j}$, the term
 corresponding to $x_{j}$ is dropped from the sum.
@@ -54,14 +54,14 @@ 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 the its gradient $\nabla u(x) \in
+the potential $u(x) \in \mathbb{C}$ and its gradient $\nabla u(x) \in
 \mathbb{C}^{2}$
 given by
 
 .. math::
    :label: llap_nbody
 
-    u(x) = \sum_{j=1}^{N} c_{j} log(\|x-x_{j}\|) - v_{j} d_{j}\cdot \nabla log(\|x-x_{j}\|)  \, , 
+    u(x) = \sum_{j=1}^{N} c_{j} \log{(\|x-x_{j}\|)} - v_{j} d_{j}\cdot \nabla \log{(\|x-x_{j}\|)}  \, , 
 
 at the source and target locations. When $x=x_{j}$, the term
 corresponding to $x_{j}$ is dropped from the sum.
@@ -76,18 +76,18 @@ $j=1,2,\ldots N$,
 denote a collection of dipole strengths.
 
 The Laplace FMM (cfmm2d) computes 
-the potential $u(z) \in \mathbb{C}$ and the its gradient $d/dz u(z) \in
+the potential $u(z) \in \mathbb{C}$ and its gradient $\frac{\mathrm{d}}{\mathrm{d}z} u(z) \in
 \mathbb{C}$
 given by
 
 .. math::
    :label: clap_nbody
 
-    u(z) = \sum_{j=1}^{N} c_{j} log(\|z-zx_{j}\|) - \frac{v_{j}}{z-zx_{j}}  \, , 
+    u(z) = \sum_{j=1}^{N} c_{j} \log{(\|z-\xi_{j}\|)} - \frac{v_{j}}{z-\xi_{j}}  \, , 
 
-at the source and target locations. When $z=zx_{j}$, the term
-corresponding to $zx_{j}$ is dropped from the sum. Here
-$z = x(1) + i x(2)$, and $zx_{j} = x_{j}(1) + ix_{j}(2)$, are the complex
+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.
 
 
@@ -103,7 +103,7 @@ Let $k\in\mathbb{C}$ denote the wave number or the Helmholtz
 parameter. 
 
 The Helmholtz FMM computes 
-the potential $u(x)$ and the its gradient $\nabla u(x)$
+the potential $u(x)$ and its gradient $\nabla u(x)$
 given by
 
 .. math::
@@ -118,13 +118,155 @@ corresponding to $x_{j}$ is dropped from the sum.
 
 Biharmonic FMM
 ***************
+Let $c_{j} = (c_{j,1}, c_{j,2})\in \mathbb{C}^2$, 
+$j=1,2,\ldots N$, 
+denote a collection of charge strengths, and 
+$v_{j} = (v_{j,1}, v_{j,2}, v_{j,3}) \in \mathbb{C}^{3}$,
+$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}})$
+given by
+
+.. math::
+   :label: biharm_nbody
+
+    u(z) &= \sum_{j=1}^{N} c_{j,1} \log{\|z - \xi_{j}\|} + 
+    c_{j,2} \frac{z-\xi_{j}}{\overline{z-\xi_{j}}} + \frac{v_{j,1}}{z - \xi_{j}} + 
+    \frac{v_{j,3}}{\overline{z-\xi_{j}}} + 
+    v_{j,2} \frac{z - \xi_{j}}{(\overline{z-\xi_{j}})^2} \, , \\
+    P_{z} \frac{\mathrm{d}}{\mathrm{d} z}u(z) &= \sum_{j=1}^{N} \frac{c_{j,1}}{z-\xi_{j}} -
+    \frac{v_{j,1}}{(z-\xi_{j})^2} \, \\
+    P_{\overline{z}} \frac{\mathrm{d}}{\mathrm{d} z} u(z) &= 
+    \sum_{j=1}^{N} \frac{c_{j,2}}{\overline{z-\xi_{j}}} + 
+    \frac{v_{j,2}}{(\overline{z-\xi_{j}})^2} \,
+    ,\\
+    \frac{\mathrm{d}}{\mathrm{d}\overline{z}} u(z) &= 
+    \sum_{j=1}^{N} \frac{c_{j,1}}{\overline{z-\xi_{j}}} - 
+    c_{j,2} \frac{z-\xi_{j}}{(\overline{z-\xi_{j}})^2} - 
+    \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
+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.
+
 
 Modified Biharmonic FMM
 ************************
 
+Let $G^{\textrm{mbh}}(x,y)$ denote the modified biharmonic
+Green's function given by
+
+.. math::
+    G^{\textrm{mbh}}(x,y) = \frac{1}{2\pi \beta^2}\left(K_{0}(\beta \|x-y \|) - \log{\|x-y\|}\right)
+
+where $K_{0}$ is the modified Bessel function of order $0$, and $\beta$ is the
+Modified Biharmonic wavenumber.
+
+
+Let $c_{j} \in \mathbb{C}$, 
+denote a collection of charge strengths, 
+$v_{j} \in \mathbb{C}$, 
+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 
+quadrupole strengths, 
+$w_{j} = (w_{j,1}, w_{j,2}, w_{j,3}) \in \mathbb{C}^{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,
+$j=1,2,\ldots N$.
+
+The Modified Biharmonic FMM computes 
+the potential $u(x)$ and its gradients, and hessians, 
+given by
+
+.. math::
+   :label: modbiharm_nbody
+
+    u(x) = \sum_{j=1}^{N} &c_{j} G^{\textrm{mbh}}(x,x_{j}) + 
+    v_{j} d_{j} \cdot \nabla_{y} G^{\textrm{mbh}}(x,x_{j}) + \\
+    &q_{j} \left(w_{j,1} \partial_{y_{1},y_{1}} G^{\textrm{mbh}}(x,x_{j}) + w_{j,2}
+    \partial_{y_{1},y_{2}} G^{\textrm{mbh}}(x,x_{j}) + w_{j,3}
+    \partial_{y_{2},y_{2}} G^{\textrm{mbh}}(x,x_{j}) \right) + \\
+    &o_{j} \big( p_{j,1} \partial_{y_{1},y_{1},y_{1}} G^{\textrm{mbh}}(x,x_{j}) +
+    p_{j,2} \partial_{y_{1},y_{1},y_{2}} G^{\textrm{mbh}}(x,x_{j}) + \\
+    &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) \, ,
+
+at the source and target locations. When $x=x_{j}$, the term
+corresponding to $x_{j}$ is dropped from the sum.
+
 Stokes FMM
 ************
 
+Let $G^{\textrm{stok}}(x,y)$ denote the Stokeslet given by
+
+.. math::
+    G^{\textrm{stok}}(x,y) = \frac{1}{2}
+    \begin{bmatrix}
+    -\log{\|x-y \|} +  \frac{(x_{1}-y_{1})^2}{\|x-y\|^2} & \frac{(x_{1}-y_{1})
+    (x_{2}-y_{2})}{\|x-y \|^2} \\
+    \frac{(x_{1}-y_{1})(x_{2}-y_{2})}{\|x-y \|^2} &
+    -\log{\|x-y \|} +  \frac{(x_{2}-y_{2})^2}{\|x-y \|^2} 
+    \end{bmatrix} \, ,
+
+$P^{\textrm{stok}}(x,y)$ denote the associated pressure tensor
+
+.. math::
+    P(x,y) = \frac{1}{\|x-y \|^2}\begin{bmatrix}
+    (x_{1}-y_{1}) &
+    (x_{2}-y_{2})
+    \end{bmatrix} \, ,
+
+$T^{\textrm{stok}}(x,y)$ denote the Stresslet whose action of a vector $v$
+is given by
+
+.. math::
+    v\cdot T^{\textrm{stok}}(x,y) = -\frac{2 v \cdot (x-y)}{\|x-y \|^4}
+    \begin{bmatrix}
+    (x_{1} - y_{1})^2 & (x_{1}-y_{1})(x_{2}-y_{2}) \\
+    (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
+
+.. math::
+    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})
+    \end{bmatrix}
+
+Let $c_{j} \in \mathbb{R}^2$, 
+$j=1,2,\ldots N$, 
+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.
+
+The Stokes FMM computes 
+the potential $u(x) \in \mathbb{R}^2$, its gradient $\nabla u(x) \in
+\mathbb{R}^{2\times 2}$, and the pressure $p$ given by
+
+.. math::
+   :label: stok_nbody
+
+    u(x) &= \sum_{j=1}^{N} G^{\textrm{stok}}(x,x_{j}) c_{j} + d_{j} \cdot
+    T^{\textrm{stok}}(x,x_{j}) \cdot v_{j} \, , \\
+    p(x) &= \sum_{j=1}^{N} c_{j} P^{\textrm{stok}}(x,x_{j}) + d_{j} \cdot
+    \Pi^{\textrm{stok}}(x,x_{j}) \cdot v_{j}^{T}
+
+at the source and target locations. When $x=x_{j}$, the term
+corresponding to $x_{j}$ is dropped from the sum.
+
 Vectorized versions   
 *******************
 The vectorized versions of the Helmholtz FMM, 

python/bhfmmexample.py

@@ -0,0 +1,34 @@
+#!/usr/bin/env python
+
+import fmm2dpy as fmm
+import numpy as np
+
+
+#
+#  This is a sample code to demonstrate how to use
+#  the fmm libraries
+#
+
+# sample with one density, sources to sources,
+# charge interactions, and potential only
+#
+n = 100000
+nd = 1
+sources = np.random.uniform(0,1,(2,n))
+eps = 10**(-5)
+
+charges = np.random.uniform(0,1,(2,n)) + 1j*np.random.uniform(0,1,(2,n))
+dipoles = np.random.uniform(0,1,(3,n)) + 1j*np.random.uniform(0,1,(3,n))
+out = fmm.bhfmm2d(eps=eps,sources=sources,charges=charges,dipoles=dipoles,pg=1)
+
+
+# sample with a vector of densities, sources to 
+# sources and targets, dipole interactions, 
+# potential and gradietns
+
+nd = 3
+nt = 1870
+targ = np.random.uniform(0,1,(2,nt))
+dipoles = np.random.uniform(0,1,(nd,3,n)) +  1j*np.random.uniform(0,1,(nd,3,n))
+out2 = fmm.bhfmm2d(eps=eps,sources=sources,dipoles=dipoles,\
+    targets=targ,nd=nd,pg=2,pgt=2)

src/biharmonic/bhfmm2d.f

@@ -51,7 +51,7 @@ subroutine bhfmm2d(nd,eps,ns,sources,ifcharge,charge,
 c       and velocity are interchangably used
 c     
 c
-c    vel(z) = \sum charge1_j *\log|z-z_j| + 
+c    vel(z) = \sum 2*charge1_j *\log|z-z_j| + 
 c                (charge2_j) (z-z_j)/(z-z_j)_bar+
 c                dip1_j/(z-z_j) + dip2_j/(z-z_j)_bar-
 c                dip3_j (z-z_j)/(z-z_j)^2_bar
@@ -65,8 +65,8 @@ subroutine bhfmm2d(nd,eps,ns,sources,ifcharge,charge,
 c    \phi(z) = \sum charge_j log(z-z_j)+
 c                     dip1/(z-z_j)
 c
-c    \psi(z) = \sum dippar2_bar/(z-z_j)+
-c               z_j_bar dippar1/(z-z_j)^2 +
+c    \psi(z) = \sum dip2_bar/(z-z_j)+
+c               z_j_bar dip1/(z-z_j)^2 +
 c
 c    In all of the sums above, the terms for which |z-z_{j}| <= 2^-51|b|
 c    where |b| is the size of the root box will be ignored.

src/modified-biharmonic/mbhfmm2d.f

@@ -36,10 +36,10 @@ subroutine mbhfmm2d(nd,eps,beta,ns,sources,ifcharge,charge,
 c     u(x) = sum_j charge(j)*G(x,y(j))
 c     + dipstr(j)*(G_{y1}(x,y(j))*dipvec(1,j) + G_{y2}*dipvec(2,j))
 c     + quadstr(j)*(G_{y1y1}(x,y(j))*quadvec(1,j)
-c             + G_{y1y2}*quadvec(2,j) + G_{y2y2}*quadvec(2,j))
+c             + G_{y1y2}*quadvec(2,j) + G_{y2y2}*quadvec(3,j))
 c     + octstr(j)*(G_{y1y1y1}(x,y(j))*octvec(1,j)
-c             + G_{y1y1y2}*octvec(2,j) + G_{y1y2y2}*octvec(2,j)
-c             + G_{y1y2y2}*octvec(2,j))
+c             + G_{y1y1y2}*octvec(2,j) + G_{y1y2y2}*octvec(3,j)
+c             + G_{y2y2y2}*octvec(4,j))
 c
 c     INPUT PARAMETERS:
 c      

src/stokes/stfmm2d.f

@@ -11,16 +11,16 @@
 c     The Stokeslet, G_{ij}, and its associated pressure tensor, P_j,
 c     (without the 1/2pi scaling) are
 c
-c   UPDATE formulae
-c
-c     G_{ij}(x,y) = (r_i r_j)/(2r^3) + delta_{ij}/(2r)
-c     P_j(x,y) = r_j/r^3
+c     G_{ij}(x,y) = (r_i r_j)/(2r^2) - delta_{ij}log(r)/(2)
+c     P_j(x,y) = r_j/r^2
 c
 c     The (Type I) stresslet, T_{ijk}, and its associated pressure
-c     tensor, PI_{jk}, (without the 1/4pi scaling) are
+c     tensor, PI_{jk}, (without the 1/2pi scaling) are
 c     
-c     T_{ijk}(x,y) = -3 r_i r_j r_k/ r^5
-c     PI_{jk} = -2 delta_{jk} + 6 r_j r_k/r^5      
+c     T_{ijk}(x,y) = -2 r_i r_j r_k/ r^4
+c     PI_{jk} = - delta_{jk}/r^2 + 2 r_j r_k/r^4      
+c
+c
 c