Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 3849)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 3850)
@@ -641,5 +641,5 @@
         """
         if name.endswith('.sww'):
-            name = name[-4:]
+            name = name[:-4]
             
         self.simulation_name = name
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 3849)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 3850)
@@ -235,5 +235,5 @@
 
     #Get variables
-    if verbose: print 'Get variables'    
+    #if verbose: print 'Get variables'    
     time = fid.variables['time'][:]    
 
@@ -318,5 +318,5 @@
                                   triangles,
                                   interpolation_points,
-                                  verbose = verbose)
+                                  verbose=verbose)
 
 
Index: /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py	(revision 3849)
+++ /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py	(revision 3850)
@@ -73,5 +73,5 @@
               a mesh origin, since geospatial has its own mesh origin.
         """
-        
+
         #Convert input to Numeric arrays
         triangles = ensure_numeric(triangles, Int)
@@ -79,7 +79,11 @@
                                              geo_reference = mesh_origin)
 
-        #Don't pass geo_reference to mesh.  It doesn't work.
+        #Don't pass geo_reference to mesh.  It doesn't work. (Still??)
+        
+        if verbose: print 'FitInterpolate: Building mesh'        
         self.mesh = Mesh(vertex_coordinates, triangles)
         self.mesh.check_integrity()
+        
+        if verbose: print 'FitInterpolate: Building quad tree'
         self.root = build_quadtree(self.mesh,
                                    max_points_per_cell = max_vertices_per_cell)
Index: /anuga_core/source/anuga/fit_interpolate/interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 3849)
+++ /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 3850)
@@ -405,9 +405,9 @@
                  time,
                  quantities,
-                 quantity_names = None,  
-                 vertex_coordinates = None,
-                 triangles = None,
-                 interpolation_points = None,
-                 verbose = False):
+                 quantity_names=None,  
+                 vertex_coordinates=None,
+                 triangles=None,
+                 interpolation_points=None,
+                 verbose=False):
         """Initialise object and build spatial interpolation if required
         """
@@ -417,10 +417,9 @@
 
 
-        #from anuga.abstract_2d_finite_volumes.util import mean, ensure_numeric
         from anuga.config import time_format
         import types
 
 
-	#Check temporal info
+	# Check temporal info
         time = ensure_numeric(time)        
         msg = 'Time must be a monotonuosly '
@@ -429,5 +428,5 @@
 
 
-        #Check if quantities is a single array only
+        # Check if quantities is a single array only
         if type(quantities) != types.DictType:
             quantities = ensure_numeric(quantities)
@@ -438,10 +437,10 @@
 
 
-        #Use keys if no names are specified
+        # Use keys if no names are specified
         if quantity_names is None:
             quantity_names = quantities.keys()
 
 
-        #Check spatial info
+        # Check spatial info
         if vertex_coordinates is None:
             self.spatial = False
@@ -454,6 +453,5 @@
 
              
-        #Save for use with statistics
-
+        # Save for use with statistics
         self.quantities_range = {}
         for name in quantity_names:
@@ -462,5 +460,4 @@
         
         self.quantity_names = quantity_names        
-        #self.quantities = quantities  #Takes too much memory      
         self.vertex_coordinates = vertex_coordinates 
         self.interpolation_points = interpolation_points
@@ -470,5 +467,5 @@
         
             
-        #Precomputed spatial interpolation if requested
+        # Precomputed spatial interpolation if requested
         if interpolation_points is not None:
             if self.spatial is False:
@@ -491,5 +488,6 @@
                 self.precomputed_values[name] = zeros((p, m), Float)
 
-            #Build interpolator
+            # Build interpolator
+            if verbose: print 'Building interpolation matrix'            
             interpol = Interpolate(vertex_coordinates,
                                    triangles,
@@ -497,5 +495,5 @@
                                    #self.interpolation_points,
                                    #alpha = 0,
-                                   verbose = verbose)
+                                   verbose=verbose)
 
             if verbose: print 'Interpolate'
Index: /anuga_core/source/anuga/fit_interpolate/search_functions.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/search_functions.py	(revision 3849)
+++ /anuga_core/source/anuga/fit_interpolate/search_functions.py	(revision 3850)
@@ -8,4 +8,5 @@
 from Numeric import dot
 
+from numerical_tools import get_machine_precision
 
 def search_tree_of_vertices(root, mesh, x):
@@ -83,15 +84,16 @@
             sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
 
-            #FIXME: Maybe move out to test or something
-            epsilon = 1.0e-6
-            assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
+            # Integrity check - machine precision is too hard, but at
+            # least check that we are within a couple of orders of magnitude
+            epsilon = get_machine_precision()
+            msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, 100*eps = %.15e'\
+                  %(abs(sigma0+sigma1+sigma2-1), 100*epsilon)
+            assert abs(sigma0 + sigma1 + sigma2 - 1.0) < 100*epsilon, msg
 
             #Check that this triangle contains the data point
             
