oldtests/gsl-intgamma/quadratures.c

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include <gsl/gsl_integration.h>
   4 #include <gsl/gsl_errno.h>
   5
   6 #ifndef EPSABS
   7 #define EPSABS 0
   8 #endif
   9 #ifndef EPSREL
  10 #define EPSREL 1e-13
  11 #endif
  12 #ifndef LIMIT
  13 #define LIMIT 30000 //???
  14 #endif
  15 #ifndef R0
  16 #define R0 8e-6
  17 #endif
  18
  19
  20 /* Relevant quadrature methods from gsl:
  21  * gsl_integration_qagiu
  22  * ... and that's probably it.
  23  */
  24
  25 //gsl_function sigma2_integrand;
  26
  27 struct sigma2_integrand_params {
  28   int n;
  29   double k, R;
  30 };
  31
  32 static inline double sq(double x) {return x * x;}
  33
  34 double sigma2_integrand(double ksi, void *params) {
  35   struct sigma2_integrand_params *p = (struct sigma2_integrand_params *) params;
  36   return exp(-sq(p->R*ksi) + sq(p->k/ksi/2)) * pow(ksi, 2*p->n);
  37 }
  38
  39
  40 int main(int argc, char **argv) {
  41   struct sigma2_integrand_params p;
  42   gsl_function F;
  43   F.function = sigma2_integrand;
  44   F.params = &p;
  45   gsl_integration_workspace *workspace = gsl_integration_workspace_alloc(LIMIT);
  46   gsl_error_handler_t * old_handler=gsl_set_error_handler_off();
  47   double eta, eta_orig, k_orig, R_orig;
  48   while (scanf("%d %lf %lf %lf", &(p.n), &k_orig, &R_orig, &eta_orig) == 4) {
  49     eta = eta_orig * R0;
  50     p.k = k_orig * R0;
  51     p.R = R_orig / R0;
  52     double result, abserr;
  53     int retval = gsl_integration_qagiu(&F, eta, EPSABS, EPSREL,
  54         LIMIT, workspace, &result, &abserr);
  55     double normfac = pow(R0, -2*p.n - 1);
  56     result *= normfac;
  57     abserr *= normfac;
  58     printf("%d %.16g %.16g %.16g %.16g %.16g %d\n", p.n, k_orig, R_orig,
  59         eta_orig, result, abserr, retval);
  60   }
  61   gsl_integration_workspace_free(workspace);
  62 }
  63
  64