pp_addpm({At=>Top},<<'EOD');
=head1 NAME

PDL::GSL::INTEG - PDL interface to numerical integration routines in GSL

=head1 DESCRIPTION

This is an interface to the numerical integration package present in the 
GNU Scientific Library, which is an implementation of QUADPACK.

Functions are named B<gslinteg_{algorithm}> where {algorithm}  
is the QUADPACK naming convention. The available functions are:

=over 3

=item gslinteg_qng: Non-adaptive Gauss-Kronrod integration

=item gslinteg_qag: Adaptive integration

=item gslinteg_qags: Adaptive integration with singularities

=item gslinteg_qagp: Adaptive integration with known singular points

=item gslinteg_qagi: Adaptive integration on infinite interval of the form (-\infty,\infty)

=item gslinteg_qagiu: Adaptive integration on infinite interval of the form (a,\infty)

=item gslinteg_qagil: Adaptive integration on infinite interval of the form (-\infty,b)

=item gslinteg_qawc: Adaptive integration for Cauchy principal values

=item gslinteg_qaws: Adaptive integration for singular functions

=item gslinteg_qawo: Adaptive integration for oscillatory functions

=item gslinteg_qawf: Adaptive integration for Fourier integrals

=back

Each algorithm computes an approximation to the integral, I, 
of the function f(x)w(x), where w(x) is a weight function 
(for general integrands w(x)=1). The user provides absolute
and relative error bounds (epsabs,epsrel) which specify
the following accuracy requirement:

|RESULT - I|  <= max(epsabs, epsrel |I|)


The routines will fail to converge if the 
error bounds are too stringent, but always return the best 
approximation obtained up to that stage

All functions return the result, and estimate of the absolute
error and an error flag (which is zero if there were no problems). 
You are responsible for checking for any errors, no warnings are issued
unless the option {Warn => 'y'} is specified in which case
the reason of failure will be printed.

You can nest integrals up to 20 levels. If you find yourself in
the unlikely situation that you need more, you can change the value
of 'max_nested_integrals' in the first line of the file 'FUNC.c' 
and recompile.

=for ref

Please check the GSL documentation for more information.

=head1 SYNOPSIS

   use PDL;
   use PDL::GSL::INTEG;

   my $a = 1.2;
   my $b = 3.7;
   my $epsrel = 0;
   my $epsabs = 1e-6;

   # Non adaptive integration
   my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&myf,$a,$b,$epsrel,$epsabs);
   # Warnings on
   my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&myf,$a,$b,$epsrel,$epsabs,{Warn=>'y'});

   # Adaptive integration with warnings on
   my $limit = 1000;
   my $key = 5;
   my ($res,$abserr,$ierr) = gslinteg_qag(\&myf,$a,$b,$epsrel,
                                     $epsabs,$limit,$key,{Warn=>'y'});

   sub myf{
     my ($x) = @_;
     return exp(-$x**2);
   }

=head1 FUNCTIONS

=head2 gslinteg_qng() -- Non-adaptive Gauss-Kronrod integration

This function applies the Gauss-Kronrod 10-point, 21-point, 43-point and 87-point 
integration rules in succession until an estimate of the integral of f over ($a,$b) 
is achieved within the desired absolute and relative error limits, $epsabs and $epsrel.
It is meant for fast integration of smooth functions. It returns an array with the 
result, an estimate of the absolute error, an error flag and the number of function
evaluations performed.

=for usage

Usage:

  ($res,$abserr,$ierr,$neval) = gslinteg_qng($function_ref,$a,$b,
                                             $epsrel,$epsabs,[{Warn => $warn}]);

=for example

Example:

   my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&f,0,1,0,1e-9);
   # with warnings on
   my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&f,0,1,0,1e-9,{Warn => 'y'});

   sub f{
     my ($x) = @_;
     return ($x**2.6)*log(1.0/$x);
   }


=head2 gslinteg_qag() -- Adaptive integration

This function applies an integration rule adaptively until an estimate of 
the integral of f over ($a,$b) is achieved within the desired absolute and 
relative error limits, $epsabs and $epsrel. On each iteration the adaptive 
integration strategy bisects the interval with the largest error estimate; 
the maximum number of allowed subdivisions is given by the parameter $limit.
The integration rule is determined by the 
value of $key, which has to be one of (1,2,3,4,5,6) and correspond to 
the 15, 21, 31, 41, 51 and 61  point Gauss-Kronrod rules respectively.
It returns an array with the result, an estimate of the absolute error 
and an error flag.

=for ref

Please check the GSL documentation for more information.

=for usage

Usage:

  ($res,$abserr,$ierr) = gslinteg_qag($function_ref,$a,$b,$epsrel,
                                      $epsabs,$limit,$key,[{Warn => $warn}]);

=for example

