1 files changed, 48 insertions, 1 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;