From b2ba263aac7332ec1c19e5258e00d93655fd76b4 Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Mon, 24 Jun 2024 02:20:21 +0000 Subject: [PATCH] Multipole: move func Integrate into integrate.cxx --- Multipole/src/integrate.cxx | 56 +++++++++++++++++++++++++++++++++++++ Multipole/src/integrate.hxx | 6 ++++ Multipole/src/multipole.cxx | 1 + Multipole/src/utils.cxx | 56 ------------------------------------- Multipole/src/utils.hxx | 6 ---- 5 files changed, 63 insertions(+), 62 deletions(-) diff --git a/Multipole/src/integrate.cxx b/Multipole/src/integrate.cxx index 1084624e..fc790312 100644 --- a/Multipole/src/integrate.cxx +++ b/Multipole/src/integrate.cxx @@ -1,4 +1,5 @@ #include "integrate.hxx" +#include "utils.hxx" #include #include @@ -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 index ece11d09..6e326d9e 100644 --- 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 index 75afb3c7..82415b5e 100644 --- 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 index bf3bcd2f..0c90b37f 100644 --- a/Multipole/src/utils.cxx +++ b/Multipole/src/utils.cxx @@ -1,4 +1,3 @@ -#include "integrate.hxx" #include "utils.hxx" #include @@ -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 index 6fa6b98c..6e242532 100644 --- 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