Index: inundation/ga/storm_surge/pyvolution/test_util.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_util.py	(revision 998)
+++ inundation/ga/storm_surge/pyvolution/test_util.py	(revision 999)
@@ -826,14 +826,4 @@
         polygon = [[0,0], [1,0], [1,1], [0,1]]
 
-
-        R = inside_polygon( (0.5, 0.5), polygon)
-        print 'Inside = ', R
-        print
-        
-        R = inside_polygon( (0.5, 0.5), polygon)
-        print 'Inside = ', R
-        print
-
-        import sys; sys.exit() 
 	assert inside_polygon( (0.5, 0.5), polygon )
 	assert not inside_polygon( (0.5, 1.5), polygon )   
@@ -953,5 +943,6 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test_inside_polygon_main')
+    #suite = unittest.makeSuite(TestCase,'test_inside_polygon_main')
+    suite = unittest.makeSuite(TestCase,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
Index: inundation/ga/storm_surge/pyvolution/util.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/util.py	(revision 998)
+++ inundation/ga/storm_surge/pyvolution/util.py	(revision 999)
@@ -87,5 +87,6 @@
 
     #FIXME (OLE): Should check origin of domain against that of file
-
+    #In fact, this is where origin should be converted to that of domain
+    #Also, check that file covers domain fully.
     
     assert type(filename) == type(''),\
@@ -974,5 +975,5 @@
     py = polygon[:,1]	
 
-    #Used for an optimisation when points are far away from polgon
+    #Used for an optimisation when points are far away from polygon
     minpx = min(px); maxpx = max(px)
     minpy = min(py); maxpy = max(py)    
@@ -1022,4 +1023,56 @@
 	     
 
+
+def inside_polygon_c(point, polygon, closed = True, verbose = False):
+    """Determine whether points are inside or outside a polygon
+
+    C-wrapper
+    """    
+
+    from Numeric import array, Float, reshape, zeros, Int
+    
+    
+    #Input checks
+    try:
+        point = array(point).astype(Float)   	   
+    except:
+        msg = 'Point %s could not be converted to Numeric array' %point
+	raise msg	 
+	
+    try:
+        polygon = array(polygon).astype(Float)   	   
+    except:
+        msg = 'Polygon %s could not be converted to Numeric array' %polygon
+	raise msg	 	
+    
+    assert len(polygon.shape) == 2,\
+       'Polygon array must be a 2d array of vertices: %s' %polygon
+
+	
+    assert polygon.shape[1] == 2,\
+       'Polygon array must have two columns: %s' %polygon    
+
+
+    if len(point.shape) == 1:
+        one_point = True
+	points = reshape(point, (1,2))
+    else:	
+        one_point = False
+	points = point
+
+    from util_ext import inside_polygon
+
+    indices = zeros( points.shape[0], Int )
+
+    count = inside_polygon(points, polygon, indices,
+                           int(closed), int(verbose))
+
+    #print 'O', point, count
+
+    if one_point:
+        return count == 1
+    else:
+        return indices[:count]
+
 class Polygon_function:
     """Create callable object f: x,y -> z, where a,y,z are vectors and 
@@ -1134,5 +1187,5 @@
     Examples:
        populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
-       will return five randomly sleected points inside the unit square 
+       will return five randomly selected points inside the unit square 
     """
 
@@ -1183,6 +1236,4 @@
 
 
-
-
 ##############################################
 #Initialise module
@@ -1190,5 +1241,6 @@
 import compile
 if compile.can_use_C_extension('util_ext.c'):
-    from util_ext import gradient, point_on_line 
+    from util_ext import gradient, point_on_line
+    inside_polygon = inside_polygon_c
 else:
     gradient = gradient_python
Index: inundation/ga/storm_surge/pyvolution/util_ext.c
===================================================================
--- inundation/ga/storm_surge/pyvolution/util_ext.c	(revision 998)
+++ inundation/ga/storm_surge/pyvolution/util_ext.c	(revision 999)
@@ -59,4 +59,86 @@
 
 
+int _inside_polygon(int M,           // Number of points
+		    int N,           // Number of polygon vertices
+		    double* points,
+		    double* polygon,
+		    int* indices,    // M-Array for storage indices 
+		    int closed, 
+		    int verbose) {
+  
+  double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j;
+  int i, j, k, count, inside;
+  
+  //Find min and max of poly used for optimisation when points 
+  //are far away from polygon 
+  
+  
+  minpx = polygon[0]; maxpx = minpx;
+  minpy = polygon[1]; maxpy = minpy;  
+  
+  for (i=0; i<N; i++) {
+    px_i = polygon[2*i];
+    py_i = polygon[2*i + 1];
+  
+    if (px_i < minpx) minpx = px_i;
+    if (px_i > maxpx) maxpx = px_i;    
+    if (py_i < minpy) minpy = py_i;
+    if (py_i > maxpy) maxpy = py_i;        
+  }
+  
+  //Begin main loop (for each point)
+  count = 0;
+  for (k=0; k<M; k++) {
+    //FIXME: Do this later
+    //if (verbose){
+    //  if (k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
+        
+    x = points[2*k];
+    y = points[2*k + 1];	
+
+    inside = 0;
+
+    //Optimisation
+    if ((x > maxpx) || (x < minpx)) continue;
+    if ((y > maxpy) || (y < minpy)) continue;        
+
+    //Check polygon 
+    for (i=0; i<N; i++) {
+      //printf("k,i=%d,%d\n", k, i);
+      j = (i+1)%N;
+      
+      px_i = polygon[2*i];
+      py_i = polygon[2*i+1];
+      px_j = polygon[2*j];
+      py_j = polygon[2*j+1];      
+      
+      //Check for case where point is contained in line segment	
+      if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) {
+	if (closed == 1) {
+	  inside = 1;
+	} else {
+	  inside = 0;
+	}
+	break;  
+      } else {	 
+        //Check if truly inside polygon		    
+	if ( ((py_i < y) && (py_j >= y)) ||
+	     ((py_j < y) && (py_i >= y)) ) {
+	  if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x)
+	    inside = 1-inside;
+	}
+      }
+    }
+    if (inside == 1) {
+      indices[count] = k;
+      count++;
+    }       
+  } // End k
+  
+  return count;
+}
+
+
+
 // Gateways to Python
 PyObject *point_on_line(PyObject *self, PyObject *args) {
@@ -81,4 +163,80 @@
   return result;
 }
