diff options
Diffstat (limited to 'gsl-1.9/integration/qmomo.c')
| -rw-r--r-- | gsl-1.9/integration/qmomo.c | 193 |
1 files changed, 193 insertions, 0 deletions
diff --git a/gsl-1.9/integration/qmomo.c b/gsl-1.9/integration/qmomo.c new file mode 100644 index 0000000..62125a9 --- /dev/null +++ b/gsl-1.9/integration/qmomo.c @@ -0,0 +1,193 @@ +/* integration/qmomo.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_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) +{ + 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; + } +} |