## Example on wrapper code

Make a loop over a NumPy array in pure C:
```#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;
}
```