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"