Example:
  
  my ($res,$abserr,$ierr) = gslinteg_qag(\&f,0,1,0,1e-10,1000,1);
  # with warnings on
  my ($res,$abserr,$ierr) = gslinteg_qag(\&f,0,1,0,1e-10,1000,1,{Warn => 'y'});

  sub f{
     my ($x) = @_;
     return ($x**2.6)*log(1.0/$x);
   }

=head2 gslinteg_qags() -- Adaptive integration with singularities

This function applies the Gauss-Kronrod 21-point integration rule 
adaptively until an estimate of the integral of f over ($a,$b) is 
achieved within the desired absolute and relative error limits, 
$epsabs and $epsrel. The algorithm is such that it
accelerates the convergence of the integral in the presence of 
discontinuities and integrable singularities. 
The maximum number of allowed subdivisions done by the adaptive
algorithm must be supplied in the parameter $limit.

=for ref

Please check the GSL documentation for more information.

=for usage

Usage:

  ($res,$abserr,$ierr) = gslinteg_qags($function_ref,$a,$b,$epsrel,
                                       $epsabs,$limit,[{Warn => $warn}]);

=for example

Example:

  my ($res,$abserr,$ierr) = gslinteg_qags(\&f,0,1,0,1e-10,1000);
  # with warnings on 
  ($res,$abserr,$ierr) = gslinteg_qags(\&f,0,1,0,1e-10,1000,{Warn => 'y'});

  sub f{
     my ($x) = @_;
     return ($x)*log(1.0/$x);
   }

=head2 gslinteg_qagp() -- Adaptive integration with known singular points

This function applies the adaptive integration algorithm used by 
gslinteg_qags taking into account the location of singular points
until an estimate of 
the integral of f over ($a,$b) is achieved within the desired absolute and 
relative error limits, $epsabs and $epsrel.
Singular points are supplied in the piddle $points, whose endpoints 
determine the integration range.
So, for example, if the function has singular points at x_1 and x_2 and the 
integral is desired from a to b (a < x_1 < x_2 < b), $points = pdl(a,x_1,x_2,b).
The maximum number of allowed subdivisions done by the adaptive
algorithm must be supplied in the parameter $limit.

=for ref

Please check the GSL documentation for more information.

=for usage

Usage:

  ($res,$abserr,$ierr) = gslinteg_qagp($function_ref,$points,$epsabs,
                                       $epsrel,$limit,[{Warn => $warn}])

=for example

Example:

  my $points = pdl(0,1,sqrt(2),3);
  my ($res,$abserr,$ierr) = gslinteg_qagp(\&f,$points,0,1e-3,1000);
  # with warnings on
  ($res,$abserr,$ierr) = gslinteg_qagp(\&f,$points,0,1e-3,1000,{Warn => 'y'});

  sub f{
    my ($x) = @_;
    my $x2 = $x**2;
    my $x3 = $x**3;
    return $x3 * log(abs(($x2-1.0)*($x2-2.0)));
  }

=head2 gslinteg_qagi() -- Adaptive integration on infinite interval

This function estimates the integral of the function f over the
infinite interval (-\infty,+\infty) within the desired absolute and 
relative error limits, $epsabs and $epsrel.
After a transformation, the algorithm
of gslinteg_qags with a 15-point Gauss-Kronrod rule is used.
The maximum number of allowed subdivisions done by the adaptive
algorithm must be supplied in the parameter $limit.

=for ref

Please check the GSL documentation for more information.

=for usage

Usage:

  ($res,$abserr,$ierr) = gslinteg_qagi($function_ref,$epsabs,
                                       $epsrel,$limit,[{Warn => $warn}]);

=for example

Example:

  my ($res,$abserr,$ierr) = gslinteg_qagi(\&myfn,1e-7,0,1000);
  # with warnings on
  ($res,$abserr,$ierr) = gslinteg_qagi(\&myfn,1e-7,0,1000,{Warn => 'y'});

  sub myfn{    
    my ($x) = @_;
    return exp(-$x - $x*$x) ;
  }


=head2 gslinteg_qagiu() -- Adaptive integration on infinite interval

This function estimates the integral of the function f over the
infinite interval (a,+\infty) within the desired absolute and 
relative error limits, $epsabs and $epsrel.
After a transformation, the algorithm
of gslinteg_qags with a 15-point Gauss-Kronrod rule is used.
The maximum number of allowed subdivisions done by the adaptive
algorithm must be supplied in the parameter $limit.

=for ref

Please check the GSL documentation for more information.

=for usage

Usage:

  ($res,$abserr,$ierr) = gslinteg_qagiu($function_ref,$a,$epsabs,
                                        $epsrel,$limit,[{Warn => $warn}]);

=for example

