Index: /anuga_work/development/anuga_1d/comp_flux_ext.c
===================================================================
--- /anuga_work/development/anuga_1d/comp_flux_ext.c	(revision 5562)
+++ /anuga_work/development/anuga_1d/comp_flux_ext.c	(revision 5563)
@@ -5,13 +5,17 @@
 const double pi = 3.14159265358979;
 
-double max(double a, double b) {
-	double z;
-	z=(a>b)?a:b;
-	return z;}
-	
-double min(double a, double b) {
-	double z;
-	z=(a<b)?a:b;
-	return z;}
+// Shared code snippets
+#include "util_ext.h"
+
+
+/* double max(double a, double b) { */
+/* 	double z; */
+/* 	z=(a>b)?a:b; */
+/* 	return z;} */
+	
+/* double min(double a, double b) { */
+/* 	double z; */
+/* 	z=(a<b)?a:b; */
+/* 	return z;} */
 	
 
@@ -154,5 +158,5 @@
 				
 				normal = normals[ki];
-				_flux_function(ql, qr, zl, zr, normal, g, epsilon, flux, &max_speed);
+				_flux_function(ql, qr, zl, zr, normal, g, epsilon, edgeflux, &max_speed);
 				flux[0] -= edgeflux[0];
 				flux[1] -= edgeflux[1];
@@ -189,13 +193,15 @@
 //=========================================================================
 PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
-  PyArrayObject *neighbours, 
+  
+    PyArrayObject 
+	*neighbours, 
 	*neighbour_vertices,
-    *normals, 
+	*normals, 
 	*areas,
-    *stage_edge_values,
-    *xmom_edge_values,
-    *bed_edge_values,
-    *stage_boundary_values,
-    *xmom_boundary_values,
+	*stage_edge_values,
+	*xmom_edge_values,
+	*bed_edge_values,
+	*stage_boundary_values,
+	*xmom_boundary_values,
 	*stage_explicit_update,
 	*xmom_explicit_update,
@@ -231,23 +237,130 @@
   // the explicit update arrays 
   timestep = _compute_fluxes_ext(timestep,
-				     epsilon,
-				     g,
-				     (long*) neighbours -> data,
-				     (long*) neighbour_vertices -> data,
-				     (double*) normals -> data,
-				     (double*) areas -> data,
-				     (double*) stage_edge_values -> data,
-				     (double*) xmom_edge_values -> data,
-				     (double*) bed_edge_values -> data,
-				     (double*) stage_boundary_values -> data,
-				     (double*) xmom_boundary_values -> data,
-					 (double*) stage_explicit_update -> data,
-					 (double*) xmom_explicit_update -> data,
-					 number_of_elements,
-					 (double*) max_speed_array -> data);
+				 epsilon,
+				 g,
+				 (long*) neighbours -> data,
+				 (long*) neighbour_vertices -> data,
+				 (double*) normals -> data,
+				 (double*) areas -> data,
+				 (double*) stage_edge_values -> data,
+				 (double*) xmom_edge_values -> data,
+				 (double*) bed_edge_values -> data,
+				 (double*) stage_boundary_values -> data,
+				 (double*) xmom_boundary_values -> data,
+				 (double*) stage_explicit_update -> data,
+				 (double*) xmom_explicit_update -> data,
+				 number_of_elements,
+				 (double*) max_speed_array -> data);
+
   // Return updated flux timestep
   return Py_BuildValue("d", timestep);
 }
 