-            #Sigmas can get negative within
-            #machine precision on some machines (e.g nautilus)
-            #Hence the small eps
-            eps = 1.0e-15
-            if sigma0 >= -eps and sigma1 >= -eps and sigma2 >= -eps:
+            # Sigmas are allowed to get negative within
+            # machine precision on some machines (e.g nautilus)
+            if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
                 element_found = True
                 break
Index: /anuga_core/source/anuga/load_mesh/loadASCII.py
===================================================================
--- /anuga_core/source/anuga/load_mesh/loadASCII.py	(revision 3849)
+++ /anuga_core/source/anuga/load_mesh/loadASCII.py	(revision 3850)
@@ -76,4 +76,7 @@
 #  Needs to be defined
 ##FIXME (DSG-DSG) if the ascii file isn't .xya give a better error message.
+
+# FIXME (Ole): Has this stuff been superseded by geospatial data?
+# Much of this code is also there
 
 
Index: /anuga_core/source/anuga/utilities/numerical_tools.py
===================================================================
--- /anuga_core/source/anuga/utilities/numerical_tools.py	(revision 3849)
+++ /anuga_core/source/anuga/utilities/numerical_tools.py	(revision 3850)
@@ -10,9 +10,13 @@
 #(this should move to somewhere central)
 #try:
-#    from scipy import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float, arange    
+#    from scipy import ArrayType, array, sum, innerproduct, ravel, sqrt,
+# searchsorted, sort, concatenate, Float, arange    
 #except:
 #    #print 'Could not find scipy - using Numeric'
-#    from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float, arange
-from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float, arange    
+#    from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt,
+#searchsorted, sort, concatenate, Float, arange
+
+from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt,\
+     searchsorted, sort, concatenate, Float, arange    
 
 # Getting an infinite number to use when using Numeric
@@ -31,12 +35,9 @@
     """
 
-    error_msg = 'Input to acos is outside allowed domain [-1.0, 1.0]. I got %.12f' %x
+    error_msg = 'Input to acos is outside allowed domain [-1.0, 1.0].'+\
+                'I got %.12f' %x
     warning_msg = 'Changing argument to acos from %.18f to %.1f' %(x, sign(x))
 
-    # FIXME(Ole): Need function to compute machine precision as this is also
-    # used in fit_interpolate/search_functions.py
-    
-    
-    eps = 1.0e-15  # Machine precision 
+    eps = get_machine_precision() # Machine precision 
     if x < -1.0:
         if x < -1.0 - eps:
@@ -50,5 +51,5 @@
             raise ValueError, errmsg
         else:
-            print 'NOTE: changing argument to acos from %.18f to 1.0' %x            
+            print 'NOTE: changing argument to acos from %.18f to 1.0' %x
             x = 1.0
 
@@ -148,7 +149,9 @@
         y = x 
 
+    x = ensure_numeric(x)
+    y = ensure_numeric(y)       
     assert(len(x)==len(y))
-    N = len(x)
- 
+
+    N = len(x) 
     cx = x - mean(x)  
     cy = y - mean(y)  
Index: /anuga_validation/okushiri_2005/create_okushiri.py
===================================================================
--- /anuga_validation/okushiri_2005/create_okushiri.py	(revision 3849)
+++ /anuga_validation/okushiri_2005/create_okushiri.py	(revision 3850)
@@ -16,8 +16,10 @@
     """
 
-    print 'Preparing time boundary from %s' %filename
+    textversion = filename[:-4] + '.txt'
+
+    print 'Preparing time boundary from %s' %textversion
     from Numeric import array
 
-    fid = open(filename)
+    fid = open(textversion)
 
     #Skip first line
@@ -43,7 +45,6 @@
     from Scientific.IO.NetCDF import NetCDFFile
 
-    outfile = filename[:-4] + '.tms'
-    print 'Writing to', outfile
-    fid = NetCDFFile(outfile, 'w')
+    print 'Writing to', filename
+    fid = NetCDFFile(filename, 'w')
 
     fid.institution = 'Geoscience Australia'
@@ -69,6 +70,9 @@
 if __name__ == "__main__":
 
+
     #Preparing points
     #(Only when using original Benchmark_2_Bathymetry.txt file)
+    # FIXME: Replace by using geospatial data
+    #
     #from anuga.abstract_2d_finite_volumes.data_manager import xya2pts
     #xya2pts(project.bathymetry_filename, verbose = True,
