qpms/tests/gsl-intgamma/quadratures.c

65 lines
1.5 KiB
C

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#ifndef EPSABS
#define EPSABS 0
#endif
#ifndef EPSREL
#define EPSREL 1e-13
#endif
#ifndef LIMIT
#define LIMIT 30000 //???
#endif
#ifndef R0
#define R0 8e-6
#endif
/* Relevant quadrature methods from gsl:
* gsl_integration_qagiu
* ... and that's probably it.
*/
//gsl_function sigma2_integrand;
struct sigma2_integrand_params {
int n;
double k, R;
};
static inline double sq(double x) {return x * x;}
double sigma2_integrand(double ksi, void *params) {
struct sigma2_integrand_params *p = (struct sigma2_integrand_params *) params;
return exp(-sq(p->R*ksi) + sq(p->k/ksi/2)) * pow(ksi, 2*p->n);
}
int main(int argc, char **argv) {
struct sigma2_integrand_params p;
gsl_function F;
F.function = sigma2_integrand;
F.params = &p;
gsl_integration_workspace *workspace = gsl_integration_workspace_alloc(LIMIT);
gsl_error_handler_t * old_handler=gsl_set_error_handler_off();
double eta, eta_orig, k_orig, R_orig;
while (scanf("%d %lf %lf %lf", &(p.n), &k_orig, &R_orig, &eta_orig) == 4) {
eta = eta_orig * R0;
p.k = k_orig * R0;
p.R = R_orig / R0;
double result, abserr;
int retval = gsl_integration_qagiu(&F, eta, EPSABS, EPSREL,
LIMIT, workspace, &result, &abserr);
double normfac = pow(R0, -2*p.n - 1);
result *= normfac;
abserr *= normfac;
printf("%d %.16g %.16g %.16g %.16g %.16g %d\n", p.n, k_orig, R_orig,
eta_orig, result, abserr, retval);
}
gsl_integration_workspace_free(workspace);
}