+
+PyObject *compute_fluxes_ext_short(PyObject *self, PyObject *args) {
+  
+   PyObject 
+	*domain,
+	*stage, 
+	*xmom, 
+	*bed;
+
+    PyArrayObject 
+	*neighbours, 
+	*neighbour_vertices,
+	*normals, 
+	*areas,
+	*stage_vertex_values,
+	*xmom_vertex_values,
+	*bed_vertex_values,
+	*stage_boundary_values,
+	*xmom_boundary_values,
+	*stage_explicit_update,
+	*xmom_explicit_update,
+	*max_speed_array;
+    
+  double timestep, epsilon, g;
+  int number_of_elements;
+
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dOOOO",
+			&timestep,
+			&domain,
+			&stage,
+			&xmom,
+			&bed)) {
+      PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext_short could not parse input");
+      return NULL;
+  }
+
+
+    epsilon           = get_python_double(domain,"epsilon");
+    g                 = get_python_double(domain,"g");
+ 
+    
+    neighbours        = get_consecutive_array(domain, "neighbours");
+    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 
+    normals           = get_consecutive_array(domain, "normals");
+    areas             = get_consecutive_array(domain, "areas");    
+    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
+    
+    stage_vertex_values = get_consecutive_array(stage, "vertex_values");    
+    xmom_vertex_values  = get_consecutive_array(xmom, "vertex_values");    
+    bed_vertex_values   = get_consecutive_array(bed, "vertex_values");    
+
+    stage_boundary_values = get_consecutive_array(stage, "boundary_values");    
+    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");    
+
+
+    stage_explicit_update = get_consecutive_array(stage, "explicit_update");    
+    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");    
+
+
+
+    number_of_elements = stage_vertex_values -> dimensions[0];
+
+
+ 
+    // Call underlying flux computation routine and update 
+    // the explicit update arrays 
+    timestep = _compute_fluxes_ext(timestep,
+				 epsilon,
+				 g,
+				 (long*) neighbours -> data,
+				 (long*) neighbour_vertices -> data,
+				 (double*) normals -> data,
+				 (double*) areas -> data,
+				 (double*) stage_vertex_values -> data,
+				 (double*) xmom_vertex_values -> data,
+				 (double*) bed_vertex_values -> data,
+				 (double*) stage_boundary_values -> data,
+				 (double*) xmom_boundary_values -> data,
+				 (double*) stage_explicit_update -> data,
+				 (double*) xmom_explicit_update -> data,
+				 number_of_elements,
+				 (double*) max_speed_array -> data);
+
+
+  Py_DECREF(neighbours);
+  Py_DECREF(neighbour_vertices);
+  Py_DECREF(normals);
+  Py_DECREF(areas);
+  Py_DECREF(stage_vertex_values);
+  Py_DECREF(xmom_vertex_values);
+  Py_DECREF(bed_vertex_values);
+  Py_DECREF(stage_boundary_values);
+  Py_DECREF(xmom_boundary_values);
+  Py_DECREF(stage_explicit_update);
+  Py_DECREF(xmom_explicit_update);
+  Py_DECREF(max_speed_array);
+
+
+
+
+  // Return updated flux timestep
+  return Py_BuildValue("d", timestep);
+}
+
  
 
@@ -259,4 +372,5 @@
 static struct PyMethodDef MethodTable[] = {
   {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
+  {"compute_fluxes_ext_short", compute_fluxes_ext_short, METH_VARARGS, "Print out"},
   {NULL, NULL}
 };
Index: /anuga_work/development/anuga_1d/domain.py
===================================================================
--- /anuga_work/development/anuga_1d/domain.py	(revision 5562)
+++ /anuga_work/development/anuga_1d/domain.py	(revision 5563)
@@ -47,4 +47,6 @@
         self.centroids = zeros(N, Float)
         self.areas     = zeros(N, Float)
+
+        self.max_speed_array = zeros(N, Float)
 
         self.normals = zeros((N, 2), Float)
Index: /anuga_work/development/anuga_1d/shallow_water_domain.py
===================================================================
--- /anuga_work/development/anuga_1d/shallow_water_domain.py	(revision 5562)
+++ /anuga_work/development/anuga_1d/shallow_water_domain.py	(revision 5563)
@@ -90,4 +90,6 @@
 
         self.quantities_to_be_stored = ['stage','xmomentum']
+
+
 
 
@@ -155,5 +157,5 @@
         #Call correct module function
         #(either from this module or C-extension)
-        compute_fluxes(self)
+        compute_fluxes_C_short(self)
 
     def compute_timestep(self):
@@ -572,5 +574,121 @@
     #print domain.quantities['stage'].centroid_values
 
+
+# Compute flux definition
+def compute_fluxes_C_long(domain):
+	from Numeric import zeros, Float
+	import sys
+
+	
+	timestep = float(sys.maxint)
+	#print 'timestep=',timestep
+	#print 'The type of timestep is',type(timestep)
+	
+	epsilon = domain.epsilon
+	#print 'epsilon=',epsilon
+	#print 'The type of epsilon is',type(epsilon)
+	
+	g = domain.g
+	#print 'g=',g
+	#print 'The type of g is',type(g)
+	
+	neighbours = domain.neighbours
+	#print 'neighbours=',neighbours
+	#print 'The type of neighbours is',type(neighbours)
+	
+	neighbour_vertices = domain.neighbour_vertices
+	#print 'neighbour_vertices=',neighbour_vertices
+	#print 'The type of neighbour_vertices is',type(neighbour_vertices)
+	
+	normals = domain.normals
+	#print 'normals=',normals
+	#print 'The type of normals is',type(normals)
+	
+	areas = domain.areas
+	#print 'areas=',areas
+	#print 'The type of areas is',type(areas)
+	
+	stage_edge_values = domain.quantities['stage'].vertex_values
+	#print 'stage_edge_values=',stage_edge_values
+	#print 'The type of stage_edge_values is',type(stage_edge_values)
+	
+	xmom_edge_values = domain.quantities['xmomentum'].vertex_values
+	#print 'xmom_edge_values=',xmom_edge_values
+	#print 'The type of xmom_edge_values is',type(xmom_edge_values)
+	
+	bed_edge_values = domain.quantities['elevation'].vertex_values
+	#print 'bed_edge_values=',bed_edge_values
+	#print 'The type of bed_edge_values is',type(bed_edge_values)
+	
+	stage_boundary_values = domain.quantities['stage'].boundary_values
+	#print 'stage_boundary_values=',stage_boundary_values
+	#print 'The type of stage_boundary_values is',type(stage_boundary_values)
+	
+	xmom_boundary_values = domain.quantities['xmomentum'].boundary_values
+	#print 'xmom_boundary_values=',xmom_boundary_values
+	#print 'The type of xmom_boundary_values is',type(xmom_boundary_values)
+	
+	stage_explicit_update = domain.quantities['stage'].explicit_update
+	#print 'stage_explicit_update=',stage_explicit_update
+	#print 'The type of stage_explicit_update is',type(stage_explicit_update)
+	
+	xmom_explicit_update = domain.quantities['xmomentum'].explicit_update
+	#print 'xmom_explicit_update=',xmom_explicit_update
+	#print 'The type of xmom_explicit_update is',type(xmom_explicit_update)
+	
+	number_of_elements = len(stage_edge_values)
+	#print 'number_of_elements=',number_of_elements
+	#print 'The type of number_of_elements is',type(number_of_elements)
+	
+	max_speed_array = domain.max_speed_array
+	#print 'max_speed_array=',max_speed_array
+	#print 'The type of max_speed_array is',type(max_speed_array)
+        
+        
+        from comp_flux_ext import compute_fluxes_ext
+	
+	domain.timestep = compute_fluxes_ext(timestep,
+                                  epsilon,
+                                  g,
+                                  neighbours,
+                                  neighbour_vertices,
+                                  normals,
+                                  areas, 
+                                  stage_edge_values,
+                                  xmom_edge_values,
+                                  bed_edge_values,
+                                  stage_boundary_values,
+                                  xmom_boundary_values,
+                                  stage_explicit_update,
+                                  xmom_explicit_update,
+                                  number_of_elements,
+                                  max_speed_array)
+
+
+# Compute flux definition
+def compute_fluxes_C_short(domain):
+	from Numeric import zeros, Float
+	import sys
+
+	
+	timestep = float(sys.maxint)
+
+        stage = domain.quantities['stage']
+        xmom = domain.quantities['xmomentum']
+        bed = domain.quantities['elevation']
+
+
+        from comp_flux_ext import compute_fluxes_ext_short
+
+ 	domain.timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
+
+
+
 ####################################
+
+
+
+
+
 
 # Module functions for gradient limiting (distribute_to_vertices_and_edges)
Index: /anuga_work/development/anuga_1d/shallow_water_try.py
===================================================================
--- /anuga_work/development/anuga_1d/shallow_water_try.py	(revision 5562)
+++ /anuga_work/development/anuga_1d/shallow_water_try.py	(revision 5563)
Index: /anuga_work/development/anuga_1d/test_shallow_water.py
===================================================================
--- /anuga_work/development/anuga_1d/test_shallow_water.py	(revision 5562)
+++ /anuga_work/development/anuga_1d/test_shallow_water.py	(revision 5563)
@@ -41,4 +41,49 @@
         assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
         assert allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
+
+
+    def test_local_flux_function(self):
+        normal = 1.0
+        ql = array([1.0, 2.0],Float)
+        qr = array([1.0, 2.0],Float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
+
+        assert allclose(array([2.0, 8.9],Float), edgeflux)
+        assert allclose(5.1304951685, maxspeed)
+
+        normal = -1.0
+        ql = array([1.0, 2.0],Float)
+        qr = array([1.0, 2.0],Float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
+
+        assert allclose(array([-2.0, -8.9],Float), edgeflux)
+        assert allclose(5.1304951685, maxspeed)
+
+    def test_domain_flux_function(self):
+        normal = 1.0
+        ql = array([1.0, 2.0],Float)
+        qr = array([1.0, 2.0],Float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
+
+        #print edgeflux
+
+        from shallow_water_domain import flux_function as domain_flux_function
+
+        domainedgeflux, domainmaxspeed = domain_flux_function(normal, ql,qr,zl,zr)
+
+        #print domainedgeflux
+
+        assert allclose(domainedgeflux, edgeflux)
+        assert allclose(domainmaxspeed, maxspeed)
+
 
 
@@ -128,4 +173,5 @@
                 #m = domain.neighbour_edges[k,i]
                 m = domain.neighbour_vertices[k,i]
+                #print i, ' ' , m
                 #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
                 qr[0] = stage[n, m]
@@ -169,4 +215,101 @@
 
 
+def flux_function(normal, ql, qr, zl, zr):
+    """Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of w = h+z using the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+
+    Conserved quantities w, uh, are stored as elements 0 and 1
+    in the numerical vectors ql an qr.
+
+    Bed elevations zl and zr.
+    """
+
+    from config import g, epsilon
+    from math import sqrt
+    from Numeric import array
+
+    #print 'ql',ql
+
+    #Align momentums with x-axis
+    #q_left  = rotate(ql, normal, direction = 1)
+    #q_right = rotate(qr, normal, direction = 1)
+    q_left = ql
+    q_left[1] = q_left[1]*normal
+    q_right = qr
+    q_right[1] = q_right[1]*normal
+
+    #z = (zl+zr)/2 #Take average of field values
+    z = 0.5*(zl+zr) #Take average of field values
+
+    w_left  = q_left[0]   #w=h+z
+    h_left  = w_left-z
+    uh_left = q_left[1]
+
+    if h_left < epsilon:
+        u_left = 0.0  #Could have been negative
+        h_left = 0.0
+    else:
+        u_left  = uh_left/h_left
+
+
+    w_right  = q_right[0]  #w=h+z
+    h_right  = w_right-z
+    uh_right = q_right[1]
+
+
+    if h_right < epsilon:
+        u_right = 0.0  #Could have been negative
+        h_right = 0.0
+    else:
+        u_right  = uh_right/h_right
+
+    #vh_left  = q_left[2]
+    #vh_right = q_right[2]
+
+    #print h_right
+    #print u_right
+    #print h_left
+    #print u_right
+
+    soundspeed_left  = sqrt(g*h_left)
+    soundspeed_right = sqrt(g*h_right)
+
+    #Maximal wave speed
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
+    
+    #Minimal wave speed
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
+    
+    #Flux computation
+
+    #flux_left  = array([u_left*h_left,
+    #                    u_left*uh_left + 0.5*g*h_left**2])
+    #flux_right = array([u_right*h_right,
+    #                    u_right*uh_right + 0.5*g*h_right**2])
+    flux_left  = array([u_left*h_left,
+                        u_left*uh_left + 0.5*g*h_left*h_left])
+    flux_right = array([u_right*h_right,
+                        u_right*uh_right + 0.5*g*h_right*h_right])
+
+    denom = s_max-s_min
+    if denom == 0.0:
+        edgeflux = array([0.0, 0.0])
+        max_speed = 0.0
+    else:
+        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
+        edgeflux += s_max*s_min*(q_right-q_left)/denom
+        
+        edgeflux[1] = edgeflux[1]*normal
+
+        max_speed = max(abs(s_max), abs(s_min))
+
+    return edgeflux, max_speed
+
 
 #-------------------------------------------------------------
Index: /anuga_work/development/anuga_1d/util_ext.h
===================================================================
--- /anuga_work/development/anuga_1d/util_ext.h	(revision 5563)
+++ /anuga_work/development/anuga_1d/util_ext.h	(revision 5563)
@@ -0,0 +1,334 @@
+// Python - C extension for finite_volumes util module.
+//
+// To compile (Python2.3):
+//  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O 
+//  gcc -shared util_ext.o  -o util_ext.so	
+//
+// See the module util.py
+//
+//
+// Ole Nielsen, GA 2004
+	
+#include "Python.h"	
+#include "Numeric/arrayobject.h"
+#include "math.h"
+
+
+double max(double x, double y) {  
+  //Return maximum of two doubles
+  
+  if (x > y) return x;
+  else return y;
+}
+
+
+double min(double x, double y) {  
+  //Return minimum of two doubles
+  
+  if (x < y) return x;
+  else return y;
+}
+
+
+int _gradient(double x0, double y0, 
+	      double x1, double y1, 
+	      double x2, double y2, 
+	      double q0, double q1, double q2, 
+	      double *a, double *b) {
+	      
+  /*Compute gradient (a,b) based on three points (x0,y0), (x1,y1) and (x2,y2) 
+  with values q0, q1 and q2.
+  
+  Extrapolation formula (q0 is selected as an arbitrary origin)
+    q(x,y) = q0 + a*(x-x0) + b*(y-y0)                    (1)
+  
+  Substituting the known values for q1 and q2 into (1) yield the 
+  equations for a and b 
+  
+      q1-q0 = a*(x1-x0) + b*(y1-y0)                      (2)
+      q2-q0 = a*(x2-x0) + b*(y2-y0)                      (3)      
+      
+  or in matrix form
+  
+  /               \  /   \   /       \  
+  |  x1-x0  y1-y0 |  | a |   | q1-q0 |
+  |               |  |   | = |       | 
+  |  x2-x0  y2-y0 |  | b |   | q2-q0 |
+  \               /  \   /   \       /
+   
+  which is solved using the standard determinant technique    
+      
+  */
+	      
+
+  double det;
+  
+  det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0);
+
+  *a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0);
+  *a /= det;
+
+  *b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0);
+  *b /= det;
+
+  return 0;
+}
+
+
+int _gradient2(double x0, double y0, 
+	       double x1, double y1, 
+	       double q0, double q1, 
+	       double *a, double *b) {
+  /*Compute gradient (a,b) between two points (x0,y0) and (x1,y1) 
+  with values q0 and q1 such that the plane is constant in the direction 
+  orthogonal to (x1-x0, y1-y0).
+  
+  Extrapolation formula
+    q(x,y) = q0 + a*(x-x0) + b*(y-y0)                    (1)
+  
+  Substituting the known values for q1 into (1) yields an 
+  under determined  equation for a and b 
+      q1-q0 = a*(x1-x0) + b*(y1-y0)                      (2)
+      
+      
+  Now add the additional requirement that the gradient in the direction 
+  orthogonal to (x1-x0, y1-y0) should be zero. The orthogonal direction 
+  is given by the vector (y0-y1, x1-x0).
+  
+  Define the point (x2, y2) = (x0 + y0-y1, y0 + x1-x0) on the orthognal line. 
+  Then we know that the corresponding value q2 should be equal to q0 in order 
+  to obtain the zero gradient, hence applying (1) again    
+      q0 = q2 = q(x2, y2) = q0 + a*(x2-x0) + b*(y2-y0)
+                          = q0 + a*(x0 + y0-y1-x0) + b*(y0 + x1-x0 - y0)
+			  = q0 + a*(y0-y1) + b*(x1-x0)
+			  
+  leads to the orthogonality constraint
+     a*(y0-y1) + b*(x1-x0) = 0                           (3) 
+     
+  which closes the system and yields
+  
+  /               \  /   \   /       \  
+  |  x1-x0  y1-y0 |  | a |   | q1-q0 |
+  |               |  |   | = |       | 
+  |  y0-y1  x1-x0 |  | b |   |   0   |
+  \               /  \   /   \       /
+   
+  which is solved using the standard determinant technique    
+      
+  */
+
+  double det, xx, yy, qq;
+  
+  xx = x1-x0;
+  yy = y1-y0;
+  qq = q1-q0;
+    
+  det = xx*xx + yy*yy;  //FIXME  catch det == 0
+  *a = xx*qq/det;
+  *b = yy*qq/det;
+        
+  return 0;
+}
+
+
+void _limit_old(int N, double beta, double* qc, double* qv, 
+	    double* qmin, double* qmax) { 
+
+  //N are the number of elements
+  int k, i, k3;
+  double dq, dqa[3], phi, r;
+  
+  //printf("INSIDE\n");
+  for (k=0; k<N; k++) {
+    k3 = k*3;
+    
+    //Find the gradient limiter (phi) across vertices  
+    phi = 1.0;
+    for (i=0; i<3; i++) {    
+      r = 1.0;
+      
+      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
+      dqa[i] = dq;              //Save dq for use in the next loop
+      
+      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
+      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;      
+  
+  
+      phi = min( min(r*beta, 1.0), phi);    
+    }
+    
+    //Then update using phi limiter
+    for (i=0; i<3; i++) {    
+      qv[k3+i] = qc[k] + phi*dqa[i];
+    }
+  }
+}
+
+
+void  print_double_array(char* name, double* array, int n, int m){
+
+    int k,i,km;
+
+    printf("%s = [",name);
+    for (k=0; k<n; k++){
+	km = m*k;
+	printf("[");
+	for (i=0; i<m ; i++){
+	    printf("%g ",array[km+i]);
+	}
+	if (k==(n-1))
+	    printf("]");
+	else
+	    printf("]\n");
+    }
+    printf("]\n");
+}
+
+void  print_int_array(char* name, int* array, int n, int m){
+
+    int k,i,km;
+
+    printf("%s = [",name);
+    for (k=0; k<n; k++){
+	km = m*k;
+	printf("[");
+	for (i=0; i<m ; i++){
+	    printf("%i ",array[km+i]);
+	}
+	if (k==(n-1))
+	    printf("]");
+	else
+	    printf("]\n");
+    }
+    printf("]\n");
+}
+
+
+void  print_long_array(char* name, long* array, int n, int m){
+
+    int k,i,km;
+
+    printf("%s = [",name);
+    for (k=0; k<n; k++){
+	km = m*k;
+	printf("[");
+	for (i=0; i<m ; i++){
+	  printf("%i ",(int) array[km+i]);
+	}
+	if (k==(n-1))
+	    printf("]");
+	else
+	    printf("]\n");
+    }
+    printf("]\n");
+}
+
+void print_numeric_array(PyArrayObject *x) {  
+  int i, j;
+  for (i=0; i<x->dimensions[0]; i++) {  
+    for (j=0; j<x->dimensions[1]; j++) {
+      printf("%f ", *(double*) (x->data + i*x->strides[0] + j*x->strides[1]));
+    }
+    printf("\n");  
+  }
+  printf("\n");    
+}
+
+void print_numeric_vector(PyArrayObject *x) {  
+  int i;
+  for (i=0; i<x->dimensions[0]; i++) {
+    printf("%f ", *(double*) (x->data + i*x->strides[0]));  
+  }
+  printf("\n");  
+}
+
+PyArrayObject *get_consecutive_array(PyObject *O, char *name) {
+  PyArrayObject *A, *B;
+  
+
+  //Get array object from attribute
+  
+  /*
+  //FIXME: THE TEST DOESN't WORK
+  printf("Err = %d\n", PyObject_HasAttrString(O, name));
+  if (PyObject_HasAttrString(O, name) == 1) {
+    B = (PyArrayObject*) PyObject_GetAttrString(O, name);
+    if (!B) return NULL;
+  } else {
+    return NULL;
+    } 
+  */
+    
+  B = (PyArrayObject*) PyObject_GetAttrString(O, name);
+  if (!B) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
+    return NULL;
+  }     
+  
+  //Convert to consecutive array
+  A = (PyArrayObject*) PyArray_ContiguousFromObject((PyObject*) B, 
+						    B -> descr -> type, 0, 0);
+  
+  Py_DECREF(B); //FIXME: Is this really needed??
+  
+  if (!A) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
+    return NULL;
+  }
+  return A;
+}
+
+double get_python_double(PyObject *O, char *name) {
+  PyObject *TObject;
+  double tmp;
+  
+
+  //Get double from attribute
+  TObject = PyObject_GetAttrString(O, name);
+  if (!TObject) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_double could not obtain double from object");
+    return 0.0;
+  }  
+  
+  tmp = PyFloat_AsDouble(TObject);
+  
+  Py_DECREF(TObject);
+  
+  return tmp;
+}
+
+int get_python_integer(PyObject *O, char *name) {
+  PyObject *TObject;
+  int tmp;
+  
+
+  //Get double from attribute
+  TObject = PyObject_GetAttrString(O, name);
+  if (!TObject) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_integer could not obtain double from object");
+    return 0;
+  }  
+  
+  tmp = PyFloat_AsDouble(TObject);
+  
+  Py_DECREF(TObject);
+  
+  return tmp;
+}
+
+
+PyObject *get_python_object(PyObject *O, char *name) {
+  PyObject *Oout;
+
+  Oout = PyObject_GetAttrString(O, name);
+  if (!Oout) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_object could not obtain object");
+    return NULL;
+  }
+
+  return Oout;
+
+}
+
+
+
