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

1 files changed, 28 insertions, 3 deletions

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");
   }
 }