Example:

  my $alfa = 1;
  my ($res,$abserr,$ierr) = gslinteg_qagiu(\&f,99.9,1e-7,0,1000);
  # with warnings on
  ($res,$abserr,$ierr) = gslinteg_qagiu(\&f,99.9,1e-7,0,1000,{Warn => 'y'});

  sub f{
    my ($x) = @_;
    if (($x==0) && ($alfa == 1)) {return 1;}
    if (($x==0) && ($alfa > 1)) {return 0;}
    return ($x**($alfa-1))/((1+10*$x)**2);
  }

=head2 gslinteg_qagil() -- Adaptive integration on infinite interval

This function estimates the integral of the function f over the
infinite interval (-\infty,b) within the desired absolute and 
relative error limits, $epsabs and $epsrel.
After a transformation, the algorithm
of gslinteg_qags with a 15-point Gauss-Kronrod rule is used.
The maximum number of allowed subdivisions done by the adaptive
algorithm must be supplied in the parameter $limit.

=for ref

Please check the GSL documentation for more information.

=for usage

Usage:

  ($res,$abserr,$ierr) = gslinteg_qagl($function_ref,$b,$epsabs,
                                       $epsrel,$limit,[{Warn => $warn}]);

=for example

Example:

  my ($res,$abserr,$ierr) = gslinteg_qagil(\&myfn,1.0,1e-7,0,1000);
  # with warnings on
  ($res,$abserr,$ierr) = gslinteg_qagil(\&myfn,1.0,1e-7,0,1000,{Warn => 'y'});

  sub myfn{
    my ($x) = @_;
    return exp($x);
  }

=head2 gslinteg_qawc() -- Adaptive integration for Cauchy principal values

This function computes the Cauchy principal value of the integral of f over (a,b), 
with a singularity at c, I = \int_a^b dx f(x)/(x - c). The integral is
estimated within the desired absolute and relative error limits, $epsabs and $epsrel.
The maximum number of allowed subdivisions done by the adaptive
algorithm must be supplied in the parameter $limit.

=for ref

Please check the GSL documentation for more information.

=for usage

Usage:

  ($res,$abserr,$ierr) = gslinteg_qawc($function_ref,$a,$b,$c,$epsabs,$epsrel,$limit)

=for example

Example:

  my ($res,$abserr,$ierr) = gslinteg_qawc(\&f,-1,5,0,0,1e-3,1000);
  # with warnings on
  ($res,$abserr,$ierr) = gslinteg_qawc(\&f,-1,5,0,0,1e-3,1000,{Warn => 'y'});

  sub f{
    my ($x) = @_;
    return 1.0 / (5.0 * $x * $x * $x + 6.0) ;
  }

=head2 gslinteg_qaws() -- Adaptive integration for singular functions

The algorithm in gslinteg_qaws is designed for integrands with algebraic-logarithmic 
singularities at the end-points of an integration region. 
Specifically, this function computes the integral given by
I = \int_a^b dx f(x) (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x).
The integral is
estimated within the desired absolute and relative error limits, $epsabs and $epsrel.
The maximum number of allowed subdivisions done by the adaptive
algorithm must be supplied in the parameter $limit.

=for ref

Please check the GSL documentation for more information.

=for usage

Usage: 

  ($res,$abserr,$ierr) = 
      gslinteg_qawc($function_ref,$alpha,$beta,$mu,$nu,$a,$b,
                    $epsabs,$epsrel,$limit,[{Warn => $warn}]);

=for example

Example:

  my ($res,$abserr,$ierr) = gslinteg_qaws(\&f,0,0,1,0,0,1,0,1e-7,1000);
  # with warnings on
  ($res,$abserr,$ierr) = gslinteg_qaws(\&f,0,0,1,0,0,1,0,1e-7,1000,{Warn => 'y'});

  sub f{
    my ($x) = @_;
    if($x==0){return 0;}
    else{
      my $u = log($x);
      my $v = 1 + $u*$u;
      return 1.0/($v*$v);
    }
  }

=head2 gslinteg_qawo() -- Adaptive integration for oscillatory functions

This function uses an adaptive algorithm to compute the integral of f over 
(a,b) with the weight function sin(omega*x) or cos(omega*x) -- which of 
sine or cosine is used is determined by the parameter $opt ('cos' or 'sin').
The integral is
estimated within the desired absolute and relative error limits, $epsabs and $epsrel.
The maximum number of allowed subdivisions done by the adaptive
algorithm must be supplied in the parameter $limit.

=for ref

Please check the GSL documentation for more information.

=for usage

Usage:

  ($res,$abserr,$ierr) = gslinteg_qawo($function_ref,$omega,$sin_or_cos,
                                $a,$b,$epsabs,$epsrel,$limit,[opt])

=for example

