/*
 * author: Pierre Schnizer
 * created: June 2003
 */
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_nan.h>
#include <pygsl/block_helpers.h>
#include <Python.h>

#define IMPORTALL 1
#ifdef PyGSL_NUMERIC
#include <Numeric/arrayobject.h>
#include <Numeric/ufuncobject.h>
#else
#error numarray ufuncs not yet supported!
#ifdef IMPORTALL
#undef IMPORTALL
#endif
#endif /* PyGSL_NUMERIC */
#include <pygsl/error_helpers.h>
/*
   #include <pygsl/utils.h>
   #include <pygsl/error_helpers.h>
*/

static PyObject* gsl_module_error;
static PyObject* module = NULL; /*used by the backtrace*/



static void 
invoke_error_handler(const char *filename, const char *func_name, int errno, int element);


static void 
invoke_error_handler(const char *filename, const char *func_name, int errno, int element)
{
     ;
     /* DEBUG_MESS(2, "Got error flag of %d for element number %d", errno, element); */
}

/* I add the evaluate functions types by hand to have an extra check */
#ifdef IMPORTALL
#include "sf__arrays.c"

#define SF_ARRAY(name, function)                                           \
static PyObject* sf_ ## name      (PyObject *self, PyObject *args)         \
{                                                                          \
     PyObject * tmp;                                                       \
     FUNC_MESS_BEGIN();                                                    \
     tmp =  PyGSL_sf_array_evaluator_ ## function (self, args, gsl_sf_ ## name); \
     if (tmp == NULL){                                                     \
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__,  __LINE__);  \
     }                                                                     \
     FUNC_MESS_END();                                                      \
     return tmp;                                                           \
}

SF_ARRAY(bessel_Jn_array, iid_ad);
SF_ARRAY(bessel_Yn_array, iid_ad);
SF_ARRAY(bessel_In_array, iid_ad);
SF_ARRAY(bessel_Kn_array, iid_ad);
SF_ARRAY(bessel_Kn_scaled_array, iid_ad);
SF_ARRAY(bessel_jl_array, id_ad);
SF_ARRAY(bessel_jl_steed_array, id_ad);
SF_ARRAY(bessel_yl_array, id_ad);
SF_ARRAY(bessel_il_scaled_array, id_ad);
SF_ARRAY(bessel_kl_scaled_array, id_ad);

#undef SF_ARRAY
#endif

static PyMethodDef sf_functions[] = {
#ifdef IMPORTALL
     {"bessel_Jn_array", (PyCFunction) sf_bessel_Jn_array, METH_VARARGS, NULL},
     {"bessel_Yn_array", (PyCFunction) sf_bessel_Yn_array, METH_VARARGS, NULL},
     {"bessel_In_array", (PyCFunction) sf_bessel_In_array, METH_VARARGS, NULL},
     {"bessel_Kn_array", (PyCFunction) sf_bessel_Kn_array, METH_VARARGS, NULL},
     {"bessel_Kn_scaled_array", (PyCFunction) sf_bessel_Kn_scaled_array, METH_VARARGS, NULL},
     {"bessel_jl_array", (PyCFunction) sf_bessel_jl_array, METH_VARARGS, NULL},
     {"bessel_jl_steed_array", (PyCFunction) sf_bessel_jl_steed_array, METH_VARARGS, NULL},
     {"bessel_yl_array", (PyCFunction) sf_bessel_yl_array, METH_VARARGS, NULL},
     {"bessel_il_scaled_array", (PyCFunction) sf_bessel_il_scaled_array, METH_VARARGS, NULL},
     {"bessel_kl_scaled_array", (PyCFunction) sf_bessel_kl_scaled_array, METH_VARARGS, NULL},
#endif
     {NULL, NULL, 0}
};

#ifdef IMPORTALL
#include "sf__evals.c" 
#include "sf__data.c"
#endif

static void * polar_to_rect_data [] = { (void *) gsl_sf_polar_to_rect, (void *) gsl_sf_polar_to_rect};
static void * rect_to_polar_data [] = { (void *) gsl_sf_rect_to_polar, (void *) gsl_sf_rect_to_polar};

static char PyGSL_sf_ufunc_qi_dd_D_one_types [] = { PyArray_FLOAT,  PyArray_FLOAT,  PyArray_CFLOAT, PyArray_CFLOAT, 
                                                    PyArray_DOUBLE, PyArray_DOUBLE, PyArray_CDOUBLE,PyArray_CDOUBLE };

static char PyGSL_sf_ufunc_qi_D_dd_one_types [] = { PyArray_CFLOAT,  PyArray_FLOAT,  PyArray_FLOAT,  PyArray_FLOAT,  PyArray_FLOAT, 
                                                    PyArray_CDOUBLE, PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE };

