Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py	(revision 7674)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py	(revision 7675)
@@ -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) # __call__ is overridden
-				
+                point_quantities = callable_sww(time, point_i) # __call__ is overridden
+                            
                 for quantity in quantities:
                     if quantity == NAN:
@@ -197,4 +197,5 @@
                                      % callable_sww.get_name)
                     else:
+                        #core quantities that are exported from the interpolator     
                         if quantity == 'stage':
                             points_list.append(point_quantities[0])
@@ -208,5 +209,6 @@
                         if quantity == 'ymomentum':
                             points_list.append(point_quantities[3])
-                            
+
+                        #derived quantities that are calculated from the core ones
                         if quantity == 'depth':
                             points_list.append(point_quantities[0] 
@@ -237,5 +239,10 @@
                             points_list.append(calc_bearing(point_quantities[2],
                                                             point_quantities[3]))
-
+                        if quantity == 'xcentroid':
+                            points_list.append(callable_sww.centroids[point_i][0])
+
+                        if quantity == 'ycentroid':
+                            points_list.append(callable_sww.centroids[point_i][1])
+							
                 points_writer[point_i].writerow(points_list)
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py	(revision 7674)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py	(revision 7675)
@@ -451,6 +451,4 @@
 
         # create a csv file containing our gauge points
-        points = [[2.0,1.0],[4.5,4.0]]
-
         points_file = tempfile.mktemp(".csv")
         file_id = open(points_file,"w")
@@ -510,4 +508,98 @@
 
 
+    def test_sww2csv_output_centroid_attribute(self):
+
+        def elevation_function(x, y):
+            return -x
+        
+        """Check sww2csv timeseries at centroid, then output the centroid coordinates.
+        
+        Test the ability to get a timeseries at the centroid of a triangle, rather
+		than the given gauge point, then output the results.
+        """
+        
+        # Create rectangular mesh
+        mesh_file = tempfile.mktemp(".tsh")    
+        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
+        m = Mesh()
+        m.add_vertices(points)
+        m.auto_segment()
+        m.generate_mesh(verbose=False)
+        m.export_mesh_file(mesh_file)
+        
+        # Create shallow water domain
+        domain = Domain(mesh_file)
+        
+        domain.default_order=2
+        
+        # This test was made before tight_slope_limiters were introduced
+        # Since were are testing interpolation values this is OK
+        domain.tight_slope_limiters = 0 
+        
+
+        # Set some field values
+        domain.set_quantity('elevation', elevation_function)
+        domain.set_quantity('friction', 0.03)
+
+        ######################
+        # Boundary conditions
+        B = Transmissive_boundary(domain)
+        domain.set_boundary( {'exterior': B})
+
+        # This call mangles the stage values.
+        domain.distribute_to_vertices_and_edges()
+        domain.set_quantity('stage', 1.0)
+
+
+        domain.set_name('datatest' + str(time.time()))
+        domain.smooth = True
+        domain.reduction = mean
+
+
+        sww = SWW_file(domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+        domain.set_quantity('stage', 10.0) # This is automatically limited
+        # so it will not be less than the elevation
+        domain.time = 2.
+        sww.store_timestep()
+
+        # create a csv file containing our gauge points
+        points_file = tempfile.mktemp(".csv")
+        file_id = open(points_file,"w")
+ 
+# These values are slightly off the centroids - will it find the centroids?
+        file_id.write("name, easting, northing, elevation \n\
+point1, 2.5, 4.25, 3.0\n")
+
+        file_id.close()
+
+        gauge_sww2csv(sww.filename, 
+                       points_file,
+                       quantities=['stage', 'xcentroid', 'ycentroid'],
+                       verbose=False,
+                       use_cache=False,
+                       output_centroids=True)
+
+        point1_answers_array = [[0.0,0.0,1.0,4.0,4.0], [2.0,2.0/3600.,10.0,4.0,4.0]]
+        point1_filename = 'gauge_point1.csv'
+        point1_handle = file(point1_filename)
+        point1_reader = reader(point1_handle)
+        point1_reader.next()
+
+        line=[]
+        for i,row in enumerate(point1_reader):
+            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4])])
+#            print 'assert line',line[i],'point1',point1_answers_array[i]
+            assert num.allclose(line[i], point1_answers_array[i])
+
+        # clean up
+        point1_handle.close()
+        os.remove(mesh_file)		
+        os.remove(sww.filename)
+        os.remove(points_file)
+        os.remove(point1_filename)
+		
+		
 
 #-------------------------------------------------------------
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 7674)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 7675)
@@ -502,5 +502,5 @@
                                    verbose=verbose,
                                    gauge_neighbour_id=gauge_neighbour_id,
-								   output_centroids=output_centroids),
+                                   output_centroids=output_centroids),
             starttime)
 
Index: /anuga_core/source/anuga/fit_interpolate/interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 7674)
+++ /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 7675)
@@ -218,5 +218,4 @@
 
 
-
     ##
     # @brief Interpolate mesh data f to determine values, z, at points.
@@ -383,5 +382,5 @@
 
         # Unpack result
-        self._A, self.inside_poly_indices, self.outside_poly_indices = X
+        self._A, self.inside_poly_indices, self.outside_poly_indices, self.centroids = X
 
         # Check that input dimensions are compatible
@@ -436,5 +435,6 @@
 	#                         point.
     # @param verbose True if this function is to be verbose.