Example:

  my $PI = 3.14159265358979323846264338328;
  my ($res,$abserr,$ierr) = PDL::GSL::INTEG::gslinteg_qawo(\&f,10*$PI,'sin',0,1,0,1e-7,1000);
  # with warnings on
  ($res,$abserr,$ierr) = PDL::GSL::INTEG::gslinteg_qawo(\&f,10*$PI,'sin',0,1,0,1e-7,1000,{Warn => 'y'});

  sub f{
    my ($x) = @_;
    if($x==0){return 0;}
    else{ return log($x);} 
  }


=head2 gslinteg_qawf() -- Adaptive integration for Fourier integrals

This function attempts to compute a Fourier integral of the function 
f over the semi-infinite interval [a,+\infty). Specifically, it attempts
tp compute I = \int_a^{+\infty} dx f(x)w(x), where w(x) is sin(omega*x)
or cos(omega*x) -- which of sine or cosine is used is determined by 
the parameter $opt ('cos' or 'sin').
The integral is
estimated within the desired absolute error limit $epsabs.
The maximum number of allowed subdivisions done by the adaptive
algorithm must be supplied in the parameter $limit.

=for ref

Please check the GSL documentation for more information.

=for usage

Usage:

  gslinteg_qawf($function_ref,$omega,$sin_or_cos,$a,$epsabs,$limit,[opt])

=for example

Example:

  my ($res,$abserr,$ierr) = gslinteg_qawf(\&f,$PI/2.0,'cos',0,1e-7,1000);
  # with warnings on
  ($res,$abserr,$ierr) = gslinteg_qawf(\&f,$PI/2.0,'cos',0,1e-7,1000,{Warn => 'y'});

  sub f{
    my ($x) = @_;
    if ($x == 0){return 0;}
    return 1.0/sqrt($x)    
  }


=head1 BUGS

