Example on wrapper code

#include <Python.h>              /*Python seen from C*/
#include <Numeric/arrayobject.h> /*NumPy  seen from C*/
#include <math.h>

static PyObject *fillarr_Numeric_C(PyObject *self, PyObject* args)
{
  PyArrayObject *array; /* C representation of Numeric array */
  int i,n;              /* number of array entries (x points) */
  double x;
  double* a;            /* C ptr to the NumPy array */

  /* parse the arguments to this function using Python tools */
  if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &array)) {
    /* extracting a NumPy array as argument was not successful */
    return NULL;  /* Error indicator */
  }

  /* check that we have a one-dimensional array */
  if (array->nd != 1) {
    /* throw a Python exception (ValueError) */
    PyErr_SetString(PyExc_ValueError,
                    "the NumPy array must be one-dimensional");
    return NULL;
  }
  /* check that the datatype is NumPy/Python float, i.e. C double */
  if (array->descr->type_num != PyArray_DOUBLE) {
    PyErr_SetString(PyExc_ValueError,
                    "the NumPy array must consist of floats");
    return NULL;
  }

  n = array->dimensions[0];    /* length of the array */
  a = (double*) array->data;   /* pointer to the NumPy data */

  /* what we came here to do: fill a with f(x) values */
  for (i=0; i<n; i++) {
    x = (double) i/(n-1);
    a[i] = f(x);
    /* more general: *(array->data+i*array->strides[0]) = f(x) */
  }

  /* return nothing; we just borrow a reference to the Python
     object */
  Py_INCREF(Py_None);  return Py_None;
}   
/* 
   the method table must always be present - it lists the 
   functions that should be callable from Python: 
*/
static PyMethodDef FillarrMethods[] = {
  {"fillarr_Numeric_C", /* name of func when called from Python */
   fillarr_Numeric_C,   /* corresponding C function */
   METH_VARARGS},
  {NULL, NULL}
};

void initfillarr_basic_C()
{
  /* 
     assign the name of the module and the name of the
     method table:
  */
  (void) Py_InitModule("fillarr_basic_C", FillarrMethods);
  import_array();   /* NumPy initialization */
}

double f (double x)
{
  double n, c, cH, A, c1, u;
  n = 0.3; c = 1.0 + 1/n; cH = pow(0.5,c);
  A = n/(1.0 + n);
  if (x <= 0.5) c1 = pow(0.5 - x,c);
  else          c1 = pow(x - 0.5,c);
  u = A*(cH-c1);
  return u;
}