The loop is straightforward

  int i,j;
  for (i = 0; i < nx; i++) {
    for (j = 0; j < ny; j++) {
      arglist = Py_BuildValue("(dd)", xcoor(i), ycoor(j));
      result = PyEval_CallObject(func1, arglist);
      a(i,j) = PyFloat_AS_DOUBLE(result);
    }
  }

  return PyArray_Return(a.getPtr());