/* integration/qawf.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 <math.h>
#include <float.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>

#include "initialise.c"
#include "append.c"
#include "qelg.c"

int
gsl_integration_qawf (gsl_function * f,
                      const double a,
                      const double epsabs,
                      const size_t limit,
                      gsl_integration_workspace * workspace,
                      gsl_integration_workspace * cycle_workspace,
                      gsl_integration_qawo_table * wf,
                      double *result, double *abserr)
{
  double area, errsum;
  double res_ext, err_ext;
  double correc, total_error = 0.0, truncation_error;

  size_t ktmin = 0;
  size_t iteration = 0;

  struct extrapolation_table table;

  double cycle;
  double omega = wf->omega;

  const double p = 0.9;
  double factor = 1;
  double initial_eps, eps;
  int error_type = 0;

  /* Initialize results */

  initialise (workspace, a, a);

  *result = 0;
  *abserr = 0;

  if (limit > workspace->limit)
    {
      GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
    }

  /* Test on accuracy */

  if (epsabs <= 0)
    {
      GSL_ERROR ("absolute tolerance epsabs must be positive", GSL_EBADTOL) ;
    }

  if (omega == 0.0)
    {
      if (wf->sine == GSL_INTEG_SINE)
        {
          /* The function sin(w x) f(x) is always zero for w = 0 */

          *result = 0;
          *abserr = 0;

          return GSL_SUCCESS;
        }
      else
        {
          /* The function cos(w x) f(x) is always f(x) for w = 0 */

          int status = gsl_integration_qagiu (f, a, epsabs, 0.0,
                                              cycle_workspace->limit,
                                              cycle_workspace,
                                              result, abserr);
          return status;
        }
    }

  if (epsabs > GSL_DBL_MIN / (1 - p))
    {
      eps = epsabs * (1 - p);
    }
  else
    {
      eps = epsabs;
    }

  initial_eps = eps;

  area = 0;
  errsum = 0;

  res_ext = 0;
  err_ext = GSL_DBL_MAX;
  correc = 0;

  cycle = (2 * floor (fabs (omega)) + 1) * M_PI / fabs (omega);

  gsl_integration_qawo_table_set_length (wf, cycle);

  initialise_table (&table);

  for (iteration = 0; iteration < limit; iteration++)
    {
      double area1, error1, reseps, erreps;

      double a1 = a + iteration * cycle;
      double b1 = a1 + cycle;

      double epsabs1 = eps * factor;

      int status = gsl_integration_qawo (f, a1, epsabs1, 0.0, limit,
                                         cycle_workspace, wf,
                                         &area1, &error1);

      append_interval (workspace, a1, b1, area1, error1);

      factor *= p;

      area = area + area1;
      errsum = errsum + error1;

      /* estimate the truncation error as 50 times the final term */

      truncation_error = 50 * fabs (area1);

      total_error = errsum + truncation_error;

      if (total_error < epsabs && iteration > 4)
        {
          goto compute_result;
        }

      if (error1 > correc)
        {
          correc = error1;
        }

      if (status)
        {
          eps = GSL_MAX_DBL (initial_eps, correc * (1.0 - p));
        }

      if (status && total_error < 10 * correc && iteration > 3)
        {
          goto compute_result;
        }

      append_table (&table, area);

      if (table.n < 2)
        {
          continue;
        }

      qelg (&table, &reseps, &erreps);

      ktmin++;

      if (ktmin >= 15 && err_ext < 0.001 * total_error)
        {
          error_type = 4;
        }

      if (erreps < err_ext)
        {
          ktmin = 0;
          err_ext = erreps;
          res_ext = reseps;

          if (err_ext + 10 * correc <= epsabs)
            break;
          if (err_ext <= epsabs && 10 * correc >= epsabs)
            break;
        }

    }

  if (iteration == limit)
    error_type = 1;

  if (err_ext == GSL_DBL_MAX)
    goto compute_result;

  err_ext = err_ext + 10 * correc;

  *result = res_ext;
  *abserr = err_ext;

  if (error_type == 0)
    {
      return GSL_SUCCESS ;
    }

  if (res_ext != 0.0 && area != 0.0)
    {
      if (err_ext / fabs (res_ext) > errsum / fabs (area))
        goto compute_result;
    }
  else if (err_ext > errsum)
    {
      goto compute_result;
    }
  else if (area == 0.0)
    {
      goto return_error;
    }

  if (error_type == 4)
    {
      err_ext = err_ext + truncation_error;
    }

  goto return_error;

compute_result:

  *result = area;
  *abserr = total_error;

return_error:

  if (error_type > 2)
    error_type--;

  if (error_type == 0)
    {
      return GSL_SUCCESS;
    }
  else if (error_type == 1)
    {
      GSL_ERROR ("number of iterations was insufficient", GSL_EMAXITER);
    }
  else if (error_type == 2)
    {
      GSL_ERROR ("cannot reach tolerance because of roundoff error",
                 GSL_EROUND);
    }
  else if (error_type == 3)
    {
      GSL_ERROR ("bad integrand behavior found in the integration interval",
                 GSL_ESING);
    }
  else if (error_type == 4)
    {
      GSL_ERROR ("roundoff error detected in the extrapolation table",
                 GSL_EROUND);
    }
  else if (error_type == 5)
    {
      GSL_ERROR ("integral is divergent, or slowly convergent",
                 GSL_EDIVERGE);
    }
  else
    {
      GSL_ERROR ("could not integrate function", GSL_EFAILED);
    }

}