Index: /inundation/ga/storm_surge/pyvolution/config.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/config.py	(revision 239)
+++ /inundation/ga/storm_surge/pyvolution/config.py	(revision 240)
@@ -61,6 +61,6 @@
 #use_extensions = False   #Do not use C-extensions
 
-#use_psyco = True  #Use psyco optimisations
-use_psyco = False  #Do not use psyco optimisations
+use_psyco = True  #Use psyco optimisations
+#use_psyco = False  #Do not use psyco optimisations
 
 
Index: /inundation/ga/storm_surge/pyvolution/doit.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/doit.py	(revision 239)
+++ /inundation/ga/storm_surge/pyvolution/doit.py	(revision 240)
@@ -1,12 +1,6 @@
-
-
 import os
 
-
-
-os.system('python compile_all.py')
-
+os.system('python compile.py')
 os.system('python test_all.py')
-
 os.system('python run_profile.py')
-
+os.system('python show_balanced_limiters.py')
Index: /inundation/ga/storm_surge/pyvolution/domain.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/domain.py	(revision 239)
+++ /inundation/ga/storm_surge/pyvolution/domain.py	(revision 240)
@@ -70,6 +70,7 @@
 
         #Defaults
-        from config import max_smallsteps, beta, newstyle
+        from config import max_smallsteps, beta, newstyle, epsilon
         self.beta = beta
+        self.epsilon = epsilon        
         self.newstyle = newstyle
         self.default_order = 1
@@ -514,47 +515,4 @@
 
 
-
-def compute_fluxes_c(domain, max_timestep):
-    """Compute all fluxes and the timestep suitable for all volumes
-    This version works directly with consecutive data structure
-    and calls C-extension.
-    Note that flux function is hardwired into C-extension.
-    """
-    
-    import sys
-
-    neighbours = Volume.neighbours
-    neighbour_faces = Volume.neighbour_faces
-    normals = Volume.normals
-    flux = Volume.explicit_update    
-
-    area = Volume.geometric[:,0]            
-    radius = Volume.geometric[:,1]            
-    edgelengths = Volume.geometric[:,2:]
-
-
-    flux[:] = 0.0 #Reset stored fluxes to zero        
-    from domain_ext import compute_fluxes
-    timestep = compute_fluxes(domain.flux_function,
-                              Volume.conserved_quantities_face0,
-                              Volume.conserved_quantities_face1,
-                              Volume.conserved_quantities_face2,
-                              Volume.field_values_face0,
-                              Volume.field_values_face1,
-                              Volume.field_values_face2,
-                              Boundary_value.conserved_quantities,
-                              Boundary_value.field_values,
-                              Vector.coordinates,
-                              neighbours,
-                              neighbour_faces,
-                              normals,
-                              flux, area, radius, edgelengths,
-                              max_timestep)
-
-    domain.timestep = timestep
-
-
-
-    
                 
 ##############################################
@@ -571,2 +529,21 @@
     #from python_versions import distribute_to_vertices_and_edges
     #from python_versions import update_conserved_quantities