-    # @return Interpolation matrix A, plus lists of the points inside and outside the mesh.
+    # @return Interpolation matrix A, plus lists of the points inside and outside the mesh
+	#         and the list of centroids, if requested.
     def _build_interpolation_matrix_A(self,
                                       point_coordinates,
@@ -487,4 +487,6 @@
         n = len(inside_poly_indices)
 
+        centroids = []
+		
         # Compute matrix elements for points inside the mesh
         if verbose: log.critical('Building interpolation matrix from %d points'
@@ -516,9 +518,10 @@
 				    # If centroids are needed, weight all 3 vertices equally
                     for j in js:
-                        A[i, j] = 1.0/3.0                    
+                        A[i, j] = 1.0/3.0
+                    centroids.append(self.mesh.centroid_coordinates[k])						
             else:
                 msg = 'Could not find triangle for point', x
                 raise Exception(msg)
-        return A, inside_poly_indices, outside_poly_indices
+        return A, inside_poly_indices, outside_poly_indices, centroids
 
 
@@ -841,4 +844,5 @@
         self.index = 0    # Initial time index
         self.precomputed_values = {}
+        self.centroids = []
 
         # Precomputed spatial interpolation if requested
@@ -994,5 +998,6 @@
                                                       self.interpolation_points,
                                                       verbose=False,
-													  output_centroids=output_centroids)
+                                                      output_centroids=output_centroids)
+                        self.centroids = interpol.centroids														  
                     elif triangles is None and vertex_coordinates is not None:
                         result = interpolate_polyline(Q,
@@ -1003,9 +1008,9 @@
 
                     #assert len(result), len(interpolation_points)
-                    self.precomputed_values[name][i, :] = result
-
+                    self.precomputed_values[name][i, :] = result									
+					
             # Report
             if verbose:
-                log.critical(self.statistics())
+                log.critical(self.statistics())			
         else:
             # Store quantitites as is
Index: /anuga_core/source/anuga/fit_interpolate/test_interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/test_interpolate.py	(revision 7674)
+++ /anuga_core/source/anuga/fit_interpolate/test_interpolate.py	(revision 7675)
@@ -104,5 +104,5 @@
 
         interp = Interpolate(points, vertices)
-        A, _, _ = interp._build_interpolation_matrix_A(data)
+        A, _, _, _ = interp._build_interpolation_matrix_A(data)
         assert num.allclose(A.todense(), [[1./3, 1./3, 1./3]])
 
@@ -285,5 +285,5 @@
                       0., 0. , 0., 0., 0., 0.]]
 
-        A,_,_ = interp._build_interpolation_matrix_A(data)
+        A,_,_,_ = interp._build_interpolation_matrix_A(data)
         assert num.allclose(A.todense(), answer)
         
@@ -294,5 +294,5 @@
                       0., 0. , 0., 0., 0., 0.]]
         
-        A,_,_ = interp._build_interpolation_matrix_A(data)        
+        A,_,_,_ = interp._build_interpolation_matrix_A(data)        
         assert num.allclose(A.todense(), answer)
 
@@ -305,5 +305,5 @@
                       0., 0. , 0., 0., 0., 0.]]
                       
-        A,_,_ = interp._build_interpolation_matrix_A(data)        
+        A,_,_,_ = interp._build_interpolation_matrix_A(data)        
         assert num.allclose(A.todense(), answer)
 
@@ -326,5 +326,5 @@
                    [0., 0., 1.]]
                    
-        A,_,_ = interp._build_interpolation_matrix_A(data)
+        A,_,_,_ = interp._build_interpolation_matrix_A(data)
         assert num.allclose(A.todense(), answer)
 
@@ -347,5 +347,5 @@
         interp = Interpolate(points, vertices)
 
-        A,_,_ = interp._build_interpolation_matrix_A(data)
+        A,_,_,_ = interp._build_interpolation_matrix_A(data)
         assert num.allclose(A.todense(), answer)
 
@@ -368,5 +368,5 @@
         interp = Interpolate(points, vertices)
 
-        A,_,_ = interp._build_interpolation_matrix_A(data)
+        A,_,_,_ = interp._build_interpolation_matrix_A(data)
         assert num.allclose(A.todense(), answer)
 
@@ -387,5 +387,5 @@
         #print "interp.get_A()", interp.get_A()
         
-        A,_,_ = interp._build_interpolation_matrix_A(data)
+        A,_,_,_ = interp._build_interpolation_matrix_A(data)
         results = A.todense()
         assert num.allclose(num.sum(results, axis=1), 1.0)
@@ -410,5 +410,5 @@
         answer = [third, third, third]
 		
-        A,_,_ = interp._build_interpolation_matrix_A(data, output_centroids=True)
+        A,_,_,_ = interp._build_interpolation_matrix_A(data, output_centroids=True)
         results = A.todense()
         assert num.allclose(results, answer)		
@@ -430,5 +430,5 @@
         interp = Interpolate(points, vertices)
         
-        A,_,_ = interp._build_interpolation_matrix_A(data)
+        A,_,_,_ = interp._build_interpolation_matrix_A(data)
         results = A.todense()
         assert num.allclose(num.sum(results, axis=1), [1,1,1,0])
@@ -461,5 +461,5 @@
 
 
-        A,_,_ = interp._build_interpolation_matrix_A(data)
+        A,_,_,_ = interp._build_interpolation_matrix_A(data)
         A = A.todense()
         for i in range(A.shape[0]):
