1 files changed, 265 insertions, 0 deletions
diff --git a/gsl-1.9/integration/qag.c b/gsl-1.9/integration/qag.c
new file mode 100644
index 0000000..f5945a8
--- /dev/null
+++ b/gsl-1.9/integration/qag.c
@@ -0,0 +1,265 @@
+/* integration/qag.c
+ * 
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
+ * 
+ * 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 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <config.h>
+#include <stdlib.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_integration.h>
+
+#include "initialise.c"
+#include "set_initial.c"
+#include "qpsrt.c"
+#include "util.c"
+
+static int
+qag (const gsl_function *f,
+     const double a, const double b,
+     const double epsabs, const double epsrel,
+     const size_t limit,
+     gsl_integration_workspace * workspace,
+     double * result, double * abserr,
+     gsl_integration_rule * q) ;
+
+int
+gsl_integration_qag (const gsl_function *f,
+                     double a, double b,
+                     double epsabs, double epsrel, size_t limit,
+                     int key,
+                     gsl_integration_workspace * workspace,
+                     double * result, double * abserr)
+{
+  int status ;
+  gsl_integration_rule * integration_rule = gsl_integration_qk15 ;
+
+  if (key < GSL_INTEG_GAUSS15)
+    {
+      key = GSL_INTEG_GAUSS15 ;
+    } 
+  else if (key > GSL_INTEG_GAUSS61) 
+    {
+      key = GSL_INTEG_GAUSS61 ;
+    }
+
+  switch (key) 
+    {
+    case GSL_INTEG_GAUSS15:
+      integration_rule = gsl_integration_qk15 ;
+      break ;
+    case GSL_INTEG_GAUSS21:
+      integration_rule = gsl_integration_qk21 ;
+      break ;
+    case GSL_INTEG_GAUSS31:
+      integration_rule = gsl_integration_qk31 ; 
+      break ;
+    case GSL_INTEG_GAUSS41:
+      integration_rule = gsl_integration_qk41 ;
+      break ;      
+    case GSL_INTEG_GAUSS51:
+      integration_rule = gsl_integration_qk51 ;
+      break ;      
+    case GSL_INTEG_GAUSS61:
+      integration_rule = gsl_integration_qk61 ;
+      break ;      
+    default:
+      GSL_ERROR("value of key does specify a known integration rule", 
+                GSL_EINVAL) ;
+    }
+
+  status = qag (f, a, b, epsabs, epsrel, limit,
+                workspace, 
+                result, abserr, 
+                integration_rule) ;
+  
+  return status ;
+}
+
+static int
+qag (const gsl_function * f,
+     const double a, const double b,
+     const double epsabs, const double epsrel,
+     const size_t limit,
+     gsl_integration_workspace * workspace,
+     double *result, double *abserr,
+     gsl_integration_rule * q)
+{
+  double area, errsum;
+  double result0, abserr0, resabs0, resasc0;
+  double tolerance;
+  size_t iteration = 0;
+  int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
+
+  double round_off;     
+
+  /* Initialize results */
+
+  initialise (workspace, a, b);
+
+  *result = 0;
+  *abserr = 0;
+
+  if (limit > workspace->limit)
+    {
+      GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
+    }
+
+  if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
+    {
+      GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
+                 GSL_EBADTOL);
+    }
+
+  /* perform the first integration */
+
+  q (f, a, b, &result0, &abserr0, &resabs0, &resasc0);
+
+  set_initial_result (workspace, result0, abserr0);
+
+  /* Test on accuracy */
+
+  tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
+
+  /* need IEEE rounding here to match original quadpack behavior */
+
+  round_off = GSL_COERCE_DBL (50 * GSL_DBL_EPSILON * resabs0);
+
+  if (abserr0 <= round_off && abserr0 > tolerance)
+    {
+      *result = result0;
+      *abserr = abserr0;
+
+      GSL_ERROR ("cannot reach tolerance because of roundoff error "
+                 "on first attempt", GSL_EROUND);
+    }
+  else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)
+    {
+      *result = result0;
+      *abserr = abserr0;
+
+      return GSL_SUCCESS;
+    }
+  else if (limit == 1)
+    {
+      *result = result0;
+      *abserr = abserr0;
+
+      GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
+    }
+
+  area = result0;
+  errsum = abserr0;
+
+  iteration = 1;
+
+  do
+    {
+      double a1, b1, a2, b2;
+      double a_i, b_i, r_i, e_i;
+      double area1 = 0, area2 = 0, area12 = 0;
+      double error1 = 0, error2 = 0, error12 = 0;
+      double resasc1, resasc2;
+      double resabs1, resabs2;
+
+      /* Bisect the subinterval with the largest error estimate */
+
+      retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
+
+      a1 = a_i; 
+      b1 = 0.5 * (a_i + b_i);
+      a2 = b1;
+      b2 = b_i;
+
+      q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);
+      q (f, a2, b2, &area2, &error2, &resabs2, &resasc2);
+
+      area12 = area1 + area2;
+      error12 = error1 + error2;
+
+      errsum += (error12 - e_i);
+      area += area12 - r_i;
+
+      if (resasc1 != error1 && resasc2 != error2)
+        {
+          double delta = r_i - area12;
+
+          if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
+            {
+              roundoff_type1++;
+            }
+          if (iteration >= 10 && error12 > e_i)
+            {
+              roundoff_type2++;
+            }
+        }
+
+      tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
+
+      if (errsum > tolerance)
+        {
+          if (roundoff_type1 >= 6 || roundoff_type2 >= 20)
+            {
+              error_type = 2;   /* round off error */
+            }
+
+          /* set error flag in the case of bad integrand behaviour at
+             a point of the integration range */
+
+          if (subinterval_too_small (a1, a2, b2))
+            {
+              error_type = 3;
+            }
+        }
+
+      update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
+
+      retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
+
+      iteration++;
+
+    }
+  while (iteration < limit && !error_type && errsum > tolerance);
+
+  *result = sum_results (workspace);
+  *abserr = errsum;
+
+  if (errsum <= tolerance)
+    {
+      return GSL_SUCCESS;
+    }
+  else if (error_type == 2)
+    {
+      GSL_ERROR ("roundoff error prevents tolerance from being achieved",
+                 GSL_EROUND);
+    }
+  else if (error_type == 3)
+    {
+      GSL_ERROR ("bad integrand behavior found in the integration interval",
+                 GSL_ESING);
+    }
+  else if (iteration == limit)
+    {
+      GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER);
+    }
+  else
+    {
+      GSL_ERROR ("could not integrate function", GSL_EFAILED);
+    }
+}
+
+