Multipole: move func Integrate into integrate.cxx

EinsteinToolkit · Jun 24, 2024 · b2ba263 · b2ba263
1 parent d625843
commit b2ba263
Show file tree

Hide file tree

Showing 5 changed files with 63 additions and 62 deletions.
diff --git a/Multipole/src/integrate.cxx b/Multipole/src/integrate.cxx
@@ -1,4 +1,5 @@
 #include "integrate.hxx"
+#include "utils.hxx"
 
 #include <assert.h>
 #include <math.h>
@@ -170,6 +171,61 @@ CCTK_REAL DriscollHealy2DIntegral(CCTK_REAL const *const f, int const nx,
   return hx * hy * integrand_sum;
 }
 
+void Integrate(int array_size, int nthetap, CCTK_REAL const array1r[],
+               CCTK_REAL const array1i[], CCTK_REAL const array2r[],
+               CCTK_REAL const array2i[], CCTK_REAL const th[],
+               CCTK_REAL const ph[], CCTK_REAL *outre, CCTK_REAL *outim) {
+  DECLARE_CCTK_PARAMETERS;
+
+  int il = Index_2d(0, 0, ntheta);
+  int iu = Index_2d(1, 0, ntheta);
+  CCTK_REAL dth = th[iu] - th[il];
+  iu = Index_2d(0, 1, ntheta);
+  CCTK_REAL dph = ph[iu] - ph[il];
+
+  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;
+  }
+
+  // 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, "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) {
+      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");
+  }
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 // 1D integrals
 ////////////////////////////////////////////////////////////////////////////////

diff --git a/Multipole/src/integrate.hxx b/Multipole/src/integrate.hxx
@@ -19,6 +19,12 @@ CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx,
 CCTK_REAL DriscollHealy2DIntegral(CCTK_REAL const *f, int nx, int ny,
                                   CCTK_REAL hx, CCTK_REAL hy);
 
+void Integrate(int array_size, int ntheta, CCTK_REAL const array1r[],
+               CCTK_REAL const array1i[], CCTK_REAL const array2r[],
+               CCTK_REAL const array2i[], CCTK_REAL const th[],
+               CCTK_REAL const pph[], CCTK_REAL out_arrayr[],
+               CCTK_REAL out_arrayi[]);
+
 } // namespace Multipole
 
 #endif // #ifndef MULTIPOLE_INTEGRATE_HXX
diff --git a/Multipole/src/multipole.cxx b/Multipole/src/multipole.cxx
@@ -1,3 +1,4 @@
+#include "integrate.hxx"
 #include "interpolate.hxx"
 #include "multipole.hxx"
 #include "sphericalharmonic.hxx"

diff --git a/Multipole/src/utils.cxx b/Multipole/src/utils.cxx
@@ -1,4 +1,3 @@
-#include "integrate.hxx"
 #include "utils.hxx"
 
 #include <errno.h>
@@ -291,59 +290,4 @@ void ScaleCartesian(int ntheta, int nphi, CCTK_REAL r, CCTK_REAL const xhat[],
   }
 }
 
-void Integrate(int array_size, int nthetap, CCTK_REAL const array1r[],
-               CCTK_REAL const array1i[], CCTK_REAL const array2r[],
-               CCTK_REAL const array2i[], CCTK_REAL const th[],
-               CCTK_REAL const ph[], CCTK_REAL *outre, CCTK_REAL *outim) {
-  DECLARE_CCTK_PARAMETERS;
-
-  int il = Index_2d(0, 0, ntheta);
-  int iu = Index_2d(1, 0, ntheta);
-  CCTK_REAL dth = th[iu] - th[il];
-  iu = Index_2d(0, 1, ntheta);
-  CCTK_REAL dph = ph[iu] - ph[il];
-
-  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;
-  }
-
-  // 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, "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) {
-      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");
-  }
-}
-
 } // namespace Multipole
diff --git a/Multipole/src/utils.hxx b/Multipole/src/utils.hxx
@@ -33,12 +33,6 @@ void ScaleCartesian(int ntheta, int nphi, CCTK_REAL r, CCTK_REAL const xhat[],
                     CCTK_REAL const yhat[], CCTK_REAL const zhat[],
                     CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[]);
 
-void Integrate(int array_size, int ntheta, CCTK_REAL const array1r[],
-               CCTK_REAL const array1i[], CCTK_REAL const array2r[],
-               CCTK_REAL const array2i[], CCTK_REAL const th[],
-               CCTK_REAL const pph[], CCTK_REAL out_arrayr[],
-               CCTK_REAL out_arrayi[]);
-
 } // namespace Multipole
 
 #endif // #ifndef MULTIPOLE_UTILS_HXX