===================================================================
RCS file: /home/cvs/OpenXM/src/ox_gsl/call_gsl.c,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -p -r1.2 -r1.3
--- OpenXM/src/ox_gsl/call_gsl.c	2018/03/30 04:43:16	1.2
+++ OpenXM/src/ox_gsl/call_gsl.c	2018/04/17 00:56:38	1.3
@@ -1,4 +1,4 @@
-/* $OpenXM$ 
+/* $OpenXM: OpenXM/src/ox_gsl/call_gsl.c,v 1.2 2018/03/30 04:43:16 takayama Exp $ 
 */
 //#include <gsl/gsl_types.h>
 //#include <gsl/gsl_sys.h>
@@ -6,12 +6,13 @@
 #include <gsl/gsl_sf_result.h>
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_sf_gamma.h>
+#include <gsl/gsl_integration.h>
 #include "ox_gsl.h"
 extern int Debug;
 // local prototype declarations
 
 void  call_gsl_sf_lngamma_complex_e() {
-  cmo *c;
+  int argc;
   double zr;
   double zi;
   gsl_sf_result lnr;
@@ -21,7 +22,12 @@ void  call_gsl_sf_lngamma_complex_e() {
   cmo *ans;
   //  gsl_set_error_handler_off();
   gsl_set_error_handler((gsl_error_handler_t *)myhandler);
-  c = pop(); // number of args
+  argc = get_i(); // number of args
+  if (argc != 2) {
+    pops(argc);
+    push(make_error2("The argc must be 2 for gsl_sf_lngamma_complex_e.",NULL,0,-1));
+    return;
+  }
   zr=get_double();
   zi=get_double();
   if (Debug) printf("gsl_sf_lngamma_complex_e(zr=%lg,zi=%lg)\n",zr,zi);
@@ -32,3 +38,55 @@ void  call_gsl_sf_lngamma_complex_e() {
   ans = (cmo *)new_cmo_list_array((void *)r,3);
   push(ans);
 }
+
+cmo *Func_x=NULL;
+double f_x(double x,void *params) {
+  double d;
+  replace(1,"x",x);
+  if (eval_cmo(Func_x,&d)==0) GSL_ERROR("eval_cmo fails in f_x",GSL_ETOL);
+  return(d);
+}
+void  call_gsl_integration_qags() {
+  int argc;
+  double a;
+  double b;
+  int status;
+  cmo *r[3];
+  cmo *ans;
+  gsl_integration_workspace * w 
+    = gsl_integration_workspace_alloc (1000);
+  double result, error;
+  gsl_function F;
+  double epsabs=0;
+  double epsrel=1e-7;
+  int limit=1000;
+
+  //  gsl_set_error_handler_off();
+  gsl_set_error_handler((gsl_error_handler_t *)myhandler);
+  argc = get_i(); // number of args
+  if (argc != 3) {
+    pops(argc);
+    push(make_error2("The argc must be 3 for gsl_integration_qags.",NULL,0,-1));
+    return;
+  }
+  Func_x = pop();
+  a = get_double();
+  b = get_double();
+
+  F.function = &f_x;
+
+  status=gsl_integration_qags (&F, a, b, epsabs, epsrel, limit,
+                               w, &result, &error); 
+
+//  printf ("result          = % .18f\n", result);
+//  printf ("estimated error = % .18f\n", error);
+//  printf ("intervals       = %zu\n", w->size);
+
+  gsl_integration_workspace_free(w);
+
+  r[0] = (cmo *)new_cmo_double(result);
+  r[1] = (cmo *)new_cmo_double(error);
+  r[2] = (cmo *)new_cmo_int32(status);
+  ans = (cmo *)new_cmo_list_array((void *)r,3);
+  push(ans);
+}