Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py	(revision 7672)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py	(revision 7673)
@@ -180,5 +180,5 @@
                                      verbose=verbose,
                                      use_cache=use_cache,
-									 output_centroids = output_centroids)
+                                     output_centroids = output_centroids)
 
         if quake_offset_time is None:
@@ -190,6 +190,6 @@
                 quake_time = time + quake_offset_time
                 points_list = [quake_time, quake_time/3600.]# fudge around SWW time bug
-                point_quantities = callable_sww(time,point_i)
-                
+                point_quantities = callable_sww(time,point_i) # __call__ is overridden
+				
                 for quantity in quantities:
                     if quantity == NAN:
@@ -271,5 +271,6 @@
                    title_on=None,
                    use_cache=False,
-                   verbose=False):
+                   verbose=False,
+				   output_centroids=False):
     """ Read sww file and plot the time series for the
     prescribed quantities at defined gauge locations and
@@ -388,5 +389,6 @@
                         title_on,
                         use_cache,
-                        verbose)
+                        verbose,
+						output_centroids = output_centroids)
     return k
 
@@ -423,5 +425,6 @@
                     title_on = None,
                     use_cache = False,
-                    verbose = False):   
+                    verbose = False,
+					output_centroids = False):   
         
     # FIXME(Ole): Shouldn't print statements here be governed by verbose?
@@ -486,5 +489,6 @@
                           time_thinning = time_thinning,
                           verbose = verbose,
-                          use_cache = use_cache)
+                          use_cache = use_cache,
+						  output_centroids = output_centroids)
 
         # determine which gauges are contained in sww file
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py	(revision 7672)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py	(revision 7673)
@@ -451,17 +451,17 @@
 
         # create a csv file containing our gauge points
-        points = [[2.0,1.0],[5.0,5.0]]
+        points = [[2.0,1.0],[4.5,4.0]]
 
         points_file = tempfile.mktemp(".csv")
         file_id = open(points_file,"w")
-# These values are where the centroids should be
+# These values are where the centroids should be		
 #        file_id.write("name, easting, northing, elevation \n\
-#point1, 1.5, 1.5, 3.0\n\
-#point2, 4.5, 4.5, 9.0\n")
+#point1, 2.0, 2.0, 3.0\n\
+#point2, 4.0, 4.0, 9.0\n")
  
- # These values are slightly off the centroids - will it find them?
+# These values are slightly off the centroids - will it find the centroids?
         file_id.write("name, easting, northing, elevation \n\
 point1, 2.0, 1.0, 3.0\n\
-point2, 5.5, 4.0, 9.0\n")
+point2, 4.5, 4.0, 9.0\n")
 
  
@@ -474,5 +474,5 @@
                        output_centroids=True)
 
-        point1_answers_array = [[0.0,0.0,1.0,2.5,-1.5,3.0,4.0], [2.0,2.0/3600.,10.0,11.5,-1.5,3.0,4.0]]
+        point1_answers_array = [[0.0,0.0,1.0,3.0,-2.0,3.0,4.0], [2.0,2.0/3600.,10.0,12.0,-2.0,3.0,4.0]]
         point1_filename = 'gauge_point1.csv'
         point1_handle = file(point1_filename)
@@ -482,5 +482,4 @@
         line=[]
         for i,row in enumerate(point1_reader):
-#           print 'i',i,'row',row
             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
                          float(row[4]),float(row[5]),float(row[6])])
@@ -488,5 +487,5 @@
             assert num.allclose(line[i], point1_answers_array[i])
 
-        point2_answers_array = [[0.0,0.0,1.0,5.5,-4.5,3.0,4.0], [2.0,2.0/3600.,10.0,14.5,-4.5,3.0,4.0]]
+        point2_answers_array = [[0.0,0.0,1.0,5.0,-4.0,3.0,4.0], [2.0,2.0/3600.,10.0,14.0,-4.0,3.0,4.0]]
         point2_filename = 'gauge_point2.csv' 
         point2_handle = file(point2_filename)
@@ -496,8 +495,7 @@
         line=[]
         for i,row in enumerate(point2_reader):
-#           print 'i',i,'row',row
             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
                          float(row[4]),float(row[5]),float(row[6])])
-#           print 'assert line',line[i],'point1',point1_answers_array[i]
+#           print i, 'assert line',line[i],'point2',point2_answers_array[i]
             assert num.allclose(line[i], point2_answers_array[i])
                          
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 7672)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 7673)
@@ -54,5 +54,5 @@
                   use_cache=False,
                   boundary_polygon=None,
-				  output_centroids=False):
+                  output_centroids=False):
     """Read time history of spatial data from NetCDF file and return
     a callable object.
