make trapzd() public

knodlseder
1 parent 784e9fe5
Showing 2 changed files with 153 additions and 139 deletions
include/GIntegral.hpp
src/numerics/GIntegral.cpp
@@ -4,10 +4,18 @@
  *  copyright (C) 2010-2011 by Jurgen Knodlseder                           *
  * ----------------------------------------------------------------------- *
  *                                                                         *
- *   This program is free software; you can redistribute it and/or modify  *
- *   it under the terms of the GNU General Public License as published by  *
- *   the Free Software Foundation; either version 2 of the License, or     *
- *   (at your option) any later version.                                   *
+ *  This program is free software: you can redistribute it and/or modify   *
+ *  it under the terms of the GNU General Public License as published by   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
  *                                                                         *
  ***************************************************************************/
 /**
@@ -64,6 +72,7 @@ public:
     void              integrand(GIntegrand* integrand) { m_integrand=integrand; }
     const GIntegrand* integrand(void) const { return m_integrand; }
     double            romb(double a, double b, int k = 5);
+    double            trapzd(double a, double b, int n = 1, double result = 0.0);
     std::string       print(void) const;
 protected:
@@ -72,7 +81,6 @@ protected:
     void   copy_members(const GIntegral& integral);
     void   free_members(void);
     double polint(double* xa, double* ya, int n, double x, double *dy);
-    double trapzd(double a, double b, int n, double result);
     // Protected data area
     GIntegrand* m_integrand;    //!< Pointer to integrand
@@ -196,63 +196,157 @@ double GIntegral::romb(double a, double b, int k)
     double result    = 0.0;
     double ss        = 0.0;
     double dss       = 0.0;
+    
+    // Continue only if integration range is valid
+    if (b > a) {
+
+        // Allocate temporal storage
+        double* s = new double[m_max_iter+2];
+        double* h = new double[m_max_iter+2];
+
+        // Initialise step size
+        h[1] = 1.0;
+        s[0] = 0.0;
+
+        // Iterative loop
+        for (m_iter = 1; m_iter <= m_max_iter; ++m_iter) {
+
+            // Integration using Trapezoid rule
+            s[m_iter] = trapzd(a, b, m_iter, s[m_iter-1]);
+
+            // Compile option: Check for NaN/Inf
+            #if defined(G_NAN_CHECK)
+            if (isnotanumber(s[m_iter]) || isinfinite(s[m_iter])) {
+                std::cout << "*** ERROR: GIntegral::romb";
+                std::cout << "(a=" << a << ", b=" << b << ", k=" << k << "):";
+                std::cout << " NaN/Inf encountered";
+                std::cout << " (s[" << m_iter << "]=" << s[m_iter] << ")";
+                std::cout << std::endl;
+            }
+            #endif
+
+            // Starting from iteration k on, use polynomial interpolation
+            if (m_iter >= k) {
+                ss = polint(&h[m_iter-k], &s[m_iter-k], k, 0.0, &dss);
+                if (std::abs(dss) <= m_eps * std::abs(ss)) {
+                    converged = true;
+                    result    = ss;
+                    break;
+                }
+            }
-    // Allocate temporal storage
-    double* s = new double[m_max_iter+2];
-    double* h = new double[m_max_iter+2];
-
-    // Initialise step size
-    h[1] = 1.0;
-    s[0] = 0.0;
+            // Reduce step size
+            h[m_iter+1]= 0.25 * h[m_iter];
-    // Loop
-    for (m_iter = 1; m_iter <= m_max_iter; ++m_iter) {
+        } // endfor: iterative loop
-        // Integration using Trapezoid rule
-        s[m_iter] = trapzd(a, b, m_iter, s[m_iter-1]);
+        // Free temporal storage
+        delete [] s;
+        delete [] h;
-        // Compile option: Check for NaN/Inf
-        #if defined(G_NAN_CHECK)
-        if (std::isnan(s[m_iter]) || std::isinf(s[m_iter])) {
-            std::cout << "*** ERROR: GIntegral::romb";
-            std::cout << "(a=" << a << ", b=" << b << ", k=" << k << "):";
-            std::cout << " NaN/Inf encountered";
-            std::cout << " (s[" << m_iter << "]=" << s[m_iter] << ")";
-            std::cout << std::endl;
-        }
-        #endif
-
-        // Starting from iteration k on, use polynomial interpolation
-        if (m_iter >= k) {
-            ss = polint(&h[m_iter-k], &s[m_iter-k], k, 0.0, &dss);
-            if (std::abs(dss) <= m_eps * std::abs(ss)) {
-                converged = true;
-                result    = ss;
-                break;
+        // Dump warning
+        if (!m_silent) {
+            if (!converged) {
+                std::cout << "*** WARNING: GIntegral::romb: ";
+                std::cout << "Integration did not converge ";
+                std::cout << "(iter=" << m_iter;
+                std::cout << ", result=" << ss;
+                std::cout << ", d=" << std::abs(dss);
+                std::cout << " > " << m_eps * std::abs(ss) << ")";
+                std::cout << std::endl;
+                throw;
             }
         }
+    
+    } // endif: integration range was valid
+
+    // Return result
+    return result;
+}
-        // Reduce step size
-        h[m_iter+1]= 0.25 * h[m_iter];
+/***********************************************************************//**
+ * @brief Perform Trapezoidal integration
+ *
+ * @param[in] a Left integration boundary.
+ * @param[in] b Right integration boundary.
+ * @param[in] n Number of steps.
+ * @param[in] result Result from a previous trapezoidal integration step.
+ *
+ * The original Numerical Recipes function had result declared as a static
+ * variable, yet this led to some untrackable integration problems. For this
+ * reason, previous results are now passed using an argument.
+ * Result initialisation is done if n=1.
+ ***************************************************************************/
+double GIntegral::trapzd(double a, double b, int n, double result)
+{
+    // Handle case of identical boundaries
+    if (a == b) {
+        result = 0.0;
     }
+    
+    // ... otherwise use trapeziodal rule
+    else {
+    
+        // Case A: Only a single step is requested
+        if (n == 1) {
+        
+            // Evaluate integrand at boundaries
+            double y_a = m_integrand->eval(a);
+            double y_b = m_integrand->eval(b);
+            
+            // Compute result
+            result = 0.5*(b-a)*(m_integrand->eval(a) + m_integrand->eval(b));
+            
+        } // endif: only a single step was requested
-    // Free temporal storage
-    delete [] s;
-    delete [] h;
-
-    // Dump warning
-    if (!m_silent) {
-        if (!converged) {
-            std::cout << "*** WARNING: GIntegral::romb: ";
-            std::cout << "Integration did not converge ";
-            std::cout << "(iter=" << m_iter;
-            std::cout << ", result=" << ss;
-            std::cout << ", d=" << std::abs(dss);
-            std::cout << " > " << m_eps * std::abs(ss) << ")";
-            std::cout << std::endl;
+        // Case B: More than a single step is requested
+        else {
+            // Compute step level 2^(n-1)
+            int it = 1;
+            for (int j = 1; j < n-1; ++j)
+                it <<= 1;
+
+            // Verify that step level is valid
+            if (it == 0) {
+                std::cout << "*** ERROR: GIntegral::trapzd(";
+                std::cout << "a=" << a << ", b=" << b << ", n=" << n;
+                std::cout << ", result=" << result << "): it=" << it << ". ";
+                std::cout << "Looks like an overflow?";
+                std::cout << std::endl;
+            }
+
+            // Set step size
+            double tnm = double(it);
+            double del = (b-a)/tnm;
+
+            // Verify that step is >0
+            if (del == 0) {
+                std::cout << "*** ERROR: GIntegral::trapzd(";
+                std::cout << "a=" << a << ", b=" << b << ", n=" << n;
+                std::cout << ", result=" << result << "): del=" << del << ". ";
+                std::cout << "Step is too small to make sense.";
+                std::cout << std::endl;
+            }
+
+            // Sum up values
+            double x   = a + 0.5*del;
+            double sum = 0.0;
+            for (int j = 0; j < it; ++j, x+=del) {
+                
+                // Evaluate integrand
+                double y = m_integrand->eval(x);
+
+                // Add integrand
+                sum += y;
+                
+            } // endfor: looped over steps
+
+            // Set result
+            result = 0.5*(result + (b-a)*sum/tnm);
         }
-    }
+        
+    } // endelse: trapeziodal rule was applied
     // Return result
     return result;
@@ -411,94 +505,6 @@ double GIntegral::polint(double* xa, double* ya, int n, double x, double *dy)
 }
-/***********************************************************************//**
- * @brief Perform Trapezoidal integration
- *
- * @param[in] a Left integration boundary.
- * @param[in] b Right integration boundary.
- * @param[in] n Number of steps.
- * @param[in] result Result from a previous trapezoidal integration step.
- *
- * The original Numerical Recipes function had result declared as a static
- * variable, yet this led to some untrackable integration problems. For this
- * reason, previous results are now passed using an argument.
- * Result initialisation is done if n=1.
- ***************************************************************************/
-double GIntegral::trapzd(double a, double b, int n, double result)
-{
-    // Handle case of identical boundaries
-    if (a == b) {
-        result = 0.0;
-    }
-    
-    // ... otherwise use trapeziodal rule
-    else {
-    
-        // Case A: Only a single step is requested
-        if (n == 1) {
-        
-            // Evaluate integrand at boundaries
-            double y_a = m_integrand->eval(a);
-            double y_b = m_integrand->eval(b);
-            
-            // Compute result
-            result = 0.5*(b-a)*(m_integrand->eval(a) + m_integrand->eval(b));
-            
-        } // endif: only a single step was requested
-
-        // Case B: More than a single step is requested
-        else {
-            // Compute step level 2^(n-1)
-            int it = 1;
-            for (int j = 1; j < n-1; ++j)
-                it <<= 1;
-
-            // Verify that step level is valid
-            if (it == 0) {
-                std::cout << "*** ERROR: GIntegral::trapzd(";
-                std::cout << "a=" << a << ", b=" << b << ", n=" << n;
-                std::cout << ", result=" << result << "): it=" << it << ". ";
-                std::cout << "Looks like an overflow?";
-                std::cout << std::endl;
-            }
-
-            // Set step size
-            double tnm = double(it);
-            double del = (b-a)/tnm;
-
-            // Verify that step is >0
-            if (it == 0) {
-                std::cout << "*** ERROR: GIntegral::trapzd(";
-                std::cout << "a=" << a << ", b=" << b << ", n=" << n;
-                std::cout << ", result=" << result << "): del=" << del << ". ";
-                std::cout << "Step is too small to make sense.";
-                std::cout << std::endl;
-            }
-
-            // Sum up values
-            double x   = a + 0.5*del;
-            double sum = 0.0;
-            for (int j = 0; j < it; ++j, x+=del) {
-                
-                // Evaluate integrand
-                double y = m_integrand->eval(x);
-
-                // Add integrand
-                sum += y;
-                
-            } // endfor: looped over steps
-
-            // Set result
-            result = 0.5*(result + (b-a)*sum/tnm);
-        }
-        
-    } // endelse: trapeziodal rule was applied
-
-    // Return result
-    return result;
-}
-
-
 /*==========================================================================
  =                                                                         =
  =                                 Friends                                 =