Feedback is welcome. Log bugs in the PDL bug database (the
database is always linked from L<http://pdl.perl.org>).

=head1 SEE ALSO

L<PDL>

The GSL documentation is online at

  http://www.gnu.org/software/gsl/manual/

=head1 AUTHOR

This file copyright (C) 2003,2005 Andres Jordan <ajordan@eso.org>
All rights reserved. There is no warranty. You are allowed to redistribute 
this software documentation under certain conditions. For details, see the file
COPYING in the PDL distribution. If this file is separated from the
PDL distribution, the copyright notice should be included in the file.

The GSL integration routines were written by Brian Gough. QUADPACK
was written by Piessens, Doncker-Kapenga, Uberhuber and Kahaner.

=cut


EOD

pp_add_exported('','gslinteg_qng gslinteg_qag gslinteg_qags gslinteg_qagp
	            gslinteg_qagi gslinteg_qagiu gslinteg_qagil gslinteg_qawc
 		    gslinteg_qaws gslinteg_qawo gslinteg_qawf');

pp_addhdr('
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
/* GSL Glue code contains it\'s own "warn" variable and doesn\'t call perl\'s warn function  */
/*    (which can cause problems when perl\'s warn is called in a separate pthread),          */
/*    so we undefine the warn redefinition                                                   */
#undef warn

#include "FUNC.c"

void my_handler (const char * reason, 
                 const char * file, 
                 int line, 
                 int gsl_errno){
printf("Warning: %s at line %d of GSL file %s\n",reason,line,file);
}

');     

pp_addpm('
sub gslinteg_qng{
  my ($opt,$warn);
  if (ref($_[$#_]) eq \'HASH\'){ $opt = pop @_; }
  else{ $opt = {Warn => \'n\'}; }
  if($$opt{Warn}=~/y/i) { $warn = 1;}
  else {$warn = 0;} 
  my ($f,$a,$b,$epsabs,$epsrel) = @_;
  barf \'Usage: gslinteg_qng($function_ref,$a,$b,$epsabs,$epsrel,[opt]) \' 
	unless ($#_ == 4);
  my ($res,$abserr,$neval,$ierr) = qng_meat($a,$b,$epsabs,$epsrel,$warn,$f);
  return ($res,$abserr,$ierr,$neval);
}
');
pp_def('qng_meat',
        Pars => 'double a(); double b(); double epsabs();
                   double epsrel(); double [o] result(); double [o] abserr(); 
                   int [o] neval(); int [o] ierr(); int warn()',
	OtherPars => 'SV* funcion;',
        Docs => undef,
        Code => '
gsl_error_handler_t * old_handler;
if ($warn() == 1) { old_handler = gsl_set_error_handler(&my_handler); }
else { gsl_set_error_handler_off ();}
current_fun++;
if (current_fun >= (max_nested_integrals)) barf("Too many nested integrals, sorry!\n");
ext_funname[current_fun] = $COMP(funcion);
F.function = &FUNC;
F.params = 0;
$ierr() = gsl_integration_qng(&F,$a(),$b(),$epsabs(),$epsrel(),$P(result),$P(abserr),(size_t *) $P(neval));
current_fun--;
');  

pp_addpm('
sub gslinteg_qag{
   my ($opt,$warn);
   if (ref($_[$#_]) eq \'HASH\'){ $opt = pop @_; }
   else{ $opt = {Warn => \'n\'}; }
   if($$opt{Warn}=~/y/i) { $warn = 1;}
   else {$warn = 0;} 
   my ($f,$a,$b,$epsabs,$epsrel,$limit,$key) = @_;
   barf \'Usage: gslinteg_qag($function_ref,$a,$b,$epsabs,$epsrel,$limit,$key,[opt]) \' 
	unless ($#_ == 6);
   my ($res,$abserr,$ierr) = qag_meat($a,$b,$epsabs,$epsrel,$limit,$key,$limit,$warn,$f);
   return ($res,$abserr,$ierr);
}
');

pp_def('qag_meat',
        Pars => 'double a(); double b(); double epsabs();double epsrel(); int limit();
	           int key(); double [o] result(); double [o] abserr();int n();int [o] ierr();int warn();',
        OtherPars => 'SV* funcion;',
        Docs => undef,
        Code =>'
gsl_error_handler_t * old_handler;
if ($warn() == 1) { old_handler = gsl_set_error_handler(&my_handler); }
else { gsl_set_error_handler_off ();}
{gsl_integration_workspace *w;
current_fun++;
if (current_fun >= (max_nested_integrals)) barf("Too many nested integrals, sorry!\n");
ext_funname[current_fun] = $COMP(funcion);
F.function = &FUNC;
F.params = 0;
w = gsl_integration_workspace_alloc((size_t) $n());
$ierr() = gsl_integration_qag(&F,$a(),$b(),$epsabs(),$epsrel(),(size_t) $limit(),$key(),w,$P(result),$P(abserr));
gsl_integration_workspace_free(w); 
current_fun--;
}
');  

pp_addpm('
sub gslinteg_qags{
  my ($opt,$warn);
  if (ref($_[$#_]) eq \'HASH\'){ $opt = pop @_; }
  else{ $opt = {Warn => \'n\'}; }
  if($$opt{Warn}=~/y/i) { $warn = 1;}
  else {$warn = 0;} 
  my ($f,$a,$b,$epsabs,$epsrel,$limit) = @_;
  barf \'Usage: gslinteg_qags($function_ref,$a,$b,$epsabs,$epsrel,$limit,[opt]) \' 
	unless ($#_ == 5);
  my ($res,$abserr,$ierr) = qags_meat($a,$b,$epsabs,$epsrel,$limit,$limit,$warn,$f);
  return ($res,$abserr,$ierr);
}
');
pp_def('qags_meat',
        Pars => 'double a(); double b(); double epsabs();double epsrel(); int limit();
	           double [o] result(); double [o] abserr();int n();int [o] ierr();int warn();',
        OtherPars => 'SV* funcion;',
        Docs => undef,
        Code =>'
gsl_error_handler_t * old_handler;
if ($warn() == 1) { old_handler = gsl_set_error_handler(&my_handler); }
else { gsl_set_error_handler_off ();}
{gsl_integration_workspace *w;
current_fun++;
if (current_fun >= (max_nested_integrals)) barf("Too many nested integrals, sorry!\n");
ext_funname[current_fun] = $COMP(funcion);
F.function = &FUNC;
F.params = 0;
w = gsl_integration_workspace_alloc((size_t) $n());
$ierr() = gsl_integration_qags(&F,$a(),$b(),$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr));
gsl_integration_workspace_free(w);
current_fun--;
}
');  

pp_addpm('
sub gslinteg_qagp{
  my ($opt,$warn);
  if (ref($_[$#_]) eq \'HASH\'){ $opt = pop @_; }
  else{ $opt = {Warn => \'n\'}; }
  if($$opt{Warn}=~/y/i) { $warn = 1;}
  else {$warn = 0;} 
  my ($f,$points,$epsabs,$epsrel,$limit) = @_;
  barf \'Usage: gslinteg_qagp($function_ref,$points,$epsabs,$epsrel,$limit,[opt]) \' 
	unless ($#_ == 4);
  my ($res,$abserr,$ierr) = qagp_meat($points,$epsabs,$epsrel,$limit,$limit,$warn,$f);
  return ($res,$abserr,$ierr);
}
');
pp_def('qagp_meat',
        Pars => 'double pts(l); double epsabs();double epsrel();int limit();
		   double [o] result(); double [o] abserr();int n();int [o] ierr();int warn();',
	OtherPars => 'SV* funcion;',
	Code =>'
gsl_error_handler_t * old_handler;
if ($warn() == 1) { old_handler = gsl_set_error_handler(&my_handler); }
else { gsl_set_error_handler_off ();}
{gsl_integration_workspace *w;
current_fun++;
if (current_fun >= (max_nested_integrals)) barf("Too many nested integrals, sorry!\n");
ext_funname[current_fun] = $COMP(funcion);
F.function = &FUNC;
F.params = 0;
w = gsl_integration_workspace_alloc((size_t) $n());
$ierr() = gsl_integration_qagp(&F,$P(pts),(size_t) $SIZE(l),$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr));
gsl_integration_workspace_free(w);
current_fun--;
}
');  

pp_addpm('
sub gslinteg_qagi{
  my ($opt,$warn);
  if (ref($_[$#_]) eq \'HASH\'){ $opt = pop @_; }
  else{ $opt = {Warn => \'n\'}; }
  if($$opt{Warn}=~/y/i) { $warn = 1;}
  else {$warn = 0;} 
  my ($f,$epsabs,$epsrel,$limit) = @_;
  barf \'Usage: gslinteg_qagi($function_ref,$epsabs,$epsrel,$limit,[opt]) \' 
	unless ($#_ == 3);
  my ($res,$abserr,$ierr) = qagi_meat($epsabs,$epsrel,$limit,$limit,$warn,$f);
  return ($res,$abserr,$ierr);
}
');
pp_def('qagi_meat',
        Pars => 'double epsabs();double epsrel(); int limit();
		   double [o] result(); double [o] abserr(); int n(); int [o] ierr();int warn();',
        OtherPars => 'SV* funcion;',
	Code =>'
gsl_error_handler_t * old_handler;
if ($warn() == 1) { old_handler = gsl_set_error_handler(&my_handler); }
else { gsl_set_error_handler_off ();}
{gsl_integration_workspace *w;
current_fun++;
if (current_fun >= (max_nested_integrals)) barf("Too many nested integrals, sorry!\n");
ext_funname[current_fun] = $COMP(funcion);
F.function = &FUNC;
F.params = 0;
w = gsl_integration_workspace_alloc((size_t) $n());
$ierr() = gsl_integration_qagi(&F,$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr));
gsl_integration_workspace_free(w);
current_fun--;
}
');  

pp_addpm('
sub gslinteg_qagiu{
  my ($opt,$warn);
  if (ref($_[$#_]) eq \'HASH\'){ $opt = pop @_; }
  else{ $opt = {Warn => \'n\'}; }
  if($$opt{Warn}=~/y/i) { $warn = 1;}
  else {$warn = 0;} 
  my ($f,$a,$epsabs,$epsrel,$limit) = @_;
  barf \'Usage: gslinteg_qagiu($function_ref,$a,$epsabs,$epsrel,$limit,[opt]) \' 
	unless ($#_ == 4);
  my ($res,$abserr,$ierr) = qagiu_meat($a,$epsabs,$epsrel,$limit,$limit,$warn,$f);
  return ($res,$abserr,$ierr);
}
');
pp_def('qagiu_meat',
        Pars => 'double a(); double epsabs();double epsrel();int limit();
		   double [o] result(); double [o] abserr();int n();int [o] ierr();int warn();',
        OtherPars => 'SV* funcion;',
        Docs => undef,
	Code =>'
gsl_error_handler_t * old_handler;
if ($warn() == 1) { old_handler = gsl_set_error_handler(&my_handler); }
else { gsl_set_error_handler_off ();}
{gsl_integration_workspace *w;
current_fun++;
if (current_fun >= (max_nested_integrals)) barf("Too many nested integrals, sorry!\n");
ext_funname[current_fun] = $COMP(funcion);
F.function = &FUNC;
F.params = 0;
w = gsl_integration_workspace_alloc((size_t) $n());
$ierr() = gsl_integration_qagiu(&F,$a(),$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr));
gsl_integration_workspace_free(w);
current_fun--;
}
');  

pp_addpm('
sub gslinteg_qagil{
  my ($opt,$warn);
  if (ref($_[$#_]) eq \'HASH\'){ $opt = pop @_; }
  else{ $opt = {Warn => \'n\'}; }
  if($$opt{Warn}=~/y/i) { $warn = 1;}
  else {$warn = 0;} 
  my ($f,$b,$epsabs,$epsrel,$limit) = @_;
  barf \'Usage: gslinteg_qagil($function_ref,$b,$epsabs,$epsrel,$limit,[opt]) \' 
	unless ($#_ == 4);
  my ($res,$abserr,$ierr) = qagil_meat($b,$epsabs,$epsrel,$limit,$limit,$warn,$f);
  return ($res,$abserr,$ierr);
}
');
pp_def('qagil_meat',
        Pars => 'double b(); double epsabs();double epsrel();int limit();
		   double [o] result(); double [o] abserr();int n();int [o] ierr();int warn();',
        OtherPars => 'SV* funcion;',
        Docs => undef,
	Code =>'
gsl_error_handler_t * old_handler;
if ($warn() == 1) { old_handler = gsl_set_error_handler(&my_handler); }
else { gsl_set_error_handler_off ();}
{gsl_integration_workspace *w;
current_fun++;
if (current_fun >= (max_nested_integrals)) barf("Too many nested integrals, sorry!\n");
ext_funname[current_fun] = $COMP(funcion);
F.function = &FUNC;
F.params = 0;
w = gsl_integration_workspace_alloc((size_t) $n());
$ierr() = gsl_integration_qagil(&F,$b(),$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr));
gsl_integration_workspace_free(w);
current_fun--;
}
');  

pp_addpm('
sub gslinteg_qawc{
  my ($opt,$warn);
  if (ref($_[$#_]) eq \'HASH\'){ $opt = pop @_; }
  else{ $opt = {Warn => \'n\'}; }
  if($$opt{Warn}=~/y/i) { $warn = 1;}
  else {$warn = 0;} 
  my ($f,$a,$b,$c,$epsabs,$epsrel,$limit) = @_;
  barf \'Usage: gslinteg_qawc($function_ref,$a,$b,$c,$epsabs,$epsrel,$limit,[opt]) \' 
	unless ($#_ == 6);
  my ($res,$abserr,$ierr) = qawc_meat($a,$b,$c,$epsabs,$epsrel,$limit,$limit,$warn,$f);
  return ($res,$abserr,$ierr);
}
');
pp_def('qawc_meat',
        Pars => 'double a(); double b(); double c(); double epsabs();double epsrel();int limit();
	           double [o] result(); double [o] abserr();int n();int [o] ierr();int warn();',
	OtherPars => 'SV* funcion;',
	Code =>'
gsl_error_handler_t * old_handler;
if ($warn() == 1) { old_handler = gsl_set_error_handler(&my_handler); }
else { gsl_set_error_handler_off ();}
{gsl_integration_workspace *w;
current_fun++;
if (current_fun >= (max_nested_integrals)) barf("Too many nested integrals, sorry!\n");
ext_funname[current_fun] = $COMP(funcion);
F.function = &FUNC;
F.params = 0;
w = gsl_integration_workspace_alloc((size_t) $n());
$ierr() = gsl_integration_qawc(&F,$a(),$b(),$c(),$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr));
gsl_integration_workspace_free(w);
current_fun--;
}
');  

pp_addpm('
sub gslinteg_qaws{
  my ($opt,$warn);
  if (ref($_[$#_]) eq \'HASH\'){ $opt = pop @_; }
  else{ $opt = {Warn => \'n\'}; }
  if($$opt{Warn}=~/y/i) { $warn = 1;}
  else {$warn = 0;} 
  my ($f,$alpha,$beta,$mu,$nu,$a,$b,$epsabs,$epsrel,$limit) = @_;
  barf \'Usage: gslinteg_qaws($function_ref,$alpha,$beta,$mu,$nu,$a,$b,$epsabs,$epsrel,$limit,[opt]) \' 
	unless ($#_ == 9);
  my ($res,$abserr,$ierr) = qaws_meat($a,$b,$epsabs,$epsrel,$limit,$limit,$alpha,$beta,$mu,$nu,$warn,$f);
  return ($res,$abserr,$ierr);
}
');
pp_def('qaws_meat',
        Pars => 'double a(); double b();double epsabs();double epsrel();int limit();
	         double [o] result(); double [o] abserr();int n();
		 double alpha(); double beta(); int mu(); int nu();int [o] ierr();int warn();',
        OtherPars => 'SV* funcion;',
	Code => '
gsl_error_handler_t * old_handler;
if ($warn() == 1) { old_handler = gsl_set_error_handler(&my_handler); }
else { gsl_set_error_handler_off ();}
{gsl_integration_qaws_table * qtab;
gsl_integration_workspace *w;

qtab = gsl_integration_qaws_table_alloc($alpha(),$beta(),$mu(),$nu());

current_fun++;
if (current_fun >= (max_nested_integrals)) barf("Too many nested integrals, sorry!\n");
ext_funname[current_fun] = $COMP(funcion);
F.function = &FUNC;
F.params = 0;
w = gsl_integration_workspace_alloc((size_t) $n());
$ierr() = gsl_integration_qaws(&F,$a(),$b(),qtab,$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr));
gsl_integration_workspace_free(w);
gsl_integration_qaws_table_free(qtab);
current_fun--;
}
');  

pp_addpm('
sub gslinteg_qawo{
  my ($opt,$warn);
  if (ref($_[$#_]) eq \'HASH\'){ $opt = pop @_; }
  else{ $opt = {Warn => \'n\'}; }
  if($$opt{Warn}=~/y/i) { $warn = 1;}
  else {$warn = 0;} 
  my ($f,$omega,$sincosopt,$a,$b,$epsabs,$epsrel,$limit) = @_;
  barf \'Usage: gslinteg_qawo($function_ref,$omega,$sin_or_cos,$a,$b,$epsabs,$epsrel,$limit,[opt]) \' 
	unless ($#_ == 7);
  my $OPTION_SIN_COS;
  if($sincosopt=~/cos/i){ $OPTION_SIN_COS = 0;}
  elsif($sincosopt=~/sin/i){ $OPTION_SIN_COS = 1;}
  else { barf("Error in argument 3 of function gslinteg_qawo: specify \'cos\' or \'sin\'\n");}

  my $L = $b - $a;
  my $nlevels = $limit;
  my ($res,$abserr,$ierr) = qawo_meat($a,$b,$epsabs,$epsrel,$limit,$limit,$OPTION_SIN_COS,$omega,$L,$nlevels,$warn,$f);
  return ($res,$abserr,$ierr);
}
');
pp_def('qawo_meat',
	Pars => 'double a(); double b();double epsabs();double epsrel();int limit();
	         double [o] result(); double [o] abserr();int n();
		 int sincosopt(); double omega(); double L(); int nlevels();int [o] ierr();int warn();',
	OtherPars =>  'SV* funcion;',
	Code => '
gsl_error_handler_t * old_handler;
if ($warn() == 1) { old_handler = gsl_set_error_handler(&my_handler); }
else { gsl_set_error_handler_off ();}
{gsl_integration_qawo_table * qtab;
gsl_integration_workspace *w;
enum gsl_integration_qawo_enum T;

T = GSL_INTEG_SINE;
if ($sincosopt() == 0){ T = GSL_INTEG_COSINE ;}

qtab = gsl_integration_qawo_table_alloc($omega(),$L(),T,(size_t) $nlevels());

current_fun++;
if (current_fun >= (max_nested_integrals)) barf("Too many nested integrals, sorry!\n");
ext_funname[current_fun] = $COMP(funcion);
F.function = &FUNC;
F.params = 0;
w = gsl_integration_workspace_alloc((size_t) $n());
$ierr() = gsl_integration_qawo(&F,$a(),$epsabs(),$epsrel(),(size_t) $limit(),w,qtab,$P(result),$P(abserr));
gsl_integration_workspace_free(w);
gsl_integration_qawo_table_free(qtab);
current_fun--;
}
');  

pp_addpm('
sub gslinteg_qawf{
  my ($opt,$warn);
  if (ref($_[$#_]) eq \'HASH\'){ $opt = pop @_; }
  else{ $opt = {Warn => \'n\'}; }
  if($$opt{Warn}=~/y/i) { $warn = 1;}
  else {$warn = 0;} 
  my ($f,$omega,$sincosopt,$a,$epsabs,$limit) = @_;
  barf \'Usage: gslinteg_qawf($function_ref,$omega,$sin_or_cos,$a,$epsabs,$limit,[opt]) \' 
	unless ($#_ == 5);
  my $OPTION_SIN_COS;
  if($sincosopt=~/cos/i){ $OPTION_SIN_COS = 0;}
  elsif($sincosopt=~/sin/i){ $OPTION_SIN_COS = 1;}
  else { barf("Error in argument 3 of function gslinteg_qawf: specify \'cos\' or \'sin\'\n");}
  my $nlevels = $limit;
  my ($res,$abserr,$ierr) = qawf_meat($a,$epsabs,$limit,$limit,$OPTION_SIN_COS,$omega,$nlevels,$warn,$f);
  return ($res,$abserr,$ierr);
}
');

pp_def('qawf_meat',
	Pars => 'double a(); double epsabs();int limit();
		 double [o] result(); double [o] abserr();int n();
		 int sincosopt(); double omega(); int nlevels();int [o] ierr();int warn();',
	OtherPars => 'SV* funcion;',
	Code => '
gsl_error_handler_t * old_handler;
if ($warn() == 1) { old_handler = gsl_set_error_handler(&my_handler); }
else { gsl_set_error_handler_off ();}
{gsl_integration_qawo_table * qtab;
gsl_integration_workspace *w;
gsl_integration_workspace *cw;
enum gsl_integration_qawo_enum T;

T = GSL_INTEG_SINE;
if ($sincosopt() == 0){ T = GSL_INTEG_COSINE ;}

qtab = gsl_integration_qawo_table_alloc($omega(),1.,T,(size_t) $nlevels());

current_fun++;
if (current_fun >= (max_nested_integrals)) barf("Too many nested integrals, sorry!\n");
ext_funname[current_fun] = $COMP(funcion);
F.function = &FUNC;
F.params = 0;

w = gsl_integration_workspace_alloc((size_t) $n());
cw = gsl_integration_workspace_alloc((size_t) $n());

$ierr() = gsl_integration_qawf(&F,$a(),$epsabs(),(size_t) $limit(),w,cw,qtab,$P(result),$P(abserr));

gsl_integration_workspace_free(w);
gsl_integration_workspace_free(cw);
gsl_integration_qawo_table_free(qtab);
current_fun--;
}
');  


pp_done();