Index: /anuga_validation/okushiri_2005/extract_timeseries.py
===================================================================
--- /anuga_validation/okushiri_2005/extract_timeseries.py	(revision 3849)
+++ /anuga_validation/okushiri_2005/extract_timeseries.py	(revision 3850)
@@ -5,128 +5,132 @@
 """
 
+from Numeric import allclose
+from Scientific.IO.NetCDF import NetCDFFile
+
+from anuga.abstract_2d_finite_volumes.util import file_function
+from anuga.utilities.numerical_tools import\
+     ensure_numeric, cov, get_machine_precision
+
+import project
+
+try:
+    from pylab import ion, hold, plot, title, legend, xlabel, ylabel, savefig
+except:
+    plotting = False
+else:
+    plotting = True
 
 
-from anuga.caching import cache
-
-
-
-
-gauges = [[0.000, 1.696]]  #Boundary
-gauges += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
-
-gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
-
+#-------------------------
+# Basic data
+#-------------------------
 
 finaltime = 22.5
 timestep = 0.05
 
-#Read reference data
+gauge_locations = [[0.000, 1.696]] # Boundary gauge
+gauge_locations += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
+gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
 
-#Input wave
-filename = 'Benchmark_2_input.tms'
-print 'Reading', filename
-from Scientific.IO.NetCDF import NetCDFFile
-from Numeric import allclose
-fid = NetCDFFile(filename, 'r')#
+validation_data = {}
+for key in gauge_names:
+    validation_data[key] = []
 
-input_time = fid.variables['time'][:]
-input_stage = fid.variables['stage'][:]
+    
+expected_covariances = {'Boundary': 5.288392008865989e-05,
+                        'ch5': 1.166748190444681e-04,
+                        'ch7': 1.121816242516758e-04,
+                        'ch9': 1.249543278366778e-04}
 
 
-#gauges
+#-------------------------
+# Read validation dataa
+#-------------------------
+
+print 'Reading', project.boundary_filename
+fid = NetCDFFile(project.boundary_filename, 'r')
+input_time = fid.variables['time'][:]
+validation_data['Boundary'] = fid.variables['stage'][:]
+
 reference_time = []
-ch5 = []
-ch7 = []
-ch9 = []
-filename = 'output_ch5-7-9.txt'
-fid = open(filename)
+fid = open(project.validation_filename)
 lines = fid.readlines()
 fid.close()
+
 for i, line in enumerate(lines[1:]):
     if i == len(input_time): break
-
+    
     fields = line.split()
 
-    reference_time.append(float(fields[0]))
-    ch5.append(float(fields[1])/100)   #cm2m
-    ch7.append(float(fields[2])/100)   #cm2m
-    ch9.append(float(fields[3])/100)   #cm2m
+    reference_time.append(float(fields[0]))    # Record reference time
+    for j, key in enumerate(gauge_names[1:]):  # Omit boundary gauge
+        value = float(fields[1:][j])           # Omit time 
+        validation_data[key].append(value/100) # Convert cm2m
 
 
-from anuga.abstract_2d_finite_volumes.util import file_function
-from anuga.utilities.numerical_tools import ensure_numeric
-gauge_values = [ensure_numeric(input_stage),
-                ensure_numeric(ch5),
-                ensure_numeric(ch7),
-                ensure_numeric(ch9)] #Reference values
-
-
-
-#Read model output
-#filename = 'output.sww'
-filename = 'lwru2.sww'
-#filename = 'lwru2_.05s.sww'
-
-#f = file_function(filename,
-#                  quantities = 'stage',
-#                  interpolation_points = gauges,
-#                  verbose = True)
-
-f = cache(file_function,filename,
-          {'quantities': 'stage',
-           'interpolation_points': gauges,
-           'verbose': True},
-          #evaluate = True,
-          dependencies = [filename],
-          verbose = True)
-
-
-
-
-#Checks
-#print reference_time
-#print input_time
+# Checks
 assert reference_time[0] == 0.0
 assert reference_time[-1] == finaltime
 assert allclose(reference_time, input_time)
 
+for key in gauge_names:
+    validation_data[key] = ensure_numeric(validation_data[key])
 
 
-#Validation
+#--------------------------------------------------
+# Read and interpolate model output
+#--------------------------------------------------
+#f = file_function('okushiri_new_limiters.sww',
+#f = file_function('okushiri_as2005_with_mxspd=0.1.sww',
+f = file_function(project.output_filename,
+                  quantities='stage',
+                  interpolation_points=gauge_locations,
+                  use_cache=True,
+                  verbose=True)
 
 
+#--------------------------------------------------
+# Compare model output to validation data
+#--------------------------------------------------
+
+eps = get_machine_precision()
 for k, name in enumerate(gauge_names):
     sqsum = 0
     denom = 0
     model = []
+    print 
     print 'Validating ' + name
+    observed_timeseries = validation_data[name]
     for i, t in enumerate(reference_time):
-        ref = gauge_values[k][i]
-        val = f(t, point_id = k)[0]
-        model.append(val)
+        model.append(f(t, point_id=k)[0])
 
-        sqsum += (ref - val)**2
-        denom += ref**2
+    res = cov(observed_timeseries, model)     
+    print 'Covariance = %.18e' %res
 
-    print sqsum
-    print sqsum/denom
+    if res < expected_covariances[name]-eps:
+        print 'Result is better than expected by: %.18e'\
+              %(res-expected_covariances[name])
+        print 'Result = %.18e' %res
+        print 'Expect = %.18e' %expected_covariances[name]    
+    elif res > expected_covariances[name]+eps:
+        print 'FAIL: %.18e' %res
+    else:
+        pass
 
-    from pylab import *
-    ion()
-    hold(False)
-    plot(reference_time, gauge_values[k], 'r.-',
-         reference_time, model, 'k-')
-    title('Gauge %s' %name)
-    xlabel('time(s)')
-    ylabel('stage (m)')    
-    legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
-    savefig(name, dpi = 300)
 
-    raw_input('Next')
-show()
+    if plotting is True:
+        ion()
+        hold(False)
+    
+        plot(reference_time, validation_data[name], 'r.-',
+             reference_time, model, 'k-')
+        title('Gauge %s' %name)
+        xlabel('time(s)')
+        ylabel('stage (m)')    
+        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
+        savefig(name, dpi = 300)        
+
+        raw_input('Next')
 
 
 
-#from pylab import *
-#plot(time, stage)
-#show()
Index: /anuga_validation/okushiri_2005/project.py
===================================================================
--- /anuga_validation/okushiri_2005/project.py	(revision 3849)
+++ /anuga_validation/okushiri_2005/project.py	(revision 3850)
@@ -2,9 +2,18 @@
 """
 