+
+
+#Optimisation with psyco
+from config import use_psyco
+if use_psyco:
+    try:
+        import psyco
+    except:
+        msg = 'WARNING: psyco (speedup) could not import'+\
+              ', you may want to consider installing it'
+        print msg
+    else:
+        psyco.bind(Domain.distribute_to_vertices_and_edges)
+        psyco.bind(Domain.update_boundary) 
+        psyco.bind(Domain.update_timestep)     #Not worth it
+        psyco.bind(Domain.update_conserved_quantities)
+
+if __name__ == "__main__":
+    pass
Index: /inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py	(revision 239)
+++ /inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py	(revision 240)
@@ -2,5 +2,4 @@
 """
 
-from util import rotate
 
 class Boundary:
Index: /inundation/ga/storm_surge/pyvolution/shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 239)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 240)
@@ -30,6 +30,7 @@
                                 conserved_quantities, other_quantities)
 
-        from config import minimum_allowed_height
+        from config import minimum_allowed_height, g
         self.minimum_allowed_height = minimum_allowed_height
+        self.g = g
 
         self.forcing_terms.append(gravity)
@@ -77,4 +78,52 @@
         distribute_to_vertices_and_edges(self)
         
+
+
+#Rotation of momentum vector
+def rotate(q, normal, direction = 1):
+    """Rotate the momentum component q (q[1], q[2])
+    from x,y coordinates to coordinates based on normal vector.
+
+    If direction is negative the rotation is inverted.
+    
+    Input vector is preserved
+
+    This function is specific to the shallow water wave equation
+    """
+
+    #FIXME: Needs to be tested
+
+    from Numeric import zeros, Float
+    
+    assert len(q) == 3,\
+           'Vector of conserved quantities must have length 3'\
+           'for 2D shallow water equation'
+
+    try:
+        l = len(normal)
+    except:
+        raise 'Normal vector must be an Numeric array'
+
+    #FIXME: Put this test into C-extension as well
+    assert l == 2, 'Normal vector must have 2 components'
+
+ 
+    n1 = normal[0]
+    n2 = normal[1]    
+    
+    r = zeros(len(q), Float) #Rotated quantities
+    r[0] = q[0]              #First quantity, height, is not rotated
+
+    if direction == -1:
+        n2 = -n2
+
+
+    r[1] =  n1*q[1] + n2*q[2]
+    r[2] = -n2*q[1] + n1*q[2]
+    
+    return r
+
+
+
 ####################################
 # Flux computation        
@@ -93,7 +142,4 @@
 
     Bed elevations zl and zr.
-    #FIXME: Remove those and pass in height directly
-    #(unless we'll include a bed elevation discontinuity later)
-    
     """
 
@@ -101,5 +147,4 @@
     from math import sqrt
     from Numeric import array
-    from util import rotate
 
     #Align momentums with x-axis
@@ -108,5 +153,4 @@
 
     z = (zl+zr)/2 #Take average of field values
-    #FIXME: Maybe allow discontinuity later
 
     w_left  = q_left[0]   #w=h+z
@@ -174,5 +218,6 @@
     Fluxes across each edge are scaled by edgelengths and summed up
     Resulting flux is then scaled by area and stored in
-    domain.explicit_update
+    explicit_update for each of the three conserved quantities 
+    level, xmomentum and ymomentum    
 
     The maximal allowable speed computed by the flux_function for each volume
@@ -187,18 +232,7 @@
     import sys
     from Numeric import zeros, Float
-    from config import max_timestep    
 
     N = domain.number_of_elements
     
-    neighbours = domain.neighbours
-    neighbour_edges = domain.neighbour_edges
-    normals = domain.normals
-
-    areas = domain.areas
-    radii = domain.radii
-    edgelengths = domain.edgelengths
-        
-    timestep = max_timestep #FIXME: Get rid of this
-
     #Shortcuts
     Level = domain.quantities['level']
@@ -220,6 +254,6 @@
 
     #Loop
+    timestep = float(sys.maxint)    
     for k in range(N):
-        optimal_timestep = float(sys.maxint)
 
         flux[:] = 0.  #Reset work array
@@ -230,11 +264,11 @@
 
             #Quantities at neighbour on nearest face
-            n = neighbours[k,i] 
+            n = domain.neighbours[k,i] 
             if n < 0:
-                m = -n-1 #Convert neg flag to index
+                m = -n-1 #Convert negative flag to index
                 qr = [level_bdry[m], xmom_bdry[m], ymom_bdry[m]]
                 zr = zl #Extend bed elevation to boundary
             else:    
-                m = neighbour_edges[k,i]
+                m = domain.neighbour_edges[k,i]
                 qr = [level[n, m], xmom[n, m], ymom[n, m]]
                 zr = bed[n, m]                
@@ -242,13 +276,13 @@
                 
             #Outward pointing normal vector   
-            normal = normals[k, 2*i:2*i+2]
+            normal = domain.normals[k, 2*i:2*i+2]
 
             #Flux computation using provided function
             edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-            flux -= edgeflux * edgelengths[k,i]
-            
+            flux -= edgeflux * domain.edgelengths[k,i]
+
             #Update optimal_timestep
             try:
-                optimal_timestep = min(optimal_timestep, radii[k]/max_speed)
+                timestep = min(timestep, domain.radii[k]/max_speed)
             except ZeroDivisionError:
                 pass
