Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 4989)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 4990)
@@ -1404,4 +1404,9 @@
         extrapolate_from_gradient(self)
 
+    def extrapolate_second_order_and_limit(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        extrapolate_second_order_and_limit(self)        
+
     def backup_centroid_values(self):
         #Call correct module function
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 4989)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 4990)
@@ -16,7 +16,9 @@
 //Shared code snippets
 #include "util_ext.h"
-//#include "quantity_ext.h" //Obsolete
-
-
+
+
+//-------------------------------------------
+// Low level routines (called from wrappers)
+//------------------------------------------
 
 int _compute_gradients(int N,
@@ -90,5 +92,4 @@
 	return 0;
 }
-
 
 
@@ -135,4 +136,98 @@
 	return 0;
 }
+
+
+int _extrapolate_and_limit_from_gradient(int N,double beta,
+					 double* centroids,
+					 long*   neighbours,
+					 double* centroid_values,
+					 double* vertex_coordinates,
+					 double* vertex_values,
+					 double* edge_values,
+					 double* x_gradient,
+					 double* y_gradient) {
+
+	int i, k, k2, k3, k6;
+	double x, y, x0, y0, x1, y1, x2, y2;
+	long n;
+	double qmin, qmax, qn, qc;
+	double dq, dqa[3], phi, r;
+
+
+	for (k=0; k<N; k++){
+		k6 = 6*k;
+		k3 = 3*k;
+		k2 = 2*k;
+
+		// Centroid coordinates
+		x = centroids[k2]; y = centroids[k2+1];
+
+		// vertex coordinates
+		// x0, y0, x1, y1, x2, y2 = X[k,:]
+		x0 = vertex_coordinates[k6 + 0];
+		y0 = vertex_coordinates[k6 + 1];
+		x1 = vertex_coordinates[k6 + 2];
+		y1 = vertex_coordinates[k6 + 3];
+		x2 = vertex_coordinates[k6 + 4];
+		y2 = vertex_coordinates[k6 + 5];
+
+		//Centroid Value
+		qc = centroid_values[k];
+
+		// Extrapolate to Vertices
+		vertex_values[k3+0] = centroid_values[k] + x_gradient[k]*(x0-x) + y_gradient[k]*(y0-y);
+		vertex_values[k3+1] = centroid_values[k] + x_gradient[k]*(x1-x) + y_gradient[k]*(y1-y);
+		vertex_values[k3+2] = centroid_values[k] + x_gradient[k]*(x2-x) + y_gradient[k]*(y2-y);
+
+		// Extrapolate to Edges (midpoints)
+		edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
+		edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
+		edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
+
+
+		phi = 1.0;
+
+		for (i=0; i<3; i++) {
+		    dq = edge_values[k3+i] - qc;      //Delta between edge and centroid values
+		    dqa[i] = dq;                      //Save dq for use in updating vertex values
+
+		    n = neighbours[k3+i];
+		    if (n >= 0) {
+			qn = centroid_values[n];     //Neighbour's centroid value
+
+			qmin = min(qc, qn);
+			qmax = max(qc, qn);
+
+
+			r = 1.0;
+      
+			if (dq > 0.0) r = (qmax - qc)/dq;
+			if (dq < 0.0) r = (qmin - qc)/dq;      
+		    
+			phi = min( min(r*beta, 1.0), phi);    
+		    }
+		}
+
+
+		//Update gradient, edge and vertex values using phi limiter
+		x_gradient[k] = x_gradient[k]*phi;
+		y_gradient[k] = y_gradient[k]*phi;
+
+		edge_values[k3+0] = qc + phi*dqa[0];
+		edge_values[k3+1] = qc + phi*dqa[1];
+		edge_values[k3+2] = qc + phi*dqa[2];
+
+
+		vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
+		vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
+		vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
+
+	}
+
+	return 0;
+
+}
+
+
 
 
@@ -542,6 +637,8 @@
 
 
-/////////////////////////////////////////////////
-// Gateways to Python
+//-----------------------------------------------------
+// Python method Wrappers 
+//-----------------------------------------------------
+
 PyObject *update(PyObject *self, PyObject *args) {
   // FIXME (Ole): It would be great to turn this text into a Python DOC string
@@ -937,30 +1034,4 @@
 	Py_DECREF(domain);
 
-	//Allocate space for return vectors a and b (don't DECREF)
-	//dimensions[0] = N;
-	//a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
-	//b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
-
-	//FIXME: Odd that I couldn't use normal arrays
-	//Allocate space for return vectors a and b (don't DECREF)
-	//a = (double*) malloc(N * sizeof(double));
-	//if (!a) return NULL;
-	//b = (double*) malloc(N * sizeof(double));
-	//if (!b) return NULL;
-
-
-/* 	err = _compute_gradients(N, */
-/* 			(double*) centroids -> data, */
-/* 			(double*) centroid_values -> data, */
-/* 			(long*) number_of_boundaries -> data, */
-/* 			(long*) surrogate_neighbours -> data, */
-/* 			(double*) x_gradient -> data, */
-/* 			(double*) y_gradient -> data); */
-
-/* 	if (err != 0) { */
-/* 	  PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); */
-/* 	  return NULL; */
-/* 	} */
-
 	err = _extrapolate_from_gradient(N,
 			(double*) centroids -> data,
@@ -993,4 +1064,109 @@
 
 	return Py_BuildValue("");
+}
+
+
+PyObject *extrapolate_second_order_and_limit(PyObject *self, PyObject *args) {
+    /* Compute edge values using second order approximation and limit values 
+       so that edge values are limited by the two corresponding centroid values
+       
+       Python Call:
+       extrapolate_second_order_and_limit(domain,quantity,beta)
+    */
+
+    PyObject *quantity, *domain;
+	
+    PyArrayObject
+	*domain_centroids,              //Coordinates at centroids
+	*domain_vertex_coordinates,     //Coordinates at vertices
+	*domain_number_of_boundaries,   //Number of boundaries for each triangle
+	*domain_surrogate_neighbours,   //True neighbours or - if one missing - self
+	*domain_neighbours,             //True neighbours, or if negative a link to boundary
+
+	*quantity_centroid_values,      //Values at centroids	
+	*quantity_vertex_values,        //Values at vertices
+	*quantity_edge_values,          //Values at edges
+	*quantity_x_gradient,           //x gradient
+	*quantity_y_gradient;           //y gradient
+    
+
+  // Local variables
+  int ntri;
+  double beta;
+  int err;
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOd", 
+			&domain,
+			&quantity,
+			&beta)) {
+      PyErr_SetString(PyExc_RuntimeError, 
+		      "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");	
+      return NULL;
+  }
+
+	
+  // Get pertinent variables
+  domain_centroids              = get_consecutive_array(domain,   "centroid_coordinates");
+  domain_surrogate_neighbours   = get_consecutive_array(domain,   "surrogate_neighbours");
+  domain_number_of_boundaries   = get_consecutive_array(domain,   "number_of_boundaries");
+  domain_vertex_coordinates     = get_consecutive_array(domain,   "vertex_coordinates");
+  domain_neighbours             = get_consecutive_array(domain,   "neighbours");
+
+  quantity_centroid_values      = get_consecutive_array(quantity, "centroid_values");
+  quantity_vertex_values        = get_consecutive_array(quantity, "vertex_values");
+  quantity_edge_values          = get_consecutive_array(quantity, "edge_values");
+  quantity_x_gradient           = get_consecutive_array(quantity, "x_gradient");
+  quantity_y_gradient           = get_consecutive_array(quantity, "y_gradient");
+
+  ntri = quantity_centroid_values -> dimensions[0];
+
+
+
+  err = _compute_gradients(ntri,
+			   (double*) domain_centroids -> data,
+			   (double*) quantity_centroid_values -> data,
+			   (long*)   domain_number_of_boundaries -> data,
+			   (long*)   domain_surrogate_neighbours -> data,
+			   (double*) quantity_x_gradient -> data,
+			   (double*) quantity_y_gradient -> data);
+
+  if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+			  "quantity_ext.c: Internal function _compute_gradient failed");
+      return NULL;
+  }
+
+
+  err = _extrapolate_and_limit_from_gradient(ntri, beta,
+					     (double*) domain_centroids -> data,
+					     (long*)   domain_neighbours -> data,
+					     (double*) quantity_centroid_values -> data,
+					     (double*) domain_vertex_coordinates -> data,
+					     (double*) quantity_vertex_values -> data,
+					     (double*) quantity_edge_values -> data,
+					     (double*) quantity_x_gradient -> data,
+					     (double*) quantity_y_gradient -> data);
+
+  if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+			  "quantity_ext.c: Internal function _extrapolate_and_limit_from_gradient failed");
+      return NULL;
+  }
+
+
+  // Release
+  Py_DECREF(domain_centroids);
+  Py_DECREF(domain_surrogate_neighbours);
+  Py_DECREF(domain_number_of_boundaries);
+  Py_DECREF(domain_vertex_coordinates);
+
+  Py_DECREF(quantity_centroid_values);
+  Py_DECREF(quantity_vertex_values);
+  Py_DECREF(quantity_edge_values);
+  Py_DECREF(quantity_x_gradient);
+  Py_DECREF(quantity_y_gradient);
+
+  return Py_BuildValue("");
 }
 