-boundary_filename = 'Benchmark_2_input.txt'
-bathymetry_filename = 'Benchmark_2_Bathymetry.txt'
+# Given boundary wave
+boundary_filename = 'Benchmark_2_input.tms'
+
+# Observed timeseries
+validation_filename = 'output_ch5-7-9.txt'
+
+# A simple conversion from txt file to the more compact NetCDF format
+bathymetry_filename = 'Benchmark_2_Bathymetry.pts'
+
+# Triangular mesh
 mesh_filename = 'Benchmark_2.msh'
 
-output_filename = 'okushiri.sww'
+# Model output
+output_filename = 'okushiri_old_limiters.sww'
 
 
Index: /anuga_validation/okushiri_2005/run_okushiri.py
===================================================================
--- /anuga_validation/okushiri_2005/run_okushiri.py	(revision 3849)
+++ /anuga_validation/okushiri_2005/run_okushiri.py	(revision 3850)
@@ -42,8 +42,8 @@
 domain.set_quantity('stage', 0.0)
 domain.set_quantity('elevation',
-                    filename = project.bathymetry_filename[:-4] + '.pts',
-                    alpha = 0.02,                    
-                    verbose = True,
-                    use_cache = True)
+                    filename=project.bathymetry_filename,
+                    alpha=0.02,                    
+                    verbose=True,
+                    use_cache=True)
 
 
@@ -51,20 +51,9 @@
 # Domain parameters 
 #-------------------------
-domain.set_name(project.output_filename)
-domain.set_default_order(2)
-domain.set_minimum_storable_height(0.001)
-domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
-
-# Set old (pre Sep 2006) defaults for limiters
-#domain.beta_w      = 0.9
-#domain.beta_w_dry  = 0.9
-#domain.beta_uh     = 0.9
-#domain.beta_uh_dry = 0.9
-#domain.beta_vh     = 0.9
-#domain.beta_vh_dry = 0.9
-
-#domain.check_integrity()
-#print domain.statistics()
-
+domain.set_name(project.output_filename)  # Name of output sww file
+domain.set_default_order(2)               # Apply second order scheme 
+domain.set_all_limiters(0.9)              # Maximal second order scheme (old lim)
+domain.set_minimum_storable_height(0.001) # Don't store w < 0.001m
+domain.set_maximum_allowed_speed(0.1)     # Allow a little runoff (0.1 is OK)
 
 
@@ -74,5 +63,5 @@
 
 # Create boundary function from timeseries provided in file
-function = file_function(project.boundary_filename[:-4] + '.tms',
+function = file_function(project.boundary_filename,
                          domain, verbose=True)
 
