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

[ian.hinder at aei.mpg.de]

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

Modified:
 /trunk/
  param.ccl
 /trunk/src/
  utils.cc

Log:
 Make trapezoidal integration method available via parameter

File Changes:

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

File [modified]: utils.cc
Delta lines: +41 -45
===================================================================
--- trunk/src/utils.cc	2013-10-09 17:56:25 UTC (rev 97)
+++ trunk/src/utils.cc	2013-10-09 17:56:35 UTC (rev 98)
@@ -345,58 +345,54 @@
   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"))
+  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)
     {
-      *outre = Midpoint2DIntegral(fr, ntheta, nphi, dth, dph);
-      *outim = Midpoint2DIntegral(fi, ntheta, nphi, dth, dph);
+      CCTK_WARN (CCTK_WARN_ABORT, "The Simpson integration method requires even ntheta and even nphi");
     }
-    else if (CCTK_Equals(integration_method, "Simpson"))
+    *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)
     {
-      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);
+      CCTK_WARN (CCTK_WARN_ABORT, "The Driscoll&Healy integration method requires even ntheta");
     }
-    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");
-    }
+    *outre = DriscollHealy2DIntegral(fr, ntheta, nphi, dth, dph);
+    *outim = DriscollHealy2DIntegral(fi, ntheta, nphi, dth, dph);
   }
   else
   {

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

File [modified]: param.ccl
Delta lines: +1 -0
===================================================================
--- trunk/param.ccl	2013-10-09 17:56:25 UTC (rev 97)
+++ trunk/param.ccl	2013-10-09 17:56:35 UTC (rev 98)
@@ -31,6 +31,7 @@
 {
   "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"