typedef int PyGSL_sf_ufunc_qi_dd_D_one(double, double, gsl_sf_result *, gsl_sf_result *);
void  PyGSL_sf_ufunc_qi_dd_D(char **args, int *dimensions, int *steps, void *func)
{
        int is0=steps[0], is1=steps[1], os0=steps[2], os1=steps[3];
        char   *ip0=args[0], *ip1=args[1], *op0=args[2], *op1=args[3];
	double *dop, *derr;

	gsl_sf_result x, y;
	int i, tmp;
	FUNC_MESS_BEGIN();
	DEBUG_MESS(2, "dimensions = %d %d %d", dimensions[0], dimensions[1],dimensions[2]);
	DEBUG_MESS(2, "steps = %d %d %d", steps[0], steps[1], steps[2]);
	DEBUG_MESS(2, "args = %p %p %p", args[0], args[1], args[2]);
        for(i = 0; i<dimensions[0]; i++, ip0+=is0, ip1+=is1, op0+=os0,op1+=os1){
	     DEBUG_MESS(2, "i = %d", i);
	    dop = ((double *) op0);
	    derr = ((double *) op1);
	    tmp = ((PyGSL_sf_ufunc_qi_dd_D_one *) func)(*((double *)ip0), *((double *)ip1), &x, &y);
	    dop[0] = x.val;
	    dop[1] = y.val;
	    derr[0] = x.err;
	    derr[1] = y.err;
        }
	FUNC_MESS_END();
}
void  PyGSL_sf_ufunc_qi_dd_D_as_ff_F(char **args, int *dimensions, int *steps, void *func)
{
        int i,  is0=steps[0], is1=steps[1], os0=steps[2], os1=steps[3], tmp;
        char   *ip0=args[0], *ip1=args[1], *op0=args[2], *op1=args[3];
	gsl_sf_result x, y;

	FUNC_MESS_BEGIN();
	DEBUG_MESS(2, "dimensions = %d %d %d", dimensions[0], dimensions[1],dimensions[2]);
	DEBUG_MESS(2, "steps = %d %d %d", steps[0], steps[1], steps[2]);
	DEBUG_MESS(2, "args = %p %p %p", args[0], args[1], args[2]);
        for(i = 0; i<dimensions[0]; i++, ip0+=is0, ip1+=is1, op0+=os0, op1+=os1){	     
	     DEBUG_MESS(2, "i = %d", i);
	     tmp = ((PyGSL_sf_ufunc_qi_dd_D_one *) func)((double) *((float *)ip0), (double) *((float *)ip1), &x, &y);
	     *((float *)op0) = (float) x.val;
	     *(((float *)op0)+1) = (float) y.val;
	     *((float *)op1) = (float) x.err;
	     *(((float *)op1)+1) = (float) y.err;
        }
	FUNC_MESS_END();
}
static PyUFuncGenericFunction PyGSL_sf_ufunc_qi_dd_D_one_data[] = {PyGSL_sf_ufunc_qi_dd_D_as_ff_F, PyGSL_sf_ufunc_qi_dd_D};

typedef int PyGSL_sf_ufunc_qi_D_dd_one(double, double, gsl_sf_result *, gsl_sf_result *);
void  PyGSL_sf_ufunc_qi_D_dd(char **args, int *dimensions, int *steps, void *func)
{
        int i,  is0=steps[0], os0=steps[1], os1=steps[2],os2=steps[3],os3=steps[4], tmp;
        char   *ip0=args[0], *op0=args[1], *op1=args[2],*op2=args[3], *op3=args[4];
	double *dip;
	gsl_sf_result radius, angle;

	FUNC_MESS_BEGIN();
	DEBUG_MESS(2, "dimensions = %d %d %d", dimensions[0], dimensions[1],dimensions[2]);
	DEBUG_MESS(2, "steps = %d %d %d", steps[0], steps[1], steps[2]);
	DEBUG_MESS(2, "args = %p %p %p", args[0], args[1], args[2]);
	DEBUG_MESS(1, "rect_to_polar %p", (void *)func);
        for(i = 0; i<dimensions[0]; i++, ip0+=is0, op0+=os0, op1+=os1, op2+=os2, op3+=os3){
	     DEBUG_MESS(2, "i = %d", i);
	    dip = ((double *) ip0);
	    tmp = ((PyGSL_sf_ufunc_qi_D_dd_one *) func)(dip[0], dip[1], &radius, &angle);
	    *((double *)op0) = radius.val;
	    *((double *)op1) = radius.err;
	    *((double *)op2) = angle.val;
	    *((double *)op3) = angle.err;
        }
	FUNC_MESS_END();
}
void  PyGSL_sf_ufunc_qi_D_dd_as_F_ff(char **args, int *dimensions, int *steps, void *func)
{
        int i,  is0=steps[0], os0=steps[1], os1=steps[2],os2=steps[3],os3=steps[4], tmp;
        char   *ip0=args[0], *op0=args[1], *op1=args[2],*op2=args[3], *op3=args[4];
	gsl_sf_result radius, angle;
	FUNC_MESS_BEGIN();
	DEBUG_MESS(2, "dimensions = %d %d %d", dimensions[0], dimensions[1],dimensions[2]);
	DEBUG_MESS(2, "steps = %d %d %d", steps[0], steps[1], steps[2]);
	DEBUG_MESS(2, "args = %p %p %p", args[0], args[1], args[2]);
	DEBUG_MESS(1, "polar_to_rect %p", (void *)func);
        for(i = 0; i<dimensions[0]; i++, ip0+=is0, op0+=os0, op1+=os1, op2+=os2, op3+=os3){
	     tmp = ((PyGSL_sf_ufunc_qi_D_dd_one *) func)((double) *((float *)ip0), (double) *(((float *)ip0)+1), &radius, &angle);
	     DEBUG_MESS(2, "i = %d", i);
	    *((float *)op0) = (float) radius.val;
	    *((float *)op1) = (float) radius.err;
	    *((float *)op2) = (float) angle.val;
	    *((float *)op3) = (float) angle.err;

	}
	FUNC_MESS_END();
}
static PyUFuncGenericFunction PyGSL_sf_ufunc_qi_D_dd_one_data[] = {PyGSL_sf_ufunc_qi_D_dd_as_F_ff, PyGSL_sf_ufunc_qi_D_dd};


