Make trapezoidal integration method available via parameter

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

2 files changed, 42 insertions, 45 deletions

diff --git a/param.ccl b/param.ccl
index b8ed861..85ec3d6 100644
--- a/param.ccl
+++ b/param.ccl
@@ -31,6 +31,7 @@ KEYWORD integration_method "How to do surface integrals" STEERABLE=always
 {
   "midpoint"      :: "Midpoint rule (1st order)"
   "Simpson"       :: "Simpson's rule (4th order) [requires even ntheta and nphi]"
+  "trapezoidal"   :: "Trapezoidal rule (2nd order)"
   "DriscollHealy" :: "Driscoll & Healy (exponentially convergent) [requires even ntheta]"
 } "midpoint"
 
diff --git a/src/utils.cc b/src/utils.cc
index dc2f879..58a08db 100644
--- a/src/utils.cc
+++ b/src/utils.cc
@@ -345,58 +345,54 @@ void Multipole_Integrate(int array_size, int nthetap,
   iu = Multipole_Index(0,1,ntheta);
   CCTK_REAL dph = ph[iu] - ph[il];
 
-  if (CCTK_Equals(integration_method, "Simpson") ||
-           CCTK_Equals(integration_method, "DriscollHealy") ||
-           CCTK_Equals(integration_method, "midpoint"))
-  {
-    static CCTK_REAL *fr = 0;
-    static CCTK_REAL *fi = 0;
-    static bool allocated_memory = false;
+  static CCTK_REAL *fr = 0;
+  static CCTK_REAL *fi = 0;
+  static bool allocated_memory = false;
 
-    // Construct an array for the real integrand
-    if (!allocated_memory)
-    {
-      fr = new CCTK_REAL[array_size];
-      fi = new CCTK_REAL[array_size];
-      allocated_memory = true;
-    }
+  // Construct an array for the real integrand
+  if (!allocated_memory)
+  {
+    fr = new CCTK_REAL[array_size];
+    fi = new CCTK_REAL[array_size];
+    allocated_memory = true;
+  }
     
-    // the below calculations take the integral of conj(array1)*array2*sin(th)
-    for (int i = 0; i < array_size; i++)
-    {
-      fr[i] = (array1r[i] * array2r[i] + 
-               array1i[i] * array2i[i] ) * sin(th[i]);
-      fi[i] = (array1r[i] * array2i[i] - 
-               array1i[i] * array2r[i] ) * sin(th[i]);
-    }
+  // the below calculations take the integral of conj(array1)*array2*sin(th)
+  for (int i = 0; i < array_size; i++)
+  {
+    fr[i] = (array1r[i] * array2r[i] + 
+             array1i[i] * array2i[i] ) * sin(th[i]);
+    fi[i] = (array1r[i] * array2i[i] - 
+             array1i[i] * array2r[i] ) * sin(th[i]);
+  }
     
-    if (CCTK_Equals(integration_method, "midpoint"))
-    {
-      *outre = Midpoint2DIntegral(fr, ntheta, nphi, dth, dph);
-      *outim = Midpoint2DIntegral(fi, ntheta, nphi, dth, dph);
-    }
-    else 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 (CCTK_Equals(integration_method, "midpoint"))
+  {
+    *outre = Midpoint2DIntegral(fr, ntheta, nphi, dth, dph);
+    *outim = Midpoint2DIntegral(fi, ntheta, nphi, dth, dph);
+  }
+  else if (CCTK_Equals(integration_method, "trapezoidal"))
+  {
+    *outre = Trapezoidal2DIntegral(fr, ntheta, nphi, dth, dph);
+    *outim = Trapezoidal2DIntegral(fi, ntheta, nphi, dth, dph);
+  }
+  else if (CCTK_Equals(integration_method, "Simpson"))
+  {
+    if (nphi % 2 != 0 || ntheta % 2 != 0)
     {
-      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);
+      CCTK_WARN (CCTK_WARN_ABORT, "The Simpson integration method requires even ntheta and even nphi");
     }
-    else
+    *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, "internal error");
+      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
   {