Implement Driscoll&Healy integration

Implements a more accurate integration over the sphere using an
algorithm by Driscoll & Healy. This algorithm uses Gaussian
integration weights, leading (almost) to exponential convergence.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Multipole/trunk@83 4f5cb9a8-4dd8-4c2d-9bbd-173fa4467843

3 files changed, 79 insertions, 4 deletions

diff --git a/src/integrate.cc b/src/integrate.cc
index 0f1c392..1a9f44b 100644
--- a/src/integrate.cc
+++ b/src/integrate.cc
@@ -107,6 +107,52 @@ CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, CC
   return ((double) 1 / (double) 9) * hx * hy * integrand_sum;
 }
 
+// See: J.R. Driscoll and D.M. Healy Jr., Computing Fourier transforms
+// and convolutions on the 2-sphere. Advances in Applied Mathematics,
+// 15(2):202–250, 1994.
+CCTK_REAL DriscollHealy2DIntegral(CCTK_REAL const *const f,
+                                  int const nx, int const ny,
+                                  CCTK_REAL const hx, CCTK_REAL const hy)
+{
+  assert(f);
+  assert(nx >= 0);
+  assert(ny >= 0);
+  assert(nx % 2 == 0);
+  
+  CCTK_REAL integrand_sum = 0.0;
+  
+  // Skip the poles (ix=0 and ix=nx), since the weight there is zero
+  // anyway
+#pragma omp parallel for reduction(+: integrand_sum)
+  for (int ix = 1; ix < nx; ++ ix)
+  {
+    
+    // These weights lead to an almost spectral convergence
+    CCTK_REAL const theta = M_PI * ix / nx;
+    CCTK_REAL weight = 0.0;
+    for (int l = 0; l < nx/2; ++ l)
+    {
+      weight += sin((2*l+1)*theta)/(2*l+1);
+    }
+    weight *= 4.0 / M_PI;
+    
+    CCTK_REAL local_sum = 0.0;
+    // Skip the last point (iy=ny), since we assume periodicity and
+    // therefor it has the same value as the first point. We don't use
+    // weights in this direction, which leads to spectral convergence.
+    // (Yay periodicity!)
+    for (int iy = 0; iy < ny; ++ iy)
+    {
+      local_sum += f[idx(ix,iy)];
+    }
+    
+    integrand_sum += weight * local_sum;
+    
+  }
+  
+  return hx * hy * integrand_sum;
+}
+
 // 1D integrals
 
 static CCTK_REAL Simpson1DIntegral(CCTK_REAL const *f, int n, CCTK_REAL h)
@@ -114,7 +160,8 @@ static CCTK_REAL Simpson1DIntegral(CCTK_REAL const *f, int n, CCTK_REAL h)
   CCTK_REAL integrand_sum = 0;
   int i = 0;
 
-  assert(n > 0); assert(f);
+  assert(f);
+  assert(n > 0);
   assert(n % 2 == 0);
 
   int p = n / 2;
diff --git a/src/integrate.hh b/src/integrate.hh
index c45917e..5bcd43f 100644
--- a/src/integrate.hh
+++ b/src/integrate.hh
@@ -6,4 +6,7 @@
 CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, 
                             CCTK_REAL hx, CCTK_REAL hy);
 
+CCTK_REAL DriscollHealy2DIntegral(CCTK_REAL const *f, int nx, int ny,
+                                  CCTK_REAL hx, CCTK_REAL hy);
+
 #endif
diff --git a/src/utils.cc b/src/utils.cc
index 74d9eef..cb25bdf 100644
--- a/src/utils.cc
+++ b/src/utils.cc
@@ -346,7 +346,8 @@ void Multipole_Integrate(int array_size, int nthetap,
       *outim = tempi;
     } 
   }
-  else if (CCTK_Equals(integration_method, "simpson"))
+  else if (CCTK_Equals(integration_method, "Simpson") ||
+           CCTK_Equals(integration_method, "DriscollHealy"))
   {
     static CCTK_REAL *fr = 0;
     static CCTK_REAL *fi = 0;
@@ -368,7 +369,31 @@ void Multipole_Integrate(int array_size, int nthetap,
                array1i[i] * array2r[i] ) * sin(th[i]);
     }
     
-    *outre = Simpson2DIntegral(fr, ntheta, nphi, dth, dph);
-    *outim = Simpson2DIntegral(fi, ntheta, nphi, dth, dph);
+    if (CCTK_Equals(integration_method, "Simpson"))
+    {
+      if (nphi % 2 != 0 || ntheta % 2 != 0)
+      {
+        CCTK_WARN (CCTK_WARN_ABORT, "The Simpson integration method requires even ntheta and even nphi");
+      }
+      *outre = Simpson2DIntegral(fr, ntheta, nphi, dth, dph);
+      *outim = Simpson2DIntegral(fi, ntheta, nphi, dth, dph);
+    }
+    else if (CCTK_Equals(integration_method, "DriscollHealy"))
+    {
+      if (ntheta % 2 != 0)
+      {
+        CCTK_WARN (CCTK_WARN_ABORT, "The Driscoll&Healy integration method requires even ntheta");
+      }
+      *outre = DriscollHealy2DIntegral(fr, ntheta, nphi, dth, dph);
+      *outim = DriscollHealy2DIntegral(fi, ntheta, nphi, dth, dph);
+    }
+    else
+    {
+      CCTK_WARN(CCTK_WARN_ABORT, "internal error");
+    }
+  }
+  else
+  {
+    CCTK_WARN(CCTK_WARN_ABORT, "internal error");
   }
 }