123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239 |
- /* integration/qc25s.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_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;
- }
|