Index: /anuga_core/source/anuga/shallow_water/test_data_manager.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 5530)
+++ /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 5531)
@@ -76,5 +76,5 @@
 
 
-        domain.distribute_to_vertices_and_edges()		
+        domain.distribute_to_vertices_and_edges()               
         self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
 
@@ -580,5 +580,5 @@
         for t in self.domain.evolve(yieldstep = 1.0, finaltime = 4.0):
             
-	    #########self.domain.write_time(track_speeds=True)
+            #########self.domain.write_time(track_speeds=True)
             stage = self.domain.quantities['stage'].vertex_values
 
@@ -6902,5 +6902,5 @@
 
     def test_file_boundary_stsIII(self):
-        """Read correct points from order file
+        """Read correct points from ordering file
         """
         from anuga.shallow_water import Domain
@@ -8490,9 +8490,9 @@
     def test_get_mesh_and_quantities_from_sww_file(self):
         """test_get_mesh_and_quantities_from_sww_file(self):
-	"""	
-	
-	# Generate a test sww file with non trivial georeference
-	
-	import time, os
+        """     
+        
+        # Generate a test sww file with non trivial georeference
+        
+        import time, os
         from Numeric import array, zeros, allclose, Float, concatenate
         from Scientific.IO.NetCDF import NetCDFFile
@@ -8524,5 +8524,5 @@
 
         
-	# Read it
+        # Read it
 
         # Get mesh and quantities from sww file
@@ -8541,10 +8541,10 @@
 
         # Check that time has been recovered
-	assert allclose(time, range(t_end+1))
+        assert allclose(time, range(t_end+1))
 
         # Check that quantities have been recovered
         # (sww files use single precision)
         z=domain.get_quantity('elevation').get_values(location='unique vertices')
-	assert allclose(quantities['elevation'], z)
+        assert allclose(quantities['elevation'], z)
 
         for q in ['stage', 'xmomentum', 'ymomentum']:
Index: /anuga_core/source/anuga/shallow_water/urs_ext.c
===================================================================
--- /anuga_core/source/anuga/shallow_water/urs_ext.c	(revision 5530)
+++ /anuga_core/source/anuga/shallow_water/urs_ext.c	(revision 5531)
@@ -116,5 +116,5 @@
         if (!fname) 
         {
-	  //printf("error\n");          
+          //printf("error\n");          
             PyErr_SetString(PyExc_ValueError, "filename not a string");
             return NULL;
@@ -256,5 +256,5 @@
 {
     FILE *fp;
-    int nsta, nsta0, i, isrc, ista;
+    int nsta, nsta0, i, isrc, ista, k;
     struct tgsrwg *mytgs=0, *mytgs0=0;
     char *muxFileName;                                                                  
@@ -299,7 +299,6 @@
     }
 
-    /* loop over remaining sources, check compatibility, and read them into 
-    *muxData */ 
-        for (isrc = 0; isrc < numSrc; isrc++)
+    /* loop over all sources, read headers and check compatibility */
+    for (isrc = 0; isrc < numSrc; isrc++)
     {
         muxFileName = muxFileNameArray[isrc];
@@ -376,5 +375,5 @@
 
     /* make array(s) to hold the demuxed data */
-    //FIXME: This can be reduced to only contain stations in order file
+    //FIXME: This can be reduced to only contain stations given in permutation file
     sts_data = (float**)malloc(nsta0*sizeof(float *));
     if (sts_data == NULL)
@@ -432,4 +431,13 @@
 
         /* loop over stations */
+        /* This is where we should only take stations with indices in the permutation array 
+           (and in that order) */
+        
+        //for i in range(len(permution)):
+        //     ista = permutation[i]
+        //     use ista with mux data
+        //     use i with the processed data to be returned         
+        i=0;
+        
         for(ista = 0; ista < nsta0; ista++)
         {               
@@ -438,14 +446,14 @@
 
             /* weight appropriately and add */
-            for(i = 0; i < mytgs0[ista].nt; i++)
+            for(k = 0; k < mytgs0[ista].nt; k++)
             {
-                if((isdata(sts_data[ista][i])) && isdata(temp_sts_data[i]))
+                if((isdata(sts_data[ista][k])) && isdata(temp_sts_data[k]))
                 {
-                    sts_data[ista][i] += temp_sts_data[i] * weights[isrc];
-                    //printf("%d,%d,%f\n",ista,i,sts_data[ista][i]);
+                    sts_data[ista][k] += temp_sts_data[k] * weights[isrc];
+                    //printf("%d,%d,%f\n",ista,k,sts_data[ista][k]);
                 }
                 else
                 {
-                    sts_data[ista][i] = NODATA;
+                    sts_data[ista][k] = NODATA;
                 }
             }
