/* integration/qc25c.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 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, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

struct fn_cauchy_params
{
  gsl_function *function;
  double singularity;
};

static double fn_cauchy (double t, void *params);

static void compute_moments (double cc, double *moment);

static void
qc25c (gsl_function * f, double a, double b, double c, 
       double *result, double *abserr, int *err_reliable);

static void
qc25c (gsl_function * f, double a, double b, double c, 
       double *result, double *abserr, int *err_reliable)
{
  double cc = (2 * c - b - a) / (b - a);

  if (fabs (cc) > 1.1)
    {
      double resabs, resasc;

      gsl_function weighted_function;
      struct fn_cauchy_params fn_params;

      fn_params.function = f;
      fn_params.singularity = c;

      weighted_function.function = &fn_cauchy;
      weighted_function.params = &fn_params;

      gsl_integration_qk15 (&weighted_function, a, b, result, abserr,
                            &resabs, &resasc);
      
      if (*abserr == resasc)
        {
          *err_reliable = 0;
        }
      else 
        {
          *err_reliable = 1;
        }

      return;
    }
  else
    {
      double cheb12[13], cheb24[25], moment[25];
      double res12 = 0, res24 = 0;
      size_t i;
      gsl_integration_qcheb (f, a, b, cheb12, cheb24);
      compute_moments (cc, moment);

      for (i = 0; i < 13; i++)
        {
          res12 += cheb12[i] * moment[i];
        }

      for (i = 0; i < 25; i++)
        {
          res24 += cheb24[i] * moment[i];
        }

      *result = res24;
      *abserr = fabs(res24 - res12) ;
      *err_reliable = 0;

      return;
    }
}

static double
fn_cauchy (double x, void *params)
{
  struct fn_cauchy_params *p = (struct fn_cauchy_params *) params;
  gsl_function *f = p->function;
  double c = p->singularity;
  return GSL_FN_EVAL (f, x) / (x - c);
}

static void
compute_moments (double cc, double *moment)
{
  size_t k;

  double a0 = log (fabs ((1.0 - cc) / (1.0 + cc)));
  double a1 = 2 + a0 * cc;

  moment[0] = a0;
  moment[1] = a1;

  for (k = 2; k < 25; k++)
    {
      double a2;

      if ((k % 2) == 0)
        {
          a2 = 2.0 * cc * a1 - a0;
        }
      else
        {
          const double km1 = k - 1.0;
          a2 = 2.0 * cc * a1 - a0 - 4.0 / (km1 * km1 - 1.0);
        }

      moment[k] = a2;

      a0 = a1;
      a1 = a2;
    }
}