@@ -256,14 +290,47 @@
         #Normalise by area and store for when all conserved
         #quantities get updated
-        flux /= areas[k]
+        flux /= domain.areas[k]
         Level.explicit_update[k] = flux[0]
         Xmom.explicit_update[k] = flux[1]
         Ymom.explicit_update[k] = flux[2]
 
-        timestep = min(timestep, optimal_timestep)
-
     domain.timestep = timestep    
 
-	
+
+def compute_fluxes_c(domain):
+    """Wrapper calling c version of compute fluxes
+    """
+
+    import sys
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+    
+    #Shortcuts
+    Level = domain.quantities['level']
+    Xmom = domain.quantities['xmomentum']
+    Ymom = domain.quantities['ymomentum']
+    Bed = domain.quantities['elevation']        
+
+    timestep = float(sys.maxint)    
+    from shallow_water_ext import compute_fluxes
+    domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
+                                     domain.neighbours,
+                                     domain.neighbour_edges,
+                                     domain.normals,
+                                     domain.edgelengths,                       
+                                     domain.radii,
+                                     domain.areas,
+                                     Level.edge_values, 
+                                     Xmom.edge_values, 
+                                     Ymom.edge_values,  
+                                     Bed.edge_values,   
+                                     Level.boundary_values,
+                                     Xmom.boundary_values,
+                                     Ymom.boundary_values,
+                                     Level.explicit_update,
+                                     Xmom.explicit_update,
+                                     Ymom.explicit_update)
+        
 
 ####################################
@@ -838,2 +905,30 @@
 #Initialise module
 
