/* integration/qmomo.c

 * 

 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2009 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.

 */

#include <config.h>

#include <stdlib.h>

#include <gsl/gsl_integration.h>

#include <gsl/gsl_errno.h>

static void

initialise (double * ri, double * rj, double * rg, double * rh,

            double alpha, double beta);

gsl_integration_qaws_table * 

gsl_integration_qaws_table_alloc (double alpha, double beta, int mu, int nu)

{

  gsl_integration_qaws_table * t;

  if (alpha < -1.0)

    {

      GSL_ERROR_VAL ("alpha must be greater than -1.0", GSL_EINVAL, 0);

    }

  if (beta < -1.0)

    {

      GSL_ERROR_VAL ("beta must be greater than -1.0", GSL_EINVAL, 0);

    }

  if (mu != 0 && mu != 1)

    {

      GSL_ERROR_VAL ("mu must be 0 or 1", GSL_EINVAL, 0);

    }

  if (nu != 0 && nu != 1)

    {

      GSL_ERROR_VAL ("nu must be 0 or 1", GSL_EINVAL, 0);

    }

  t = (gsl_integration_qaws_table *) 

    malloc(sizeof(gsl_integration_qaws_table));

  if (t == 0)

    {

      GSL_ERROR_VAL ("failed to allocate space for qaws_table struct",

                        GSL_ENOMEM, 0);

    }

  t->alpha = alpha;

  t->beta = beta;

  t->mu = mu;

  t->nu = nu;

  initialise (t->ri, t->rj, t->rg, t->rh, alpha, beta);

  return t;

}

int

gsl_integration_qaws_table_set (gsl_integration_qaws_table * t,

                                double alpha, double beta, int mu, int nu)

{

  if (alpha < -1.0)

    {

      GSL_ERROR ("alpha must be greater than -1.0", GSL_EINVAL);

    }

  if (beta < -1.0)

    {

      GSL_ERROR ("beta must be greater than -1.0", GSL_EINVAL);

    }

  if (mu != 0 && mu != 1)

    {

      GSL_ERROR ("mu must be 0 or 1", GSL_EINVAL);

    }

  if (nu != 0 && nu != 1)

    {

      GSL_ERROR ("nu must be 0 or 1", GSL_EINVAL);

    }

  t->alpha = alpha;

  t->beta = beta;

  t->mu = mu;

  t->nu = nu;

  initialise (t->ri, t->rj, t->rg, t->rh, alpha, beta);

  return GSL_SUCCESS;

}

void

gsl_integration_qaws_table_free (gsl_integration_qaws_table * t)

{

  RETURN_IF_NULL (t);

  free (t);

}

static void

initialise (double * ri, double * rj, double * rg, double * rh,

            double alpha, double beta)

{

  const double alpha_p1 = alpha + 1.0;

  const double beta_p1 = beta + 1.0;

  const double alpha_p2 = alpha + 2.0;

  const double beta_p2 = beta + 2.0;

  const double r_alpha = pow (2.0, alpha_p1);

  const double r_beta = pow (2.0, beta_p1);

  size_t i;

  double an, anm1;

  ri[0] = r_alpha / alpha_p1;

  ri[1] = ri[0] * alpha / alpha_p2;

  an = 2.0;

  anm1 = 1.0;

  for (i = 2; i < 25; i++)

    {

      ri[i] = -(r_alpha + an * (an - alpha_p2) * ri[i - 1])

        / (anm1 * (an + alpha_p1));

      anm1 = an;

      an = an + 1.0;

    }

  rj[0] = r_beta / beta_p1;

  rj[1] = rj[0] * beta / beta_p2;

  an = 2.0;

  anm1 = 1.0;

  for (i = 2; i < 25; i++)

    {

      rj[i] = -(r_beta + an * (an - beta_p2) * rj[i - 1])

        / (anm1 * (an + beta_p1));

      anm1 = an;

      an = an + 1.0;

    }

  rg[0] = -ri[0] / alpha_p1;

  rg[1] = -rg[0] - 2.0 * r_alpha / (alpha_p2 * alpha_p2);

  an = 2.0;

  anm1 = 1.0;

  for (i = 2; i < 25; i++)

    {

      rg[i] = -(an * (an - alpha_p2) * rg[i - 1] - an * ri[i - 1]

                + anm1 * ri[i]) / (anm1 * (an + alpha_p1));

      anm1 = an;

      an = an + 1.0;

    }

  rh[0] = -rj[0] / beta_p1;

  rh[1] = -rh[0] - 2.0 * r_beta / (beta_p2 * beta_p2);

  an = 2.0;

  anm1 = 1.0;

  for (i = 2; i < 25; i++)

    {

      rh[i] = -(an * (an - beta_p2) * rh[i - 1] - an * rj[i - 1]

                + anm1 * rj[i]) / (anm1 * (an + beta_p1));

      anm1 = an;

      an = an + 1.0;

    }

  for (i = 1; i < 25; i += 2)

    {

      rj[i] *= -1;

      rh[i] *= -1;

    }

}