1 #include <gsl/gsl_integration.h>
2 #include <pygsl/error_helpers.h>
3 #include <pygsl/function_helpers.h>
4 #include <setjmp.h>
5
6 PyObject *module = NULL;
7 typedef struct{
8 PyObject * callback;
9 PyObject * args;
10 jmp_buf buffer;
11 }pygsl_function_args;
12
13
14 static double
pygsl_function_callback(double x,void * p)15 pygsl_function_callback(double x, void *p)
16 {
17 double value;
18 int flag;
19 pygsl_diff_args *pargs = NULL;
20
21 pargs = (pygsl_diff_args *) p;
22
23 assert(pargs->callback);
24 assert(pargs->args);
25 flag = PyGSL_function_wrap_helper(x, &value, NULL, pargs->callback,
26 pargs->args, (char *)__FUNCTION__);
27
28 if(GSL_SUCCESS != flag){
29 longjmp(pargs->buffer, flag);
30 return gsl_nan();
31 }
32 return value;
33 }
34
35
36 PyObject *
integrate(PyObject * self,PyObject * pyargs)37 integrate(PyObject *self, PyObject *pyargs)
38 {
39 PyObject *callback = NULL, *args = NULL;
40 pygsl_function_args cb_args;
41 gsl_function f;
42 int flag;
43 double a, b, epsabs, epsrel, result, abserr;
44 size_t neval;
45
46 if(!PyArg_ParseTuple(args, "Odddd|O", &callback, &a, &b, &epsabs, &epsrel, &args))
47 goto fail;
48
49 if(!Py_Callable(callback)){
50 pygsl_error("First argument must be a callable object!", GSL_EINVAL);
51 goto fail;
52 }
53
54 if(args == NULL){
55 args = Py_NONE;
56 }
57
58 f.f = pygsl_function_callback;
59 f.params = &cb_args;
60 cb_args.args = args;
61 cb_args.callback = callback;
62 flag = gsl_integration_qng (const *F, a, b, epsabs, epsrel, &result, &abserr, &neval)
63
64
65 if(PyGSL_ERROR_FLAG(flag) != GSL_SUCCESS)
66 goto fail;
67
68 return Py_BuildValue("(ddi)", result, abserr, neval);
69
70 fail:
71 Py_XDECREF(callback);
72 Py_XDECREF(args);
73 return NULL;
74 }
75 void
initintegrate(void)76 initintegrate(void)
77 {
78 PyObject *m;
79
80 m=Py_InitModule("integrate", mMethods);
81 module = m;
82 assert(m);
83 }
84