From 167b7353013f12c35d935a88a91c7c961ae5a2db Mon Sep 17 00:00:00 2001
From: eschnett <eschnett@4f5cb9a8-4dd8-4c2d-9bbd-173fa4467843>
Date: Sat, 16 Jun 2012 15:39:07 +0000
Subject: 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
---
 src/utils.cc | 31 ++++++++++++++++++++++++++++---
 1 file changed, 28 insertions(+), 3 deletions(-)

(limited to 'src/utils.cc')

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");
   }
 } 
-- 
cgit v1.2.3