+
+import compile
+if compile.can_use_C_extension('shallow_water_ext.c'):
+    #Replace python version with c implementations
+    from shallow_water_ext import rotate
+    compute_fluxes = compute_fluxes_c
+    #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c
+    #update_conserved_quantities = update_conserved_quantities_c    
+
+
+#Optimisation with psyco
+from config import use_psyco
+if use_psyco:
+    try:
+        import psyco
+    except:
+        msg = 'WARNING: psyco (speedup) could not import'+\
+              ', you may want to consider installing it'
+        print msg
+    else:
+        psyco.bind(distribute_to_vertices_and_edges)
+        psyco.bind(compute_fluxes_c)  
+        #psyco.bind(update_boundary_values) 
+        #psyco.bind(Domain.update_timestep)     #Not worth it
+        #psyco.bind(update_conserved_quantities)
+
+if __name__ == "__main__":
+    pass
Index: /inundation/ga/storm_surge/pyvolution/shallow_water_ext.c
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water_ext.c	(revision 240)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water_ext.c	(revision 240)
@@ -0,0 +1,424 @@
+// Python - C extension module for shallow_water.py
+//
+// To compile (Python2.3):
+//  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O 
+//  gcc -shared domain_ext.o  -o domain_ext.so  
+//
+// or use python compile.py
+//
+// See the module shallow_water.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;
+}
+
+
+// Computational function for rotation
+// FIXME: Should probably go into shallow_water_ext.c
+int _rotate(double *q, double n1, double n2) {
+  /*Rotate the momentum component q (q[1], q[2])
+    from x,y coordinates to coordinates based on normal vector (n1, n2).
+    
+    Result is returned in array 3x1 r     
+    To rotate in opposite direction, call rotate with (q, n1, -n2)
+    
+    Contents of q are changed by this function */    
+
+
+  double q1, q2;
+  
+  //Shorthands
+  q1 = q[1];  //uh momentum
+  q2 = q[2];  //vh momentum
+
+  //Rotate
+  q[1] =  n1*q1 + n2*q2;
+  q[2] = -n2*q1 + n1*q2; 
+
+  return 0;
+}  
+
+
+	
+// Computational function for flux computation (using stage w=z+h)
+int flux_function(double *q_left, double *q_right, 
+		  double z_left, double z_right, 
+		  double n1, double n2, 
+		  double epsilon, double g, 
+		  double *edgeflux, double *max_speed) {
+  
+  /*Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of the 'stage', 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 
+  */
+       
+  int i;
+       
+  double w_left, h_left, uh_left, vh_left, u_left;
+  double w_right, h_right, uh_right, vh_right, u_right;
+  double s_min, s_max, soundspeed_left, soundspeed_right;
+  double denom, z;
+  double q_left_copy[3], q_right_copy[3];    
+  double flux_right[3], flux_left[3];
+  
+  //Copy conserved quantities to protect from modification
+  for (i=0; i<3; i++) {
+    q_left_copy[i] = q_left[i];
+    q_right_copy[i] = q_right[i];
+  } 
+    
+  //Align x- and y-momentum with x-axis
+  _rotate(q_left_copy, n1, n2);
+  _rotate(q_right_copy, n1, n2);    
+
+  z = (z_left+z_right)/2; //Take average of field values
+
+  //Compute speeds in x-direction
+  w_left = q_left_copy[0];              // h+z
+  h_left = w_left-z;
+  uh_left = q_left_copy[1];
+
+  if (h_left < epsilon) {
+    h_left = 0.0;  //Could have been negative
+    u_left = 0.0;
+  } else { 
+    u_left = uh_left/h_left;
+  }
+  
+  w_right = q_right_copy[0];
+  h_right = w_right-z;
+  uh_right = q_right_copy[1];
+
+  if (h_right < epsilon) {
+    h_right = 0.0; //Could have been negative
+    u_right = 0.0;
+  } else { 
+    u_right = uh_right/h_right;
+  }
+
+  //Momentum in y-direction              
+  vh_left  = q_left_copy[2];
+  vh_right = q_right_copy[2];   
+       
+
+  //Maximal and minimal wave speeds
+  soundspeed_left  = sqrt(g*h_left); 
+  soundspeed_right = sqrt(g*h_right);
+    
+  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
+  if (s_max < 0.0) s_max = 0.0; 
+  
+  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
+  if (s_min > 0.0) s_min = 0.0;   
+  
+  //Flux formulas  
+  flux_left[0] = u_left*h_left;
+  flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
+  flux_left[2] = u_left*vh_left;
+  
+  flux_right[0] = u_right*h_right;
+  flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
+  flux_right[2] = u_right*vh_right;
+    
+
+  //Flux computation    
+  denom = s_max-s_min;
+  if (denom == 0.0) {
+    for (i=0; i<3; i++) edgeflux[i] = 0.0;
+    *max_speed = 0.0;
+  } else {    
+    for (i=0; i<3; i++) {
+      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
+      edgeflux[i] += s_max*s_min*(q_right_copy[i]-q_left_copy[i]);
+      edgeflux[i] /= denom;
+    }  
+        
+    //Maximal wavespeed
+    *max_speed = max(fabs(s_max), fabs(s_min));
+    
+    //Rotate back       
+    _rotate(edgeflux, n1, -n2);
+  }
+  return 0;
+}
+	
+
+///////////////////////////////////////////////////////////////////
+// Gateways to Python
+
+
+
+
+PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
+  //
+  // r = rotate(q, normal, direction=1)
+  //
+  // Where q is assumed to be a Float numeric array of length 3 and
+  // normal a Float numeric array of length 2.
+
+  
+  PyObject *Q, *Normal;
+  PyArrayObject *q, *r, *normal;
+  
+  static char *argnames[] = {"q", "normal", "direction", NULL};
+  int dimensions[1], i, direction=1;
+  double n1, n2;
+
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames, 
+				   &Q, &Normal, &direction))  
+    return NULL;  
+
+  //Input checks (convert sequences into numeric arrays)
+  q = (PyArrayObject *) 
+    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
+  normal = (PyArrayObject *) 
+    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
+  
+  //Allocate space for return vector r (don't DECREF) 
+  dimensions[0] = 3;
+  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
+
+  //Copy
+  for (i=0; i<3; i++) {
+    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];  
+  }
+  
+  //Get normal and direction
+  n1 = ((double *) normal -> data)[0]; 
+  n2 = ((double *) normal -> data)[1];   
+  if (direction == -1) n2 = -n2;
+
+  //Rotate
+  _rotate((double *) r -> data, n1, n2);
+  
+  //return result using PyArray to avoid memory leak
+  return PyArray_Return(r);
+}    
+
+
+
+
+
+PyObject *compute_fluxes(PyObject *self, PyObject *args) {
+  /*Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities 
+    level, xmomentum and ymomentum
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+    
+    Python call:
+    domain.timestep = compute_fluxes(timestep, 
+                                     domain.epsilon, 
+				     domain.g,
+                                     domain.neighbours,
+                                     domain.neighbour_edges,
+                                     domain.normals,
+                                     domain.edgelengths,                       
+                                     domain.radii,
+                                     domain.areas,
+                                     Level.edge_values, 
+                                     Xmom.edge_values, 
+                                     Ymom.edge_values,  
+                                     Bed.edge_values,   
+                                     Level.boundary_values,
+                                     Xmom.boundary_values,
+                                     Ymom.boundary_values,
+                                     Level.explicit_update,
+                                     Xmom.explicit_update,
+                                     Ymom.explicit_update)
+        
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes. 
+
+            
+  */
+
+  
+  PyArrayObject *neighbours, *neighbour_edges,
+    *normals, *edgelengths, *radii, *areas,
+    *level_edge_values, 
+    *xmom_edge_values, 
+    *ymom_edge_values,  
+    *bed_edge_values,   
+    *level_boundary_values,
+    *xmom_boundary_values,
+    *ymom_boundary_values,
+    *level_explicit_update,
+    *xmom_explicit_update,
+    *ymom_explicit_update;
+
+    
+  //Local variables  
+  double timestep, max_speed, epsilon, g;
+  double normal[2], ql[3], qr[3], zl, zr;
+  double flux[3], edgeflux[3]; //Work arrays for summing up fluxes
+
+  int number_of_elements, k, i, j, m, n;
+  int ki, nm, ki2; //Index shorthands
+  
+ 
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOO",
+			&timestep,
+			&epsilon,
+			&g,
+			&neighbours, 
+			&neighbour_edges,
+			&normals, 
+			&edgelengths, &radii, &areas,
+			&level_edge_values, 
+			&xmom_edge_values, 
+			&ymom_edge_values,  
+			&bed_edge_values,   
+			&level_boundary_values,
+			&xmom_boundary_values,
+			&ymom_boundary_values,
+			&level_explicit_update,
+			&xmom_explicit_update,
+			&ymom_explicit_update)) {
+    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
+    return NULL;
+  }
+
+  number_of_elements = level_edge_values -> dimensions[0];
+  
+    
+  for (k=0; k<number_of_elements; k++) {
+
+    //Reset work array
+    for (j=0; j<3; j++) flux[j] = 0.0;
+			  
+    //Loop through neighbours and compute edge flux for each
+    for (i=0; i<3; i++) {
+      ki = k*3+i;
+      ql[0] = ((double *) level_edge_values -> data)[ki];
+      ql[1] = ((double *) xmom_edge_values -> data)[ki];
+      ql[2] = ((double *) ymom_edge_values -> data)[ki];            
+      zl =    ((double *) bed_edge_values -> data)[ki];                  
+      
+      //Quantities at neighbour on nearest face
+      n = ((int *) neighbours -> data)[ki];
+      if (n < 0) {
+	m = -n-1; //Convert negative flag to index
+	qr[0] = ((double *) level_boundary_values -> data)[m];	
+	qr[1] = ((double *) xmom_boundary_values -> data)[m];	
+	qr[2] = ((double *) ymom_boundary_values -> data)[m];	
+	zr = zl; //Extend bed elevation to boundary
+      } else {    
+	m = ((int *) neighbour_edges -> data)[ki];
+	
+	nm = n*3+m;	
+	qr[0] = ((double *) level_edge_values -> data)[nm];
+	qr[1] = ((double *) xmom_edge_values -> data)[nm];
+	qr[2] = ((double *) ymom_edge_values -> data)[nm];            
+	zr =    ((double *) bed_edge_values -> data)[nm];                  
+      }
+      
+      // Outward pointing normal vector   
+      // normal = domain.normals[k, 2*i:2*i+2]
+      ki2 = 2*ki; //k*6 + i*2
+      normal[0] = ((double *) normals -> data)[ki2];
+      normal[1] = ((double *) normals -> data)[ki2+1];      
+
+      //Edge flux computation
+      flux_function(ql, qr, zl, zr, 
+		    normal[0], normal[1],
+		    epsilon, g, 
+		    edgeflux, &max_speed);
+
+		    
+      //flux -= edgeflux * edgelengths[k,i]
+      for (j=0; j<3; j++) { 
+	flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
+      }
+      
+      //Update timestep
+      //timestep = min(timestep, domain.radii[k]/max_speed)
+      if (max_speed > epsilon) {
+	timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
+      }   
+    } // end for i
+    
+    //Normalise by area and store for when all conserved
+    //quantities get updated
+    // flux /= areas[k]
+    for (j=0; j<3; j++) { 
+      flux[j] /= ((double *) areas -> data)[k];
+    }
+
+    ((double *) level_explicit_update -> data)[k] = flux[0];
+    ((double *) xmom_explicit_update -> data)[k] = flux[1];
+    ((double *) ymom_explicit_update -> data)[k] = flux[2];        
+
+  } //end for k
+
+  return Py_BuildValue("d", timestep);
+}    
+
+
+
+
+
+
+
+
+//////////////////////////////////////////	
+// Method table for python module
+static struct PyMethodDef MethodTable[] = {
+  /* The cast of the function is necessary since PyCFunction values
+   * only take two PyObject* parameters, and rotate() takes
+   * three.
+   */
+  
+  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
+  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},    
+  //{"distribute_to_vertices_and_edges", 
+  // distribute_to_vertices_and_edges, METH_VARARGS},    
+  //{"update_conserved_quantities", 
+  // update_conserved_quantities, METH_VARARGS},       
+  //{"set_initialcondition", 
+  // set_initialcondition, METH_VARARGS},    
+  {NULL, NULL}
+};
+	
+// Module initialisation   
+void initshallow_water_ext(void){
+  Py_InitModule("shallow_water_ext", MethodTable);
+  
+  import_array();     //Necessary for handling of NumPY structures  
+}
Index: /inundation/ga/storm_surge/pyvolution/test_shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 239)
+++ /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 240)
@@ -17,4 +17,18 @@
         pass
         #print "  Tearing down"