@@ -1046,16 +1222,4 @@
 	Py_DECREF(domain);
 
-	//Allocate space for return vectors a and b (don't DECREF)
-	//dimensions[0] = N;
-	//a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
-	//b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
-
-	//FIXME: Odd that I couldn't use normal arrays
-	//Allocate space for return vectors a and b (don't DECREF)
-	//a = (double*) malloc(N * sizeof(double));
-	//if (!a) return NULL;
-	//b = (double*) malloc(N * sizeof(double));
-	//if (!b) return NULL;
-
 
 	err = _compute_gradients(N,
@@ -1071,20 +1235,4 @@
 	  return NULL;
 	}
-
-/* 	err = _extrapolate(N, */
-/* 			(double*) centroids -> data, */
-/* 			(double*) centroid_values -> data, */
-/* 			(double*) vertex_coordinates -> data, */
-/* 			(double*) vertex_values -> data, */
-/* 			(double*) edge_values -> data, */
-/* 			(double*) x_gradient -> data, */
-/* 			(double*) y_gradient -> data); */
-
-
-/* 	if (err != 0) { */
-/* 	  PyErr_SetString(PyExc_RuntimeError, */
-/* 			  "Internal function _extrapolate failed"); */
-/* 	  return NULL; */
-/* 	} */
 
 
@@ -1558,4 +1706,6 @@
 	{"extrapolate_from_gradient", extrapolate_from_gradient,
 		METH_VARARGS, "Print out"},
+	{"extrapolate_second_order_and_limit", extrapolate_second_order_and_limit,
+		METH_VARARGS, "Print out"},
 	{"interpolate_from_vertices_to_edges",
 		interpolate_from_vertices_to_edges,
