Index: inundation/utilities/__init__.py
===================================================================
--- inundation/utilities/__init__.py	(revision 1910)
+++ inundation/utilities/__init__.py	(revision 1910)
@@ -0,0 +1,5 @@
+"""Make directory available as a Python package
+"""
+
+pass
+
Index: inundation/utilities/compile.py
===================================================================
--- inundation/utilities/compile.py	(revision 1908)
+++ inundation/utilities/compile.py	(revision 1910)
@@ -240,5 +240,5 @@
     """
 
-    from config import use_extensions
+    #from config import use_extensions
 
     from os.path import splitext
@@ -247,34 +247,34 @@
     
     C=False
-    if use_extensions:
+
+    try:
+        s = 'import %s' %root
+        #print s
+        exec(s)
+    except:
         try:
-            s = 'import %s' %root
-            #print s
-            exec(s)
+            open(filename)
         except:
+            msg = 'C extension %s cannot be opened' %filename
+            print msg                
+        else:    
+            print '------- Trying to compile c-extension %s' %filename
+        
             try:
-                open(filename)
+                compile(filename)
             except:
-                msg = 'C extension %s cannot be opened' %filename
-                print msg                
-            else:    
-                print '------- Trying to compile c-extension %s' %filename
-	    
+                print 'WARNING: Could not compile C-extension %s'\
+                      %filename
+            else:
                 try:
-                    compile(filename)
+                    exec('import %s' %root)
                 except:
-                    print 'WARNING: Could not compile C-extension %s'\
-                          %filename
+                    msg = 'C extension %s seems to compile OK, '
+                    msg += 'but it can still not be imported.'
+                    raise msg
                 else:
-                    try:
-                        exec('import %s' %root)
-                    except:
-                        msg = 'C extension %s seems to compile OK, '
-                        msg += 'but it can still not be imported.'
-                        raise msg
-                    else:
-                        C=True
-        else:
-            C=True
+                    C=True
+    else:
+        C=True
             
     if not C:
Index: inundation/utilities/numerical_tools.py
===================================================================
--- inundation/utilities/numerical_tools.py	(revision 1910)
+++ inundation/utilities/numerical_tools.py	(revision 1910)
@@ -0,0 +1,35 @@
+#!/usr/bin/env python
+"""Auxiliary numerical tools
+
+"""
+
+	
+def ensure_numeric(A, typecode = None):
+    """Ensure that sequence is a Numeric array.
+    Inputs:
+        A: Sequence. If A is already a Numeric array it will be returned
+                     unaltered
+                     If not, an attempt is made to convert it to a Numeric
+                     array
+        typecode: Numeric type. If specified, use this in the conversion.
+                                If not, let Numeric decide
+
+    This function is necessary as array(A) can cause memory overflow.
+    """
+
+    from Numeric import ArrayType, array
+
+    if typecode is None:
+        if type(A) == ArrayType:
+            return A
+        else:
+            return array(A)
+    else:
+        if type(A) == ArrayType:
+            if A.typecode == typecode:
+                return array(A)  #FIXME: Shouldn't this just return A?
+            else:
+                return A.astype(typecode)
+        else:
+            return array(A).astype(typecode)
+
Index: inundation/utilities/polygon.py
===================================================================
--- inundation/utilities/polygon.py	(revision 1908)
+++ inundation/utilities/polygon.py	(revision 1910)
@@ -1,12 +1,8 @@
-
-
-
-
-#####################################
-#POLYGON STUFF
-#
-#FIXME: All these should be put into new module polygon.py
-
-
+#!/usr/bin/env python
+"""Polygon manipulations
+
+"""
+
+from utilities.numerical_tools import ensure_numeric
 
 def point_on_line(x, y, x0, y0, x1, y1):
@@ -284,5 +280,5 @@
     M = points.shape[0]  #Number of points
 
-    from util_ext import separate_points_by_polygon
+    from polygon_ext import separate_points_by_polygon
 
     if verbose: print 'Allocating array for indices'
@@ -351,5 +347,4 @@
 
     def __call__(self, x, y):
-        from util import inside_polygon
 	from Numeric import ones, Float, concatenate, array, reshape, choose
 
Index: inundation/utilities/polygon_ext.c
===================================================================
--- inundation/utilities/polygon_ext.c	(revision 1908)
+++ inundation/utilities/polygon_ext.c	(revision 1910)
@@ -245,48 +245,4 @@
 
 
-PyObject *gradient(PyObject *self, PyObject *args) {
-  //
-  // a,b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
-  //
-
-  double x0, y0, x1, y1, x2, y2, q0, q1, q2, a, b;
-  PyObject *result;
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "ddddddddd", &x0, &y0, &x1, &y1, &x2, &y2,
-			&q0, &q1, &q2))
-    return NULL;
-
-
-  // Call underlying routine
-  _gradient(x0, y0, x1, y1, x2, y2,
-	    q0, q1, q2, &a, &b);
-
-  // Return values a and b
-  result = Py_BuildValue("dd", a, b);
-  return result;
-}
-
-PyObject *gradient2(PyObject *self, PyObject *args) {
-  //
-  // a,b = gradient2(x0, y0, x1, y1, q0, q1)
-  //
-
-  double x0, y0, x1, y1, q0, q1, a, b;
-  PyObject *result;
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dddddd", &x0, &y0, &x1, &y1, &q0, &q1))
-    return NULL;
-
-
-  // Call underlying routine
-  _gradient2(x0, y0, x1, y1, q0, q1, &a, &b);
-
-  // Return values a and b
-  result = Py_BuildValue("dd", a, b);
-  return result;
-}
-
 // Method table for python module
 static struct PyMethodDef MethodTable[] = {
@@ -296,7 +252,4 @@
    */
 
-  //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
-  {"gradient", gradient, METH_VARARGS, "Print out"},
-  {"gradient2", gradient2, METH_VARARGS, "Print out"},  
   {"point_on_line", point_on_line, METH_VARARGS, "Print out"},
   {"separate_points_by_polygon", separate_points_by_polygon, 
@@ -308,5 +261,5 @@
 
 // Module initialisation
-void initutil_ext(void){
+void initpolygon_ext(void){
   Py_InitModule("polygon_ext", MethodTable);
 
Index: inundation/utilities/test_numerical_tools.py
===================================================================
--- inundation/utilities/test_numerical_tools.py	(revision 1910)
+++ inundation/utilities/test_numerical_tools.py	(revision 1910)
@@ -0,0 +1,75 @@
+#!/usr/bin/env python
+
+
+import unittest
+from Numeric import zeros, array, allclose
+from math import sqrt, pi
+
+from numerical_tools import *
+
+def test_function(x, y):
+    return x+y
+
+class Test_Numerical_Tools(unittest.TestCase):
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+
+
+	
+    def test_ensure_numeric(self):
+        from numerical_tools import ensure_numeric
+        from Numeric import ArrayType, Float, array
+
+        A = [1,2,3,4]
+        B = ensure_numeric(A)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'l'
+        assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4
+
+
+        A = [1,2,3.14,4]
+        B = ensure_numeric(A)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'd'
+        assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4
+
+
+        A = [1,2,3,4]
+        B = ensure_numeric(A, Float)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'd'
+        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
+
+
+        A = [1,2,3,4]
+        B = ensure_numeric(A, Float)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'd'
+        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
+
+
+        A = array([1,2,3,4])
+        B = ensure_numeric(A)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'l'        
+        assert A == B    
+        assert A is B   #Same object
+
+
+        A = array([1,2,3,4])
+        B = ensure_numeric(A, Float)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'd'        
+        assert A == B    
+        assert A is not B   #Not the same object        
+
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Numerical_Tools,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: inundation/utilities/test_polygon.py
===================================================================
--- inundation/utilities/test_polygon.py	(revision 1910)
+++ inundation/utilities/test_polygon.py	(revision 1910)
@@ -0,0 +1,350 @@
+#!/usr/bin/env python
+
+
+import unittest
+from Numeric import zeros, array, allclose
+from math import sqrt, pi
+
+from polygon import *
+
+
+def test_function(x, y):
+    return x+y
+
+class Test_Polygon(unittest.TestCase):
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+
+    def test_that_C_extension_compiles(self):
+        FN = 'polygon_ext.c'
+        try:
+            import polygon_ext
+        except:
+            from compile import compile
+
+            try:
+                compile(FN)
+            except:
+                raise 'Could not compile %s' %FN
+            else:
+                import polygon_ext
+
+
+    #Polygon stuff
+    def test_polygon_function_constants(self):
+        p1 = [[0,0], [10,0], [10,10], [0,10]]
+        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
+
+	f = Polygon_function( [(p1, 1.0)] )
+	z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
+	assert allclose(z, [1,1,0,0])
+
+
+	f = Polygon_function( [(p2, 2.0)] )
+	z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
+	assert allclose(z, [2,0,0,2])
+
+
+	#Combined
+	f = Polygon_function( [(p1, 1.0), (p2, 2.0)] )
+	z = f([5, 5, 27, 35], [5, 9, 8, -5])
+	assert allclose(z, [2,1,0,2])
+
+
+    def test_polygon_function_callable(self):
+        """Check that values passed into Polygon_function can be callable
+	themselves.
+	"""
+        p1 = [[0,0], [10,0], [10,10], [0,10]]
+        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
+
+	f = Polygon_function( [(p1, test_function)] )
+	z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
+	assert allclose(z, [10,14,0,0])
+
+	#Combined
+	f = Polygon_function( [(p1, test_function), (p2, 2.0)] )
+	z = f([5, 5, 27, 35], [5, 9, 8, -5])
+	assert allclose(z, [2,14,0,2])
+
+
+	#Combined w default
+	f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14)
+	z = f([5, 5, 27, 35], [5, 9, 8, -5])
+	assert allclose(z, [2,14,3.14,2])
+
+
+	#Combined w default func
+	f = Polygon_function( [(p1, test_function), (p2, 2.0)],
+			      default = test_function)
+	z = f([5, 5, 27, 35], [5, 9, 8, -5])
+	assert allclose(z, [2,14,35,2])
+
+
+    def test_point_on_line(self):
+
+	#Endpoints first
+	assert point_on_line( 0, 0, 0,0, 1,0 )
+	assert point_on_line( 1, 0, 0,0, 1,0 )
+
+	#Then points on line
+	assert point_on_line( 0.5, 0, 0,0, 1,0 )
+	assert point_on_line( 0, 0.5, 0,1, 0,0 )
+	assert point_on_line( 1, 0.5, 1,1, 1,0 )
+	assert point_on_line( 0.5, 0.5, 0,0, 1,1 )
+
+	#Then points not on line
+	assert not point_on_line( 0.5, 0, 0,1, 1,1 )
+	assert not point_on_line( 0, 0.5, 0,0, 1,1 )
+
+	#From real example that failed
+	assert not point_on_line( 40,50, 40,20, 40,40 )
+
+
+	#From real example that failed
+	assert not point_on_line( 40,19, 40,20, 40,40 )
+
+
+
+
+    def test_inside_polygon_main(self):
+
+
+        #Simplest case: Polygon is the unit square
+        polygon = [[0,0], [1,0], [1,1], [0,1]]
+
+	assert inside_polygon( (0.5, 0.5), polygon )
+	assert not inside_polygon( (0.5, 1.5), polygon )
+	assert not inside_polygon( (0.5, -0.5), polygon )
+	assert not inside_polygon( (-0.5, 0.5), polygon )
+	assert not inside_polygon( (1.5, 0.5), polygon )
+
+	#Try point on borders
+	assert inside_polygon( (1., 0.5), polygon, closed=True)
+	assert inside_polygon( (0.5, 1), polygon, closed=True)
+	assert inside_polygon( (0., 0.5), polygon, closed=True)
+	assert inside_polygon( (0.5, 0.), polygon, closed=True)
+
+	assert not inside_polygon( (0.5, 1), polygon, closed=False)
+	assert not inside_polygon( (0., 0.5), polygon, closed=False)
+	assert not inside_polygon( (0.5, 0.), polygon, closed=False)
+	assert not inside_polygon( (1., 0.5), polygon, closed=False)
+
+
+
+        #From real example (that failed)
+	polygon = [[20,20], [40,20], [40,40], [20,40]]
+	points = [ [40, 50] ]
+	res = inside_polygon(points, polygon)
+	assert len(res) == 0
+
+	polygon = [[20,20], [40,20], [40,40], [20,40]]
+        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
+	res = inside_polygon(points, polygon)
+	assert len(res) == 2
+	assert allclose(res, [0,1])
+
+
+
+	#More convoluted and non convex polygon
+        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
+	assert inside_polygon( (0.5, 0.5), polygon )
+	assert inside_polygon( (1, -0.5), polygon )
+	assert inside_polygon( (1.5, 0), polygon )
+
+	assert not inside_polygon( (0.5, 1.5), polygon )
+	assert not inside_polygon( (0.5, -0.5), polygon )
+
+
+	#Very convoluted polygon
+        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
+	assert inside_polygon( (5, 5), polygon )
+	assert inside_polygon( (17, 7), polygon )
+	assert inside_polygon( (27, 2), polygon )
+	assert inside_polygon( (35, -5), polygon )
+	assert not inside_polygon( (15, 7), polygon )
+	assert not inside_polygon( (24, 3), polygon )
+	assert not inside_polygon( (25, -10), polygon )
+
+
+
+	#Another combination (that failed)
+        polygon = [[0,0], [10,0], [10,10], [0,10]]
+	assert inside_polygon( (5, 5), polygon )
+	assert inside_polygon( (7, 7), polygon )
+	assert not inside_polygon( (-17, 7), polygon )
+	assert not inside_polygon( (7, 17), polygon )
+	assert not inside_polygon( (17, 7), polygon )
+	assert not inside_polygon( (27, 8), polygon )
+	assert not inside_polygon( (35, -5), polygon )
+
+
+
+
+	#Multiple polygons
+
+        polygon = [[0,0], [1,0], [1,1], [0,1], [0,0],
+		   [10,10], [11,10], [11,11], [10,11], [10,10]]
+        assert inside_polygon( (0.5, 0.5), polygon )
+        assert inside_polygon( (10.5, 10.5), polygon )
+
+	#FIXME: Fails if point is 5.5, 5.5
+        assert not inside_polygon( (0, 5.5), polygon )
+
+	#Polygon with a hole
+        polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
+	           [0,0], [1,0], [1,1], [0,1], [0,0]]
+
+        assert inside_polygon( (0, -0.5), polygon )
+        assert not inside_polygon( (0.5, 0.5), polygon )
+
+    def test_inside_polygon_vector_version(self):
+	#Now try the vector formulation returning indices
+        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
+	points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+	res = inside_polygon( points, polygon, verbose=False )
+
+	assert allclose( res, [0,1,2] )
+
+    def test_outside_polygon(self):
+        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square    
+
+        assert not outside_polygon( [0.5, 0.5], U )
+        #evaluate to False as the point 0.5, 0.5 is inside the unit square
+        
+        assert outside_polygon( [1.5, 0.5], U )
+        #evaluate to True as the point 1.5, 0.5 is outside the unit square
+        
+        indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
+        assert allclose( indices, [1] )
+        
+        #One more test of vector formulation returning indices
+        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
+	points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+	res = outside_polygon( points, polygon )
+
+	assert allclose( res, [3, 4] )
+
+
+
+        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
+	points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+	res = outside_polygon( points, polygon )
+
+	assert allclose( res, [0, 4, 5] )        
+     
+    def test_outside_polygon2(self):
+        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square    
+   
+        assert not outside_polygon( [0.5, 1.0], U, closed = True )
+        #evaluate to False as the point 0.5, 1.0 is inside the unit square
+        
+        assert outside_polygon( [0.5, 1.0], U, closed = False )
+        #evaluate to True as the point 0.5, 1.0 is outside the unit square
+
+    def test_separate_points_by_polygon(self):
+        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square    
+
+        indices, count = separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
+        assert allclose( indices, [0,2,1] )
+        assert count == 2
+        
+        #One more test of vector formulation returning indices
+        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
+	points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+	res, count = separate_points_by_polygon( points, polygon )
+
+	assert allclose( res, [0,1,2,4,3] )
+	assert count == 3
+
+
+        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
+	points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+	res, count = separate_points_by_polygon( points, polygon )
+
+	assert allclose( res, [1,2,3,5,4,0] )        
+     	assert count == 3
+	
+
+    def test_populate_polygon(self):
+
+        polygon = [[0,0], [1,0], [1,1], [0,1]]
+        points = populate_polygon(polygon, 5)
+
+        assert len(points) == 5
+        for point in points:
+            assert inside_polygon(point, polygon)
+
+
+	#Very convoluted polygon
+        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
+
+        points = populate_polygon(polygon, 5)
+
+        assert len(points) == 5
+        for point in points:
+            assert inside_polygon(point, polygon)
+
+
+    def test_populate_polygon_with_exclude(self):
+        
+
+        polygon = [[0,0], [1,0], [1,1], [0,1]]
+        ex_poly = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
+        points = populate_polygon(polygon, 5, exclude = [ex_poly])
+
+        assert len(points) == 5
+        for point in points:
+            assert inside_polygon(point, polygon)
+            assert not inside_polygon(point, ex_poly)            
+
+
+        #overlap
+        polygon = [[0,0], [1,0], [1,1], [0,1]]
+        ex_poly = [[-1,-1], [0.5,0], [0.5, 0.5], [-1,0.5]]
+        points = populate_polygon(polygon, 5, exclude = [ex_poly])
+
+        assert len(points) == 5
+        for point in points:
+            assert inside_polygon(point, polygon)
+            assert not inside_polygon(point, ex_poly)                        
+        
+        #Multiple
+        polygon = [[0,0], [1,0], [1,1], [0,1]]
+        ex_poly1 = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
+        ex_poly2 = [[0.5,0.5], [0.5,1], [1, 1], [1,0.5]] #NE quarter        
+        
+        points = populate_polygon(polygon, 20, exclude = [ex_poly1, ex_poly2])
+
+        assert len(points) == 20
+        for point in points:
+            assert inside_polygon(point, polygon)
+            assert not inside_polygon(point, ex_poly1)
+            assert not inside_polygon(point, ex_poly2)                                
+        
+
+	#Very convoluted polygon
+        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
+        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
+        points = populate_polygon(polygon, 20, exclude = [ex_poly])
+        
+        assert len(points) == 20
+        for point in points:
+            assert inside_polygon(point, polygon)
+            assert not inside_polygon(point, ex_poly), '%s' %str(point)                        
+
+
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Polygon,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
+
+
+
+
