% prepend_to_slang_load_path (".");
% set_import_module_path (path_dirname(__FILE__) + "/.." + ":" + get_import_module_path());

require ("gslinteg");

private define func (x)
{
   % Int from 0->1 is exp(0) - exp(1)
   return exp(-abs(x));
}

private define polyfunc (x, parms)
{
   return x^parms[0];
}

private define funct1 (x, parms)
{
   variable m = parms[0];
   return x^m + 1.0;
}

define slsh_main ()
{
   variable s = integration_qng (&func, -5, 5, 1e-12, 1e-12);
   print (s);
   s = integration_qag (&func, -5, 5, 1e-12, 1e-12, 32, 6);
   print (s);
   s = integration_qags (&func, -5, 5, 1e-12, 1e-12, 32);
   print (s);
   s = integration_qagp (&func, [-5:5:#10], 1e-12, 1e-12, 32);
   print (s);
   s = integration_qagi (&func, 1e-12, 1e-12, 32);
   print (s);
   s = integration_qagiu (&func, -5, 1e-12, 1e-12, 32);
   print (s);
   s = integration_qagil (&func, 5, 1e-12, 1e-12, 32);
   print (s);
   s = integration_qawc (&func, -5, 5, 0, 1e-12, 1e-12, 32);
   print (s);

   variable gl = integration_glfixed_alloc (3);
   s = integration_glfixed (&polyfunc, {3}, 0, 1, gl);
   % Result is int_0^1 x^3 dx = 1^4/4 = 0.25
   print (s);
   print (gl);

   variable m = 10, n = 6;
   variable ft = integration_fixed_alloc ("hermite", n, 0, 1, 0, 0);
   s = integration_fixed (&funct1, {m}, ft);
   print (s);
   print (ft);

   variable q = integration_qaws_alloc (0, 0, 0, 0);
   s = integration_qaws (&func, -5, 5, 1e-12, 1e-12, 64, q);
   print (s);
   print (q);

   q = integration_qawo_alloc (0, 5.0, GSL_INTEG_COSINE, 5);
   s = integration_qawo (&func, 0, 1e-12, 1e-12, 64, q);
   print (s);
   print (q);
   s = integration_qawf (&func, 0, 1e-12, 64, q);
   print (s);

   s = integration_cquad (&func, -5, 5, 1e-12, 1e-12, 100);
   print (s);

   s = integration_romberg (&func, -5, 5, 1e-12, 1e-12, 10);
   print (s);

}


   define example_func (x, parms)
   {
      variable alpha = parms[0];
      return log (alpha*x)/sqrt(x);
   }
   define slsh_main ()
   {
      variable alpha = 1.0;
      variable parms = {alpha};
      variable res
         = integration_qags (&example_func, parms, 0, 1, 0, 1e-7, 1000);

      vmessage ("Result = %S", res.result);
      vmessage ("Function calls: %S", res.neval);
      vmessage ("Status: %S", res.status);
      vmessage ("Abserror: %S", res.abserr);
   }