@@ -140,5 +140,5 @@
               'verbose': verbose,
               'boundary_polygon': boundary_polygon,
-			  'output_centroids': output_centroids}
+              'output_centroids': output_centroids}
 
     # Call underlying engine with or without caching
@@ -232,5 +232,5 @@
                                         verbose=verbose,
                                         boundary_polygon=boundary_polygon,
-										output_centroids=output_centroids)
+                                        output_centroids=output_centroids)
     else:
         # FIXME (Ole): Could add csv file here to address Ted Rigby's
@@ -260,5 +260,5 @@
                              verbose=False,
                              boundary_polygon=None,
-							 output_centroids=False):
+                             output_centroids=False):
     """Read time history of spatial data from NetCDF sww file and
     return a callable object f(t,x,y)
@@ -501,5 +501,6 @@
                                    time_thinning=time_thinning,
                                    verbose=verbose,
-                                   gauge_neighbour_id=gauge_neighbour_id),
+                                   gauge_neighbour_id=gauge_neighbour_id,
+								   output_centroids=output_centroids),
             starttime)
 
Index: /anuga_core/source/anuga/fit_interpolate/interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 7672)
+++ /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 7673)
@@ -69,5 +69,6 @@
                 start_blocking_len=500000,
                 use_cache=False,
-                verbose=False):
+                verbose=False,
+				output_centroids=False):
     """Interpolate vertex_values to interpolation points.
 
@@ -153,5 +154,6 @@
     result = I.interpolate_block(vertex_values, interpolation_points,
                                  use_cache=use_cache,
-                                 verbose=verbose)
+                                 verbose=verbose,
+								 output_centroids=output_centroids)
 
     return result
@@ -228,5 +230,6 @@
                     point_coordinates=None,
                     start_blocking_len=500000,
-                    verbose=False):
+                    verbose=False,
+					output_centroids=False):
         """Interpolate mesh data f to determine values, z, at points.
 
@@ -288,5 +291,5 @@
                 self._A_can_be_reused = True
                 z = self.interpolate_block(f, point_coordinates,
-                                           verbose=verbose)
+                                           verbose=verbose, output_centroids=output_centroids)
             else:
                 # Handle blocking
@@ -305,5 +308,5 @@
                                  start_blocking_len):
                     t = self.interpolate_block(f, point_coordinates[start:end],
-                                               verbose=verbose)
+                                               verbose=verbose, output_centroids=output_centroids)
                     z = num.concatenate((z, t), axis=0)    #??default#
                     start = end
@@ -311,5 +314,5 @@
                 end = len(point_coordinates)
                 t = self.interpolate_block(f, point_coordinates[start:end],
-                                           verbose=verbose)
+                                           verbose=verbose, output_centroids=output_centroids)
                 z = num.concatenate((z, t), axis=0)    #??default#
         return z
@@ -317,12 +320,12 @@
 
     ##
-    # @brief Control whether blocking occurs or not.
-    # @param f ??
-    # @param point_coordinates ??
-    # @param use_cache ??
+    # @brief Interpolate a block of vertices
+    # @param f Array of arbitrary data to be interpolated
+    # @param point_coordinates List of vertices to intersect with the mesh
+    # @param use_cache True if caching should be used to accelerate the calculations
     # @param verbose True if this function is verbose.
     # @return ??
     def interpolate_block(self, f, point_coordinates,
-                          use_cache=False, verbose=False):
+                          use_cache=False, verbose=False, output_centroids=False):
         """
         Call this if you want to control the blocking or make sure blocking
@@ -354,5 +357,5 @@
 
                 X = cache(self._build_interpolation_matrix_A,
-                          args=(point_coordinates,),
+                          args=(point_coordinates, output_centroids),
                           kwargs={'verbose': verbose},
                           verbose=verbose)
@@ -372,8 +375,9 @@
                 if reuse_A is False:
                     X = self._build_interpolation_matrix_A(point_coordinates,
+                                                           output_centroids,
                                                            verbose=verbose)
                     self.interpolation_matrices[key] = (X, point_coordinates)
         else:
-            X = self._build_interpolation_matrix_A(point_coordinates,
+            X = self._build_interpolation_matrix_A(point_coordinates, output_centroids,
                                                    verbose=verbose)
 
@@ -403,8 +407,11 @@
 
     ##
-    # @brief ??
-    # @param f ??
+    # @brief Get interpolated data at given points.
+	#        Applies a transform to all points to calculate the
+	#        interpolated values. Points outside the mesh are returned as NaN.
+	# @note self._A matrix must be valid
+    # @param f Array of arbitrary data
     # @param verbose True if this function is to be verbose.
-    # @return ??
+    # @return f transformed by interpolation matrix (f')
     def _get_point_data_z(self, f, verbose=False):
         """
@@ -424,9 +431,13 @@
     ##
     # @brief Build NxM interpolation matrix.
-    # @param point_coordinates ??
+    # @param point_coordinates Points to sample at
+	# @param output_centroids set to True to always sample from the centre
+	#                         of the intersected triangle, instead of the intersection
+	#                         point.
     # @param verbose True if this function is to be verbose.
-    # @return ??
+    # @return Interpolation matrix A, plus lists of the points inside and outside the mesh.
     def _build_interpolation_matrix_A(self,
                                       point_coordinates,
+                                      output_centroids=False,
                                       verbose=False):
         """Build n x m interpolation matrix, where
@@ -495,10 +506,15 @@
                 j1 = self.mesh.triangles[k,1] # Global vertex id for sigma1
                 j2 = self.mesh.triangles[k,2] # Global vertex id for sigma2
-
-                sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
-                js     = [j0, j1, j2]
-
-                for j in js:
-                    A[i, j] = sigmas[j]
+                js = [j0, j1, j2]
+				
+                if output_centroids is False:
+                    # Weight each vertex according to its distance from x
+                    sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
+                    for j in js:
+                        A[i, j] = sigmas[j]
+                else:
+				    # If centroids are needed, weight all 3 vertices equally
+                    for j in js:
+                        A[i, j] = 1.0/3.0                    
             else:
                 msg = 'Could not find triangle for point', x
@@ -751,5 +767,6 @@
                  time_thinning=1,
                  verbose=False,
-                 gauge_neighbour_id=None):
+                 gauge_neighbour_id=None,
+				 output_centroids=False):
         """Initialise object and build spatial interpolation if required
 
@@ -976,5 +993,6 @@
                                                       point_coordinates=\
                                                       self.interpolation_points,
-                                                      verbose=False) # No clutter
+                                                      verbose=False,
+													  output_centroids=output_centroids)
                     elif triangles is None and vertex_coordinates is not None:
                         result = interpolate_polyline(Q,
@@ -1002,5 +1020,5 @@
 
     ##
-    # @brief Override object() method.
+    # @brief Evaluate interpolation function
     # @param t Model time - must lie within existing timesteps.
     # @param point_id Index of one of the preprocessed points.
Index: /anuga_core/source/anuga/fit_interpolate/test_interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/test_interpolate.py	(revision 7672)
+++ /anuga_core/source/anuga/fit_interpolate/test_interpolate.py	(revision 7673)
@@ -391,4 +391,28 @@
         assert num.allclose(num.sum(results, axis=1), 1.0)
 
+		
+    def test_arbitrary_datapoints_return_centroids(self):
+        #Try arbitrary datapoints, making sure they return
+        #an interpolation matrix for the intersected triangle's
+        #centroid.
+        
+        a = [1.0, 0.0]
+        b = [0.0, 3.0]
+        c = [4.0,0.0]
+        points = [a, b, c]
+        vertices = [ [1,0,2] ]
+
+        data = [ [1.2, 1.5], [1.123, 1.768], [2.43, 0.44] ]
+
+        interp = Interpolate(points, vertices)
+        
+        third = [1.0/3.0, 1.0/3.0, 1.0/3.0]
+        answer = [third, third, third]
+		
+        A,_,_ = interp._build_interpolation_matrix_A(data, output_centroids=True)
+        results = A.todense()
+        assert num.allclose(results, answer)		
+		
+		
     def test_arbitrary_datapoints_some_outside(self):
         #Try arbitrary datapoints one outside the triangle.