static char *PyGSL_polar_to_rect_doc = "Convert r, theta into a complex representation";
static char *PyGSL_rect_to_polar_doc = "Convert a complex into r and theta. Theta will be in the range [-pi,pi]";
static char * sf_module_doc = 
"GSL Special Functions wrapped mostly as Numeric ufuncs. \n\
Use print sf.<function>.__doc__ to get more information about the ufunc.\n\
Help does not display the information stored there!\n";

DL_EXPORT(void) initsf(void)
{
  PyObject* gsl_errors_module=NULL;
  PyObject* sf_module=NULL, *sf_dict=NULL;
  PyObject* f=NULL, *item=NULL;

  fprintf(stderr, "Module compiled at %s %s\n", __DATE__, __TIME__);
  sf_module=Py_InitModule("sf", sf_functions);
  module = sf_module;

  import_array();
  import_ufunc();
  init_pygsl();

  gsl_errors_module=PyImport_ImportModule ("pygsl.errors");
  if(gsl_errors_module == NULL)
       return;


  sf_dict=PyModule_GetDict(sf_module);
  if (!(item = PyString_FromString(sf_module_doc))){
       PyErr_SetString(PyExc_ImportError, 
		       "I could not generate module doc string!");
       goto fail;
  }
  if (PyDict_SetItemString(sf_dict, "__doc__", item) != 0){
       PyErr_SetString(PyExc_ImportError, 
		       "I could not add doc string to module dict!");
       goto fail;
  }

  PyDict_SetItemString(sf_dict, "PREC_DOUBLE", PyInt_FromLong((long) GSL_PREC_DOUBLE));
  PyDict_SetItemString(sf_dict, "PREC_SINGLE", PyInt_FromLong((long) GSL_PREC_SINGLE));
  PyDict_SetItemString(sf_dict, "PREC_APPROX", PyInt_FromLong((long) GSL_PREC_APPROX));

  DEBUG_MESS(1, "polar_to_rect %p", (void *)gsl_sf_polar_to_rect);
  f = PyUFunc_FromFuncAndData(PyGSL_sf_ufunc_qi_dd_D_one_data,
			      polar_to_rect_data,
			      PyGSL_sf_ufunc_qi_dd_D_one_types,
			      2, /* number of supported types */
			      2, /* in args */
			      2, /* out args */
			      PyUFunc_None,
			      "polar_to_rect",
			      PyGSL_polar_to_rect_doc,
			      0 /*check return*/);   
  PyDict_SetItemString(sf_dict, "polar_to_rect", f);
  /* Py_DECREF(f); */

  DEBUG_MESS(1, "rect_to_polar %p", (void *)gsl_sf_rect_to_polar);
  f = PyUFunc_FromFuncAndData(PyGSL_sf_ufunc_qi_D_dd_one_data,
			      rect_to_polar_data,
			      PyGSL_sf_ufunc_qi_D_dd_one_types,
			       2, /* number of supported types */
			       1, /* in args */
			       4, /* out args */
			       PyUFunc_None,
			       "rect_to_polar",
			       PyGSL_rect_to_polar_doc,
			       0 /*check return*/);   
  PyDict_SetItemString(sf_dict, "rect_to_polar", f);
  /* Py_DECREF(f); */

#ifdef IMPORTALL
#include "sf__objects.c"
#endif

  return;

 fail:
     if(!PyErr_Occurred()){
	  PyErr_SetString(PyExc_ImportError, "I could not init sf module!");
     }
}