+
+    def test_rotate(self):
+        normal = array([0.0,-1.0])        
+
+        q = array([1.0,2.0,3.0])
+      
+        r = rotate(q, normal, direction = 1)
+        assert r[0] == 1
+        assert r[1] == -3
+        assert r[2] == 2
+
+        w = rotate(r, normal, direction = -1)        
+        assert allclose(w, q)
+
 
 
@@ -189,5 +203,5 @@
         
 
-    def test_compute_fluxes(self):
+    def test_compute_fluxes0(self):
         #Do a full triangle and check that fluxes cancel out for
         #the constant level case
Index: /inundation/ga/storm_surge/pyvolution/test_util.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 239)
+++ /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 240)
@@ -27,17 +27,4 @@
         assert zy == 0.0
 
-
-    def test_rotate(self):
-        normal = array([0.0,-1.0])        
-
-        q = array([1.0,2.0,3.0])
-      
-        r = rotate(q, normal, direction = 1)
-        assert r[0] == 1
-        assert r[1] == -3
-        assert r[2] == 2
-
-        w = rotate(r, normal, direction = -1)        
-        assert allclose(w, q)
 
 
Index: /inundation/ga/storm_surge/pyvolution/util.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/util.py	(revision 239)
+++ /inundation/ga/storm_surge/pyvolution/util.py	(revision 240)
@@ -82,48 +82,4 @@
 
 
-def rotate_python(q, normal, direction = 1):
-    """Rotate the momentum component q (q[1], q[2])
-    from x,y coordinates to coordinates based on normal vector.
-
-    If direction is negative the rotation is inverted.
-    
-    Input vector is preserved
-
-    This function is specific to the shallow water wave equation
-    """
-
-    #FIXME: Needs to be tested
-
-    from Numeric import zeros, Float
-    
-    assert len(q) == 3,\
-           'Vector of conserved quantities must have length 3'\
-           'for 2D shallow water equation'
-
-    try:
-        l = len(normal)
-    except:
-        raise 'Normal vector must be an Numeric array'
-
-    #FIXME: Put this test into C-extension as well
-    assert l == 2, 'Normal vector must have 2 components'
-
- 
-    n1 = normal[0]
-    n2 = normal[1]    
-    
-    r = zeros(len(q), Float) #Rotated quantities
-    r[0] = q[0]              #First quantity, height, is not rotated
-
-    if direction == -1:
-        n2 = -n2
-
-
-    r[1] =  n1*q[1] + n2*q[2]
-    r[2] = -n2*q[1] + n1*q[2]
-    
-    return r
-
-
 
 
@@ -133,10 +89,7 @@
 import compile
 if compile.can_use_C_extension('util_ext.c'):
