[Commits] [svn:einsteintoolkit] Multipole/trunk/src/ (Rev. 93)

[ian.hinder at aei.mpg.de]

User: hinder
Date: 2013/10/09 12:55 PM

Modified:
 /trunk/src/
  utils.cc

Log:
 utils.cc: Treat all integration methods uniformly
 
 All tests pass

File Changes:

Directory: /trunk/src/
======================

File [modified]: utils.cc
Delta lines: +10 -19
===================================================================
--- trunk/src/utils.cc	2013-10-09 17:55:39 UTC (rev 92)
+++ trunk/src/utils.cc	2013-10-09 17:55:48 UTC (rev 93)
@@ -345,25 +345,10 @@
   iu = Multipole_Index(0,1,ntheta);
   CCTK_REAL dph = ph[iu] - ph[il];
 
-  if (CCTK_Equals(integration_method, "midpoint"))
+  if (CCTK_Equals(integration_method, "Simpson") ||
+           CCTK_Equals(integration_method, "DriscollHealy") ||
+           CCTK_Equals(integration_method, "midpoint"))
   {
-    CCTK_REAL tempr = 0.0;
-    CCTK_REAL tempi = 0.0;
-      
-    for (int i=0; i < array_size; i++) {
-      // the below calculations take the integral of conj(array1)*array2
-      tempr += ( array1r[i]*array2r[i] + array1i[i]*array2i[i] )
-        *sin(th[i])*dth*dph;
-      tempi += ( array1r[i]*array2i[i] - array1i[i]*array2r[i] )
-        *sin(th[i])*dth*dph;
-      
-      *outre = tempr;
-      *outim = tempi;
-    } 
-  }
-  else if (CCTK_Equals(integration_method, "Simpson") ||
-           CCTK_Equals(integration_method, "DriscollHealy"))
-  {
     static CCTK_REAL *fr = 0;
     static CCTK_REAL *fi = 0;
     static bool allocated_memory = false;
@@ -376,6 +361,7 @@
       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] + 
@@ -384,8 +370,13 @@
                array1i[i] * array2r[i] ) * sin(th[i]);
     }
     
-    if (CCTK_Equals(integration_method, "Simpson"))
+    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");