diff options
Diffstat (limited to 'gsl-1.9/integration/qc25s.c')
| -rw-r--r-- | gsl-1.9/integration/qc25s.c | 238 |
1 files changed, 238 insertions, 0 deletions
diff --git a/gsl-1.9/integration/qc25s.c b/gsl-1.9/integration/qc25s.c new file mode 100644 index 0000000..196e30d --- /dev/null +++ b/gsl-1.9/integration/qc25s.c @@ -0,0 +1,238 @@ +/* integration/qc25s.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. + */ + +struct fn_qaws_params +{ + gsl_function *function; + double a; + double b; + gsl_integration_qaws_table *table; +}; + +static double fn_qaws (double t, void *params); +static double fn_qaws_L (double x, void *params); +static double fn_qaws_R (double x, void *params); + +static void +compute_result (const double * r, const double * cheb12, const double * cheb24, + double * result12, double * result24); + + +static void +qc25s (gsl_function * f, double a, double b, double a1, double b1, + gsl_integration_qaws_table * t, + double *result, double *abserr, int *err_reliable); + +static void +qc25s (gsl_function * f, double a, double b, double a1, double b1, + gsl_integration_qaws_table * t, + double *result, double *abserr, int *err_reliable) +{ + gsl_function weighted_function; + struct fn_qaws_params fn_params; + + fn_params.function = f; + fn_params.a = a; + fn_params.b = b; + fn_params.table = t; + + weighted_function.params = &fn_params; + + if (a1 == a && (t->alpha != 0.0 || t->mu != 0)) + { + double cheb12[13], cheb24[25]; + + double factor = pow(0.5 * (b1 - a1), t->alpha + 1.0); + + weighted_function.function = &fn_qaws_R; + + gsl_integration_qcheb (&weighted_function, a1, b1, cheb12, cheb24); + + if (t->mu == 0) + { + double res12 = 0, res24 = 0; + double u = factor; + + compute_result (t->ri, cheb12, cheb24, &res12, &res24); + + *result = u * res24; + *abserr = fabs(u * (res24 - res12)); + } + else + { + double res12a = 0, res24a = 0; + double res12b = 0, res24b = 0; + + double u = factor * log(b1 - a1); + double v = factor; + + compute_result (t->ri, cheb12, cheb24, &res12a, &res24a); + compute_result (t->rg, cheb12, cheb24, &res12b, &res24b); + + *result = u * res24a + v * res24b; + *abserr = fabs(u * (res24a - res12a)) + fabs(v * (res24b - res12b)); + } + + *err_reliable = 0; + + return; + } + else if (b1 == b && (t->beta != 0.0 || t->nu != 0)) + { + double cheb12[13], cheb24[25]; + double factor = pow(0.5 * (b1 - a1), t->beta + 1.0); + + weighted_function.function = &fn_qaws_L; + + gsl_integration_qcheb (&weighted_function, a1, b1, cheb12, cheb24); + + if (t->nu == 0) + { + double res12 = 0, res24 = 0; + double u = factor; + + compute_result (t->rj, cheb12, cheb24, &res12, &res24); + + *result = u * res24; + *abserr = fabs(u * (res24 - res12)); + } + else + { + double res12a = 0, res24a = 0; + double res12b = 0, res24b = 0; + + double u = factor * log(b1 - a1); + double v = factor; + + compute_result (t->rj, cheb12, cheb24, &res12a, &res24a); + compute_result (t->rh, cheb12, cheb24, &res12b, &res24b); + + *result = u * res24a + v * res24b; + *abserr = fabs(u * (res24a - res12a)) + fabs(v * (res24b - res12b)); + } + + *err_reliable = 0; + + return; + } + else + { + double resabs, resasc; + + weighted_function.function = &fn_qaws; + + gsl_integration_qk15 (&weighted_function, a1, b1, result, abserr, + &resabs, &resasc); + + if (*abserr == resasc) + { + *err_reliable = 0; + } + else + { + *err_reliable = 1; + } + + return; + } + +} + +static double +fn_qaws (double x, void *params) +{ + struct fn_qaws_params *p = (struct fn_qaws_params *) params; + gsl_function *f = p->function; + gsl_integration_qaws_table *t = p->table; + + double factor = 1.0; + + if (t->alpha != 0.0) + factor *= pow(x - p->a, t->alpha); + + if (t->beta != 0.0) + factor *= pow(p->b - x, t->beta); + + if (t->mu == 1) + factor *= log(x - p->a); + + if (t->nu == 1) + factor *= log(p->b - x); + + return factor * GSL_FN_EVAL (f, x); +} + +static double +fn_qaws_L (double x, void *params) +{ + struct fn_qaws_params *p = (struct fn_qaws_params *) params; + gsl_function *f = p->function; + gsl_integration_qaws_table *t = p->table; + + double factor = 1.0; + + if (t->alpha != 0.0) + factor *= pow(x - p->a, t->alpha); + + if (t->mu == 1) + factor *= log(x - p->a); + + return factor * GSL_FN_EVAL (f, x); +} + +static double +fn_qaws_R (double x, void *params) +{ + struct fn_qaws_params *p = (struct fn_qaws_params *) params; + gsl_function *f = p->function; + gsl_integration_qaws_table *t = p->table; + + double factor = 1.0; + + if (t->beta != 0.0) + factor *= pow(p->b - x, t->beta); + + if (t->nu == 1) + factor *= log(p->b - x); + + return factor * GSL_FN_EVAL (f, x); +} + + +static void +compute_result (const double * r, const double * cheb12, const double * cheb24, + double * result12, double * result24) +{ + size_t i; + double res12 = 0; + double res24 = 0; + + for (i = 0; i < 13; i++) + { + res12 += r[i] * cheb12[i]; + } + + for (i = 0; i < 25; i++) + { + res24 += r[i] * cheb24[i]; + } + + *result12 = res12; + *result24 = res24; +} |