-    from util_ext import gradient, rotate
+    from util_ext import gradient
 else:
     gradient = gradient_python
-    rotate = rotate_python
-
-
 
 
Index: /inundation/ga/storm_surge/pyvolution/util_ext.c
===================================================================
--- /inundation/ga/storm_surge/pyvolution/util_ext.c	(revision 239)
+++ /inundation/ga/storm_surge/pyvolution/util_ext.c	(revision 240)
@@ -36,39 +36,4 @@
 
 
-// Computational function for rotation
-int _rotate(double *q, double *r, double *normal, int direction) {
-  /*Rotate the momentum component q (q[1], q[2])
-    from x,y coordinates to coordinates based on normal vector.
-    
-    To rotate in opposite direction, call rotate with direction=1
-    Result is returned in r 
-    
-    q and r are assumed to have three components each*/
-
-
-  double n0, n1, q1, q2;
-  
-  //Determine normal and direction
-  n0 = normal[0];
-  if (direction == 1) {
-    n1 = normal[1];    
-  } else {
-    n1 = -normal[1];
-  }
-
-  //Shorthands
-  q1 = q[1];  //uh momentum
-  q2 = q[2];  //vh momentum
-
-  //Rotate
-  r[0] = q[0];
-  r[1] = n0*q1 + n1*q2;
-  r[2] = n0*q2 - n1*q1;
-
-  return 0;
-}  
-
-
-
 // Gateway to Python
 PyObject *gradient(PyObject *self, PyObject *args) {
@@ -96,61 +61,4 @@
 
 
-// Gateway to Python
-PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
-  //
-  // r = rotate(q, normal, direction=0)
-  //
-  // Where q is assumed to be a Float numeric array of length 3 and
-  // normal a Float numeric array of length 2.
-
-  //FIXME: Input checks and unit tests
-  
-  PyObject *Q, *Normal;
-  PyArrayObject *q, *r, *normal;
-  
-  static char *argnames[] = {"q", "normal", "direction", NULL};
-  int dimensions[1];
-  int N, direction=0;
-
-  // Convert Python arguments to C  
-  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames, 
-				   &Q, &Normal, &direction))  
-    return NULL;  
-
-  //Input checks (convert sequenced into numeric arrays)
-  
-  q = (PyArrayObject *) 
-    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
-  normal = (PyArrayObject *) 
-    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
-  
-  //Allocate space for return vector r (don't DECREF) 
-  N = q -> dimensions[0];
-  dimensions[0] = N;
-  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
-
-  //Rotate
-  _rotate((double *) q -> data,
-	  (double *) r -> data,
-	  (double *) normal -> data,	  
-	  direction);	  
-
-  //Build result and DECREF r - returning r directly causes memory leak
-  //result = Py_BuildValue("O", r);	  
-  //Py_DECREF(r);
-  
-  //return result;
-  return PyArray_Return(r);
-}    
-
-
-// Method table for python module
-//static struct PyMethodDef MethodTable[] = {
-//  {"gradient", gradient, METH_VARARGS},  
-//  {NULL, NULL}
-//};
-
-
-
 // Method table for python module
 static struct PyMethodDef MethodTable[] = {
@@ -160,5 +68,5 @@
    */
   
-  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},  
+  //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},  
   {"gradient", gradient, METH_VARARGS, "Print out"},    
   {NULL, NULL, 0, NULL}   /* sentinel */