+
+
+
+PyObject *inside_polygon(PyObject *self, PyObject *args) {
+  //def inside_polygon(point, polygon, closed, verbose, one_point):  
+  //  """Determine whether points are inside or outside a polygon
+  //  
+  //  Input:
+  //     point - Tuple of (x, y) coordinates, or list of tuples
+  //     polygon - list of vertices of polygon
+  //     one_poin - If True Boolean value should be returned 
+  //                If False, indices of points inside returned
+  //     closed - (optional) determine whether points on boundary should be 
+  //     regarded as belonging to the polygon (closed = True) 
+  //     or not (closed = False) 
+  
+  //  
+  //  Output:
+  //     If one point is considered, True or False is returned.
+  //     If multiple points are passed in, the function returns indices 
+  //     of those points that fall inside the polygon     
+  //     
+  //  Examples:
+  //     inside_polygon( [0.5, 0.5], [[0,0], [1,0], [1,1], [0,1]] )
+  //     will evaluate to True as the point 0.5, 0.5 is inside the unit square 
+  //     
+  //     inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
+  //     will return the indices [0, 2] as only the first and the last point 
+  //     is inside the unit square
+  //     
+  //  Remarks:
+  //     The vertices may be listed clockwise or counterclockwise and
+  //     the first point may optionally be repeated.    
+  //     Polygons do not need to be convex.
+  //     Polygons can have holes in them and points inside a hole is 
+  //     regarded as being outside the polygon. 	       
+  //	       
+  //     
+  //  Algorithm is based on work by Darel Finley, 
+  //  http://www.alienryderflex.com/polygon/ 
+  //
+  //
+  
+  PyArrayObject 
+    *point, 
+    *polygon, 
+    *indices;
+    
+  int closed, verbose; //Flags     
+  int count, M, N;
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "OOOii", 
+			&point, 
+			&polygon, 
+			&indices, 
+			&closed, 
+			&verbose))  
+    return NULL;
+
+  M = point -> dimensions[0];    //Number of points          
+  N = polygon -> dimensions[0];  //Number of vertices in polygon	
+
+  //printf("M=%d, N=%d\n", M, N);
+  // Call underlying routine
+  count = _inside_polygon(M, N, 
+			  (double*) point -> data,
+			  (double*) polygon -> data,
+			  (int*) indices -> data,
+			  closed, verbose); 
+			   
+  //NOTE: return number of points inside..
+  //printf("COunt=%d\n", count);
+  return Py_BuildValue("i", count);
+}
+
 
 
@@ -117,4 +275,5 @@
   {"gradient", gradient, METH_VARARGS, "Print out"},    
   {"point_on_line", point_on_line, METH_VARARGS, "Print out"},      
+  {"inside_polygon", inside_polygon, METH_VARARGS, "Print out"},      
   {NULL, NULL, 0, NULL}   /* sentinel */
 };
