Index: /trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py
===================================================================
--- /trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py	(revision 7938)
+++ /trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py	(revision 7939)
@@ -1480,9 +1480,8 @@
                 self.evolve_one_rk3_step(yieldstep, finaltime)
 
-
             #==========================================
             # Apply other fractional steps
             #==========================================
-            #self.apply_fractional_steps()
+            self.apply_fractional_steps()
 
             #==========================================
@@ -1840,4 +1839,11 @@
         msg = 'Method compute_fluxes must be overridden by Domain subclass'
         raise Exception, msg
+
+
+    ##
+    # @brief apply_fractional_step.
+    # Method should be overwritten in subclass, here it does nothing
+    def apply_fractional_steps(self):
+        pass
 
     ##
Index: /trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 7938)
+++ /trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 7939)
@@ -166,8 +166,10 @@
         self.forcing_terms.append(gravity)
 
+
+        self.fractional_step_operators = []
+
         # Stored output
         self.store = True
         self.set_store_vertices_uniquely(False)
-
         self.quantities_to_be_stored = {'elevation': 1, 
                                         'stage': 2, 
@@ -529,4 +531,10 @@
             (either from this module or C-extension)"""
         compute_fluxes(self)
+
+
+    def apply_fractional_steps(self):
+        """Loop throughand apply all extra fractional steps"""
+        for operator in self.fractional_step_operators:
+            operator.apply()
 
     def distribute_to_vertices_and_edges(self):
@@ -1132,4 +1140,5 @@
     elevation = domain.quantities['elevation']
 
+    #FIXME SR Should avoid allocating memory!
     height = stage.centroid_values - elevation.centroid_values
     elevation = elevation.vertex_values
Index: /trunk/anuga_core/source/anuga/structures/culvert_class.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/culvert_class.py	(revision 7939)
+++ /trunk/anuga_core/source/anuga/structures/culvert_class.py	(revision 7939)
@@ -0,0 +1,1627 @@
+import sys
+
+from anuga.shallow_water.forcing import Inflow, General_forcing
+from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
+from anuga.utilities.system_tools import log_to_file
+from anuga.geometry.polygon import inside_polygon, is_inside_polygon, plot_polygons
+
+from anuga.utilities.numerical_tools import mean
+from anuga.utilities.numerical_tools import ensure_numeric, sign
+        
+from anuga.config import g, epsilon
+from anuga.config import minimum_allowed_height, velocity_protection        
+import anuga.utilities.log as log
+
+import numpy as num
+from math import sqrt
+from math import sqrt
+
+class Below_interval(Exception): pass 
+class Above_interval(Exception): pass
+
+# FIXME(Ole): Take a good hard look at logging here
+
+
+# FIXME(Ole): Write in C and reuse this function by similar code 
+# in interpolate.py
+def interpolate_linearly(x, xvec, yvec):
+
+    msg = 'Input to function interpolate_linearly could not be converted '
+    msg += 'to numerical scalar: x = %s' % str(x)
+    try:
+        x = float(x)
+    except:
+        raise Exception, msg
+
+
+    # Check bounds
+    if x < xvec[0]: 
+        msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
+            % (x, xvec[0])
+        raise Below_interval, msg
+        
+    if x > xvec[-1]: 
+        msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
+            %(x, xvec[-1])
+        raise Above_interval, msg        
+        
+        
+    # Find appropriate slot within bounds            
+    i = 0
+    while x > xvec[i]: i += 1
+
+    
+    x0 = xvec[i-1]
+    x1 = xvec[i]            
+    alpha = (x - x0)/(x1 - x0) 
+            
+    y0 = yvec[i-1]
+    y1 = yvec[i]                        
+    y = alpha*y1 + (1-alpha)*y0
+    
+    return y 
+            
+
+            
+def read_culvert_description(culvert_description_filename):            
+    
+    # Read description file
+    fid = open(culvert_description_filename)
+    
+    read_rating_curve_data = False        
+    rating_curve = []
+    for i, line in enumerate(fid.readlines()):
+        
+        if read_rating_curve_data is True:
+            fields = line.split(',')
+            head_difference = float(fields[0].strip())
+            flow_rate = float(fields[1].strip())                
+            barrel_velocity = float(fields[2].strip())
+            
+            rating_curve.append([head_difference, flow_rate, barrel_velocity]) 
+        
+        if i == 0:
+            # Header
+            continue
+        if i == 1:
+            # Metadata
+            fields = line.split(',')
+            label=fields[0].strip()
+            type=fields[1].strip().lower()
+            assert type in ['box', 'pipe']
+            
+            width=float(fields[2].strip())
+            height=float(fields[3].strip())
+            length=float(fields[4].strip())
+            number_of_barrels=int(fields[5].strip())
+            #fields[6] refers to losses
+            description=fields[7].strip()                
+                
+        if line.strip() == '': continue # Skip blanks
+
+        if line.startswith('Rating'):
+            read_rating_curve_data = True
+            # Flow data follows
+                
+    fid.close()
+    
+    return label, type, width, height, length, number_of_barrels, description, rating_curve
+    
+
+    
+
+class Culvert_flow_general:
+    """Culvert flow - transfer water from one hole to another
+    
+    This version will work with either rating curve file or with culvert
+    routine.
+    
+    Input: Two points, pipe_size (either diameter or width, height), 
+    mannings_rougness,
+    """	
+
+    def __init__(self,
+                 domain,
+                 culvert_description_filename=None,
+                 culvert_routine=None,
+                 end_point0=None, 
+                 end_point1=None,
+                 enquiry_point0=None, 
+                 enquiry_point1=None,
+                 type='box',
+                 width=None,
+                 height=None,
+                 length=None,
+                 number_of_barrels=1,
+                 number_of_smoothing_steps=2000,
+                 trigger_depth=0.01, # Depth below which no flow happens
+                 manning=None,          # Mannings Roughness for Culvert 
+                 sum_loss=None,
+                 use_velocity_head=False, # FIXME(Ole): Get rid of - always True
+                 use_momentum_jet=False, # FIXME(Ole): Not yet implemented
+                 label=None,
+                 description=None,
+                 update_interval=None,
+                 log_file=False,
+                 discharge_hydrograph=False,
+                 verbose=False):
+        
+
+        
+        # Input check
+        
+        if height is None: height = width
+        
+        assert number_of_barrels >= 1
+        assert use_velocity_head is True or use_velocity_head is False
+        
+        #msg = 'Momentum jet not yet moved to general culvert'
+        #assert use_momentum_jet is False, msg
+        self.use_momentum_jet = use_momentum_jet
+ 
+        self.culvert_routine = culvert_routine        
+        self.culvert_description_filename = culvert_description_filename
+        if culvert_description_filename is not None:
+            label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
+            self.rating_curve = ensure_numeric(rating_curve)            
+
+        self.height = height
+        self.width = width
+
+        
+        self.domain = domain
+        self.trigger_depth = trigger_depth        
+                
+        if manning is None: 
+            self.manning = 0.012   # Default roughness for pipe
+        
+        if sum_loss is None:
+            self.sum_loss = 0.0
+            
+        
+                        
+        # Store culvert information
+        self.label = label
+        self.description = description
+        self.culvert_type = type
+        self.number_of_barrels = number_of_barrels
+        
+        # Store options        
+        self.use_velocity_head = use_velocity_head
+
+        if label is None: label = 'culvert_flow'
+        label += '_' + str(id(self)) 
+        self.label = label
+        
+        # File for storing discharge_hydrograph
+        if discharge_hydrograph is True:
+            self.timeseries_filename = label + '_timeseries.csv'
+            fid = open(self.timeseries_filename, 'w')
+            fid.write('time, discharge\n')
+            fid.close()
+
+        # Log file for storing general textual output
+        if log_file is True:        
+            self.log_filename = label + '.log'         
+            log_to_file(self.log_filename, self.label)        
+            log_to_file(self.log_filename, description)
+            log_to_file(self.log_filename, self.culvert_type)        
+        else:
+            self.log_filename = None
+
+
+        # Create the fundamental culvert polygons from polygon
+        P = create_culvert_polygons(end_point0,
+                                    end_point1,
+                                    width=width,   
+                                    height=height,
+                                    number_of_barrels=number_of_barrels)
+        self.culvert_polygons = P
+        
+        # Select enquiry points
+        if enquiry_point0 is None:
+            enquiry_point0 = P['enquiry_point0']
+            
+        if enquiry_point1 is None:
+            enquiry_point1 = P['enquiry_point1']            
+            
+        if verbose is True:
+            pass
+            #plot_polygons([[end_point0, end_point1],
+            #               P['exchange_polygon0'],
+            #               P['exchange_polygon1'],
+            #               [enquiry_point0, 1.005*enquiry_point0],
+            #               [enquiry_point1, 1.005*enquiry_point1]],
+            #              figname='culvert_polygon_output')
+
+            
+            
+        self.enquiry_points = [enquiry_point0, enquiry_point1]
+        self.enquiry_indices = self.get_enquiry_indices()                  
+        self.check_culvert_inside_domain()
+        
+                    
+        # Create inflow object at each end of the culvert. 
+        self.openings = []
+        self.openings.append(Inflow(domain,
+                                    polygon=P['exchange_polygon0']))
+        self.openings.append(Inflow(domain,
+                                    polygon=P['exchange_polygon1']))
+                                    
+        # Assume two openings for now: Referred to as 0 and 1
+        assert len(self.openings) == 2
+
+        # Establish initial values at each enquiry point
+        dq = domain.quantities 
+        for i, opening in enumerate(self.openings):
+            idx = self.enquiry_indices[i]
+            elevation = dq['elevation'].get_values(location='centroids',
+                                                   indices=[idx])[0]
+            stage = dq['stage'].get_values(location='centroids',
+                                           indices=[idx])[0]
+            opening.elevation = elevation
+            opening.stage = stage
+            opening.depth = stage-elevation
+
+            
+                            
+        # Determine initial pipe direction.
+        # This may change dynamically based on the total energy difference     
+        # Consequently, this may be superfluous
+        delta_z = self.openings[0].elevation - self.openings[1].elevation
+        if delta_z > 0.0:
+            self.inlet = self.openings[0]
+            self.outlet = self.openings[1]
+        else:                
+            self.outlet = self.openings[0]
+            self.inlet = self.openings[1]        
+        
+        
+        # Store basic geometry 
+        self.end_points = [end_point0, end_point1]
+        self.vector = P['vector']
+        self.length = P['length']; assert self.length > 0.0
+        if culvert_description_filename is not None:
+            if not num.allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
+                msg = 'WARNING: barrel length specified in "%s" (%.2f m)'\
+                    % (culvert_description_filename, 
+                       length)
+                msg += ' does not match distance between specified'
+                msg += ' end points (%.2f m)' %self.length
+                log.critical(msg)
+        
+        self.verbose = verbose
+
+        # Circular index for flow averaging in culvert
+        self.N = N = number_of_smoothing_steps
+        self.Q_list = [0]*N
+        self.i = i
+        
+        # For use with update_interval                        
+        self.last_update = 0.0
+        self.update_interval = update_interval
+        
+        # Create objects to update momentum (a bit crude at this stage). This is used with the momentum jet.
+        xmom0 = General_forcing(domain, 'xmomentum',
+                                polygon=P['exchange_polygon0'])
+
+        xmom1 = General_forcing(domain, 'xmomentum',
+                                polygon=P['exchange_polygon1'])
+
+        ymom0 = General_forcing(domain, 'ymomentum',
+                                polygon=P['exchange_polygon0'])
+
+        ymom1 = General_forcing(domain, 'ymomentum',
+                                polygon=P['exchange_polygon1'])
+
+        self.opening_momentum = [[xmom0, ymom0], [xmom1, ymom1]]
+
+
+
+        # Print some diagnostics to log if requested
+        if self.log_filename is not None:
+            s = 'Culvert Effective Length = %.2f m' %(self.length)
+            log_to_file(self.log_filename, s)
+   
+            s = 'Culvert Direction is %s\n' %str(self.vector)
+            log_to_file(self.log_filename, s)        
+        
+        
+        
+        
+        
+    def __call__(self, domain):
+
+        # Time stuff
+        time = domain.get_time()
+        
+        
+        update = False
+        if self.update_interval is None:
+            # Use next timestep as has been computed in domain.py        
+            delta_t = domain.timestep            
+            update = True
+        else:    
+            # Use update interval 
+            delta_t = self.update_interval            
+            if time - self.last_update > self.update_interval or time == 0.0:
+                update = True
+            
+        if self.log_filename is not None:        
+            s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
+            log_to_file(self.log_filename, s)
+                
+                                
+        if update is True:
+            self.compute_rates(delta_t)
+            
+    
+        # Execute flow term for each opening
+        # This is where Inflow objects are evaluated using the last rate 
+        # that has been calculated
+        # 
+        # This will take place at every internal timestep and update the domain
+        for opening in self.openings:
+            opening(domain)
+            
+
+
+    def get_enquiry_indices(self):
+        """Get indices for nearest centroids to self.enquiry_points
+        """
+        
+        domain = self.domain
+        
+        enquiry_indices = []                  
+        for point in self.enquiry_points:
+            # Find nearest centroid 
+            N = len(domain)    
+            points = domain.get_centroid_coordinates(absolute=True)
+
+            # Calculate indices in exchange area for this forcing term
+            
+            triangle_id = min_dist = sys.maxint
+            for k in range(N):
+                x, y = points[k,:] # Centroid
+
+                c = point                                
+                distance = (x-c[0])**2+(y-c[1])**2
+                if distance < min_dist:
+                    min_dist = distance
+                    triangle_id = k
+
+                    
+            if triangle_id < sys.maxint:
+                msg = 'found triangle with centroid (%f, %f)'\
+                    %tuple(points[triangle_id, :])
+                msg += ' for point (%f, %f)' %tuple(point)
+                
+                enquiry_indices.append(triangle_id)
+            else:
+                msg = 'Triangle not found for point (%f, %f)' %point
+                raise Exception, msg
+        
+        return enquiry_indices
+
+        
+    def check_culvert_inside_domain(self):
+        """Check that all polygons and enquiry points lie within the mesh.
+        """
+        bounding_polygon = self.domain.get_boundary_polygon()
+        P = self.culvert_polygons
+        for key in P.keys():
+            if key in ['exchange_polygon0', 
+                       'exchange_polygon1']:
+                for point in list(P[key]) + self.enquiry_points:
+                    msg = 'Point %s in polygon %s for culvert %s did not'\
+                        %(str(point), key, self.label)
+                    msg += 'fall within the domain boundary.'
+                    assert is_inside_polygon(point, bounding_polygon), msg
+            
+
+    def adjust_flow_for_available_water_at_inlet(self, Q, delta_t):
+        """Adjust Q downwards depending on available water at inlet
+
+           This is a critical step in modelling bridges and Culverts
+           the predicted flow through a structure based on an abstract
+           algorithm can at times request for water that is simply not
+           available due to any number of constrictions that limit the
+           flow approaching the structure In order to ensure that
+           there is adequate flow available certain checks are
+           required There needs to be a check using the Static Water
+           Volume sitting infront of the structure, In addition if the
+           water is moving the available water will be larger than the
+           static volume
+           
+           NOTE To temporarily switch this off for Debugging purposes
+           rem out line in function def compute_rates below
+        """
+    
+        if delta_t < epsilon:
+            # No need to adjust if time step is very small or zero
+            # In this case the possible flow will be very large
+            # anyway.
+            return Q
+        
+        # Short hands
+        domain = self.domain        
+        dq = domain.quantities                
+        time = domain.get_time()
+        I = self.inlet
+        idx = I.exchange_indices    
+
+        # Find triangle with the smallest depth
+        stage = dq['stage'].get_values(location='centroids', 
+                                               indices=[idx])
+        elevation = dq['elevation'].get_values(location='centroids', 
+                                               indices=[idx])        
+        depth = stage-elevation
+        min_depth = min(depth.flat)  # This may lead to errors if edge of area is at a higher level !!!!
+        avg_depth = mean(depth.flat) # Yes, but this one violates the conservation unit tests
+
+
+
+        # FIXME (Ole): If you want these, use log.critical() and
+        # make the statements depend on verbose
+        #print I.depth
+        #print I.velocity
+        #print self.width
+
+        # max_Q Based on Volume Calcs
+
+
+        depth_term = min_depth*I.exchange_area/delta_t
+        if min_depth < 0.2:
+            # Only add velocity term in shallow waters (< 20 cm)
+            # This is a little ad hoc, but maybe it is reasonable
+            velocity_term = self.width*min_depth*I.velocity
+        else:
+            velocity_term = 0.0
+
+        # This one takes approaching water into account    
+        max_Q = max(velocity_term, depth_term)
+
+        # This one preserves Volume
+        #max_Q = depth_term
+
+
+        if self.verbose is True:
+            log.critical('Max_Q = %f' % max_Q)            
+            msg = 'Width = %.2fm, Depth at inlet = %.2f m, Velocity = %.2f m/s.      ' % (self.width, I.depth, I.velocity)
+            msg += 'Max Q = %.2f m^3/s' %(max_Q)
+            log.critical(msg)
+
+        if self.log_filename is not None:                
+            log_to_file(self.log_filename, msg)
+        # New Procedure for assessing the flow available to the Culvert 
+        # This routine uses the GET FLOW THROUGH CROSS SECTION
+        #   Need to check Several Polyline however
+        # Firstly 3 sides of the exchange Poly
+        # then only the Line Directly infront of the Polygon
+        # Access polygon Points from   self.inlet.polygon
+      
+        #  The Following computes the flow crossing over 3 sides of the exchange polygon for the structure
+        # Clearly the flow in the culvert can not be more than that flowing toward it through the exhange polygon
+        
+        #q1 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][1:3]) # First Side Segment
+        #q2 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][2:])   # Second Face Segment
+        #q3 =domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'].take([3,0], axis=0)) # Third Side Segment
+        # q4 = domain.get_flow_through_cross_section([self.culvert_polygons['exchange_polygon0'][1:4]][0])
+        #q4=max(q1,0.0)+max(q2,0.0)+max(q3,0.0)
+        # To use only the Flow crossing the 3 sides of the Exchange Polygon use the following Line Only
+        #max_Q=max(q1,q2,q3,q4)
+        # Try Simple Smoothing using Average of 2 approaches
+        #max_Q=(max(q1,q2,q3,q4)+max_Q)/2.0
+        # Calculate the minimum in absolute terms of
+        # the requsted flow and the possible flow
+        Q_reduced = sign(Q)*min(abs(Q), abs(max_Q))
+        if self.verbose is True:
+            msg = 'Initial Q Reduced = %.2f m3/s.      ' % (Q_reduced)
+            log.critical(msg)
+
+        if self.log_filename is not None:                
+            log_to_file(self.log_filename, msg)
+        # Now Keep Rolling Average of Computed Discharge to Reduce / Remove Oscillations
+        #  can use delta_t if we want to averageover a time frame for example
+        # N = 5.0/delta_t  Will provide the average over 5 seconds
+
+        self.i=(self.i+1)%self.N
+        self.Q_list[self.i]=Q_reduced
+        Q_reduced = sum(self.Q_list)/len(self.Q_list)
+
+        if self.verbose is True:
+            msg = 'Final Q Reduced = %.2f m3/s.      ' % (Q_reduced)
+            log.critical(msg)
+
+        if self.log_filename is not None:                
+            log_to_file(self.log_filename, msg)
+
+
+        if abs(Q_reduced) < abs(Q): 
+            msg = '%.2fs: Requested flow is ' % time
+            msg += 'greater than what is supported by the smallest '
+            msg += 'depth at inlet exchange area:\n        '
+            msg += 'inlet exchange area: %.2f '% (I.exchange_area) 
+            msg += 'velocity at inlet :%.2f '% (I.velocity)
+            msg += 'Vel* Exch Area = : %.2f '% (I.velocity*avg_depth*self.width)
+            msg += 'h_min*inlet_area/delta_t = %.2f*%.2f/%.2f '\
+                % (avg_depth, I.exchange_area, delta_t)
+            msg += ' = %.2f m^3/s\n        ' % Q_reduced
+            msg += 'Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
+            msg += 'Note calculate max_Q from V %.2f m^3/s ' % (max_Q)
+            if self.verbose is True:
+                log.critical(msg)
+                
+            if self.log_filename is not None:                
+                log_to_file(self.log_filename, msg)
+        
+        return Q_reduced    
+                        
+            
+    def compute_rates(self, delta_t):
+        """Compute new rates for inlet and outlet
+        """
+
+        # Short hands
+        domain = self.domain        
+        dq = domain.quantities                
+        
+        # Time stuff
+        time = domain.get_time()
+        self.last_update = time
+
+            
+        if hasattr(self, 'log_filename'):
+            log_filename = self.log_filename
+            
+        # Compute stage, energy and velocity at the 
+        # enquiry points at each end of the culvert
+        openings = self.openings
+        for i, opening in enumerate(openings):
+            idx = self.enquiry_indices[i]                
+            
+            stage = dq['stage'].get_values(location='centroids',
+                                           indices=[idx])[0]
+            depth = h = stage-opening.elevation
+                                                           
+            
+            # Get velocity                                 
+            xmomentum = dq['xmomentum'].get_values(location='centroids',
+                                                   indices=[idx])[0]
+            ymomentum = dq['xmomentum'].get_values(location='centroids',
+                                                   indices=[idx])[0]
+
+            if h > minimum_allowed_height:
+                u = xmomentum/(h + velocity_protection/h)
+                v = ymomentum/(h + velocity_protection/h)
+            else:
+                u = v = 0.0
+                
+            v_squared = u*u + v*v
+            
+            if self.use_velocity_head is True:
+                velocity_head = 0.5*v_squared/g    
+            else:
+                velocity_head = 0.0
+            
+            opening.total_energy = velocity_head + stage
+            opening.specific_energy = velocity_head + depth
+            opening.stage = stage
+            opening.depth = depth
+            opening.velocity = sqrt(v_squared)
+            
+
+        # We now need to deal with each opening individually
+        # Determine flow direction based on total energy difference
+        delta_total_energy = openings[0].total_energy - openings[1].total_energy
+        if delta_total_energy > 0:
+            inlet = openings[0]
+            outlet = openings[1]
+
+            # FIXME: I think this whole momentum jet thing could be a bit more elegant
+            inlet.momentum = self.opening_momentum[0]
+            outlet.momentum = self.opening_momentum[1]
+        else:
+            inlet = openings[1]
+            outlet = openings[0]
+            
+            inlet.momentum = self.opening_momentum[1]
+            outlet.momentum = self.opening_momentum[0]
+
+            delta_total_energy = -delta_total_energy
+
+        self.inlet = inlet
+        self.outlet = outlet
+            
+        msg = 'Total energy difference is negative'
+        assert delta_total_energy >= 0.0, msg
+
+        # Recompute slope and issue warning if flow is uphill
+        # These values do not enter the computation
+        delta_z = inlet.elevation - outlet.elevation
+        culvert_slope = (delta_z/self.length)
+        if culvert_slope < 0.0:
+            # Adverse gradient - flow is running uphill
+            # Flow will be purely controlled by uphill outlet face
+            if self.verbose is True:
+                log.critical('%.2fs - WARNING: Flow is running uphill.' % time)
+            
+        if self.log_filename is not None:
+            s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
+                %(time, self.inlet.stage, self.outlet.stage)
+            log_to_file(self.log_filename, s)
+            s = 'Delta total energy = %.3f' %(delta_total_energy)
+            log_to_file(log_filename, s)
+
+            
+        # Determine controlling energy (driving head) for culvert
+        if inlet.specific_energy > delta_total_energy:
+            # Outlet control
+            driving_head = delta_total_energy
+        else:
+            # Inlet control
+            driving_head = inlet.specific_energy
+            
+
+
+        if self.inlet.depth <= self.trigger_depth:
+            Q = 0.0
+        else:
+            # Calculate discharge for one barrel and 
+            # set inlet.rate and outlet.rate
+            
+            if self.culvert_description_filename is not None:
+                try:
+                    Q = interpolate_linearly(driving_head, 
+                                             self.rating_curve[:,0], 
+                                             self.rating_curve[:,1]) 
+                except Below_interval, e:
+                    Q = self.rating_curve[0,1]             
+                    msg = '%.2fs: ' % time 
+                    msg += 'Delta head smaller than rating curve minimum: '
+                    msg += str(e)
+                    msg += '\n        '
+                    msg += 'I will use minimum discharge %.2f m^3/s ' % Q
+                    msg += 'for culvert "%s"' % self.label
+                    
+                    if hasattr(self, 'log_filename'):                    
+                        log_to_file(self.log_filename, msg)
+                except Above_interval, e:
+                    Q = self.rating_curve[-1,1]          
+                    msg = '%.2fs: ' % time                 
+                    msg += 'Delta head greater than rating curve maximum: '
+                    msg += str(e)
+                    msg += '\n        '
+                    msg += 'I will use maximum discharge %.2f m^3/s ' % Q
+                    msg += 'for culvert "%s"' % self.label 
+                    
+                    if self.log_filename is not None:                    
+                        log_to_file(self.log_filename, msg)
+            else:
+                # User culvert routine
+                Q, barrel_velocity, culvert_outlet_depth =\
+                    self.culvert_routine(inlet.depth,
+                                         outlet.depth,
+                                         inlet.velocity,
+                                         outlet.velocity,
+                                         inlet.specific_energy, 
+                                         delta_total_energy, 
+                                         g,
+                                         culvert_length=self.length,
+                                         culvert_width=self.width,
+                                         culvert_height=self.height,
+                                         culvert_type=self.culvert_type,
+                                         manning=self.manning,
+                                         sum_loss=self.sum_loss,
+                                         log_filename=self.log_filename)
+            
+            
+        
+        # Adjust discharge for multiple barrels
+        Q *= self.number_of_barrels
+
+        # Adjust discharge for available water at the inlet
+        Q = self.adjust_flow_for_available_water_at_inlet(Q, delta_t)
+        
+        self.inlet.rate = -Q
+        self.outlet.rate = Q
+
+
+        # Momentum jet stuff
+        if self.use_momentum_jet is True:
+
+
+            # Compute barrel momentum
+            barrel_momentum = barrel_velocity*culvert_outlet_depth
+
+            if self.log_filename is not None:                                    
+                s = 'Barrel velocity = %f' %barrel_velocity
+                log_to_file(self.log_filename, s)
+
+            # Compute momentum vector at outlet
+            outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
+                
+            if self.log_filename is not None:                
+                s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
+                log_to_file(self.log_filename, s)
+
+
+            # Update momentum        
+            if delta_t > 0.0:
+                xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
+                xmomentum_rate /= delta_t
+                    
+                ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
+                ymomentum_rate /= delta_t
+                        
+                if self.log_filename is not None:                
+                    s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
+                    log_to_file(self.log_filename, s)                    
+            else:
+                xmomentum_rate = ymomentum_rate = 0.0
+
+
+            # Set momentum rates for outlet jet
+            outlet.momentum[0].rate = xmomentum_rate
+            outlet.momentum[1].rate = ymomentum_rate
+
+            # Remember this value for next step (IMPORTANT)
+            outlet.momentum[0].value = outlet_mom_x
+            outlet.momentum[1].value = outlet_mom_y                    
+
+            if int(domain.time*100) % 100 == 0:
+
+                if self.log_filename is not None:                
+                    s = 'T=%.5f, Culvert Discharge = %.3f f'\
+                        %(time, Q)
+                    s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
+                        %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
+                    s += ' Momentum rate: (%.4f, %.4f)'\
+                        %(xmomentum_rate, ymomentum_rate)                    
+                    s+='Outlet Vel= %.3f'\
+                        %(barrel_velocity)
+                    log_to_file(self.log_filename, s)
+
+
+            # Execute momentum terms
+            # This is where Inflow objects are evaluated and update the domain
+                self.outlet.momentum[0](domain)
+                self.outlet.momentum[1](domain)        
+            
+
+            
+        # Log timeseries to file
+        try:
+            fid = open(self.timeseries_filename, 'a')        
+        except:
+            pass
+        else:    
+            fid.write('%.2f, %.2f\n' %(time, Q))
+            fid.close()
+
+        # Store value of time
+        self.last_time = time
+
+
+            
+
+
+
+                            
+# OBSOLETE (Except for momentum jet in Culvert_flow_energy)    
+class Culvert_flow_rating:
+    """Culvert flow - transfer water from one hole to another
+    
+
+    Input: Two points, pipe_size (either diameter or width, height), 
+    mannings_rougness,
+    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
+    top-down_blockage_factor and bottom_up_blockage_factor
+    
+    """	
+
+    def __init__(self,
+                 domain,
+                 culvert_description_filename=None,
+                 end_point0=None, 
+                 end_point1=None,
+                 enquiry_point0=None, 
+                 enquiry_point1=None,
+                 update_interval=None,
+                 log_file=False,
+                 discharge_hydrograph=False,
+                 verbose=False):
+        
+
+        
+        label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
+        
+                
+        # Store culvert information
+        self.label = label
+        self.description = description
+        self.culvert_type = type
+        self.rating_curve = ensure_numeric(rating_curve)
+        self.number_of_barrels = number_of_barrels
+
+        if label is None: label = 'culvert_flow'
+        label += '_' + str(id(self)) 
+        self.label = label
+        
+        # File for storing discharge_hydrograph
+        if discharge_hydrograph is True:
+            self.timeseries_filename = label + '_timeseries.csv'
+            fid = open(self.timeseries_filename, 'w')
+            fid.write('time, discharge\n')
+            fid.close()
+
+        # Log file for storing general textual output
+        if log_file is True:        
+            self.log_filename = label + '.log'         
+            log_to_file(self.log_filename, self.label)        
+            log_to_file(self.log_filename, description)
+            log_to_file(self.log_filename, self.culvert_type)        
+
+
+        # Create the fundamental culvert polygons from POLYGON
+        #if self.culvert_type == 'circle':
+        #    # Redefine width and height for use with create_culvert_polygons
+        #    width = height = diameter
+        
+        P = create_culvert_polygons(end_point0,
+                                    end_point1,
+                                    width=width,   
+                                    height=height,
+                                    number_of_barrels=number_of_barrels)
+        
+        # Select enquiry points
+        if enquiry_point0 is None:
+            enquiry_point0 = P['enquiry_point0']
+            
+        if enquiry_point1 is None:
+            enquiry_point1 = P['enquiry_point1']            
+            
+        if verbose is True:
+            pass
+            #plot_polygons([[end_point0, end_point1],
+            #               P['exchange_polygon0'],
+            #               P['exchange_polygon1'],
+            #               [enquiry_point0, 1.005*enquiry_point0],
+            #               [enquiry_point1, 1.005*enquiry_point1]],                           
+            #              figname='culvert_polygon_output')
+
+            
+            
+        self.enquiry_points = [enquiry_point0, enquiry_point1]                           
+
+        self.enquiry_indices = []                  
+        for point in self.enquiry_points:
+            # Find nearest centroid 
+            N = len(domain)    
+            points = domain.get_centroid_coordinates(absolute=True)
+
+            # Calculate indices in exchange area for this forcing term
+            
+            triangle_id = min_dist = sys.maxint
+            for k in range(N):
+                x, y = points[k,:] # Centroid
+
+                c = point                                
+                distance = (x-c[0])**2+(y-c[1])**2
+                if distance < min_dist:
+                    min_dist = distance
+                    triangle_id = k
+
+                    
+            if triangle_id < sys.maxint:
+                msg = 'found triangle with centroid (%f, %f)'\
+                    %tuple(points[triangle_id, :])
+                msg += ' for point (%f, %f)' %tuple(point)
+                
+                self.enquiry_indices.append(triangle_id)
+            else:
+                msg = 'Triangle not found for point (%f, %f)' %point
+                raise Exception, msg
+        
+                          
+
+        # Check that all polygons lie within the mesh.
+        bounding_polygon = domain.get_boundary_polygon()
+        for key in P.keys():
+            if key in ['exchange_polygon0', 
+                       'exchange_polygon1']:
+                for point in list(P[key]) + self.enquiry_points:
+                    msg = 'Point %s in polygon %s for culvert %s did not'\
+                        %(str(point), key, self.label)
+                    msg += 'fall within the domain boundary.'
+                    assert is_inside_polygon(point, bounding_polygon), msg
+        
+                    
+        # Create inflow object at each end of the culvert. 
+        self.openings = []
+        self.openings.append(Inflow(domain,
+                                    polygon=P['exchange_polygon0']))
+
+        self.openings.append(Inflow(domain,
+                                    polygon=P['exchange_polygon1']))                                    
+
+
+
+        dq = domain.quantities                                            
+        for i, opening in enumerate(self.openings):                            
+            elevation = dq['elevation'].get_values(location='centroids',
+                                                   indices=[self.enquiry_indices[i]])            
+            opening.elevation = elevation
+            opening.stage = elevation # Simple assumption that culvert is dry initially
+
+        # Assume two openings for now: Referred to as 0 and 1
+        assert len(self.openings) == 2
+                            
+        # Determine pipe direction     
+        self.delta_z = delta_z = self.openings[0].elevation - self.openings[1].elevation
+        if delta_z > 0.0:
+            self.inlet = self.openings[0]
+            self.outlet = self.openings[1]
+        else:                
+            self.outlet = self.openings[0]
+            self.inlet = self.openings[1]        
+        
+        
+        # Store basic geometry 
+        self.end_points = [end_point0, end_point1]
+        self.vector = P['vector']
+        self.length = P['length']; assert self.length > 0.0
+        if not num.allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
+            msg = 'WARNING: barrel length specified in "%s" (%.2f m)' %(culvert_description_filename, length)
+            msg += ' does not match distance between specified'
+            msg += ' end points (%.2f m)' %self.length
+            log.critical(msg)
+        
+        self.verbose = verbose
+        self.last_update = 0.0 # For use with update_interval        
+        self.last_time = 0.0                
+        self.update_interval = update_interval
+        
+
+        # Print something
+        if hasattr(self, 'log_filename'):
+            s = 'Culvert Effective Length = %.2f m' %(self.length)
+            log_to_file(self.log_filename, s)
+   
+            s = 'Culvert Direction is %s\n' %str(self.vector)
+            log_to_file(self.log_filename, s)        
+        
+        
+        
+        
+        
+    def __call__(self, domain):
+
+        # Time stuff
+        time = domain.get_time()
+        
+        
+        update = False
+        if self.update_interval is None:
+            update = True
+            delta_t = domain.timestep # Next timestep has been computed in domain.py
+        else:    
+            if time - self.last_update > self.update_interval or time == 0.0:
+                update = True
+            delta_t = self.update_interval
+            
+        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
+        if hasattr(self, 'log_filename'):            
+            log_to_file(self.log_filename, s)
+                
+                                
+        if update is True:
+            self.last_update = time
+       
+            dq = domain.quantities
+                        
+            # Get average water depths at each opening        
+            openings = self.openings   # There are two Opening [0] and [1]
+            for i, opening in enumerate(openings):
+                
+                # Compute mean values of selected quantitites in the 
+                # enquiry area in front of the culvert
+                
+                stage = dq['stage'].get_values(location='centroids',
+                                               indices=[self.enquiry_indices[i]])
+                
+                # Store current average stage and depth with each opening object
+                opening.depth = stage - opening.elevation
+                opening.stage = stage
+
+                
+
+            #################  End of the FOR loop ################################################
+
+            # We now need to deal with each opening individually
+                
+            # Determine flow direction based on total energy difference
+
+            delta_w = self.inlet.stage - self.outlet.stage
+            
+            if hasattr(self, 'log_filename'):
+                s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f' %(time, 
+                                                                           self.inlet.stage,
+                                                                           self.outlet.stage)
+                log_to_file(self.log_filename, s)
+
+
+            if self.inlet.depth <= 0.01:
+                Q = 0.0
+            else:
+                # Calculate discharge for one barrel and set inlet.rate and outlet.rate
+                
+                try:
+                    Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1]) 
+                except Below_interval, e:
+                    Q = self.rating_curve[0,1]             
+                    msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time
+                    msg += str(e)
+                    msg += '\n        I will use minimum discharge %.2f m^3/s for culvert "%s"'\
+                        %(Q, self.label)
+                    if hasattr(self, 'log_filename'):                    
+                        log_to_file(self.log_filename, msg)
+                except Above_interval, e:
+                    Q = self.rating_curve[-1,1]          
+                    msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
+                    msg += str(e)
+                    msg += '\n        I will use maximum discharge %.2f m^3/s for culvert "%s"'\
+                        %(Q, self.label)
+                    if hasattr(self, 'log_filename'):                    
+                        log_to_file(self.log_filename, msg)                    
+
+                
+                
+            
+            # Adjust discharge for multiple barrels
+            Q *= self.number_of_barrels
+            
+
+            # Adjust Q downwards depending on available water at inlet
+            stage = self.inlet.get_quantity_values(quantity_name='stage')
+            elevation = self.inlet.get_quantity_values(quantity_name='elevation')
+            depth = stage-elevation
+            
+            
+            V = 0
+            for i, d in enumerate(depth):
+                V += d * domain.areas[i]
+            
+            dt = delta_t            
+            if Q*dt > V:
+            
+                Q_reduced = 0.9*V/dt # Reduce with safety factor
+                
+                msg = '%.2fs: Computed extraction for this time interval (Q*dt) is ' % time
+                msg += 'greater than current volume (V) at inlet.\n'
+                msg += ' Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
+                
+                if self.verbose is True:
+                    log.critical(msg)
+                if hasattr(self, 'log_filename'):                    
+                    log_to_file(self.log_filename, msg)                                        
+                
+                Q = Q_reduced
+        
+            self.inlet.rate = -Q
+            self.outlet.rate = Q
+
+            # Log timeseries to file
+            try:
+                fid = open(self.timeseries_filename, 'a')        
+            except:
+                pass
+            else:    
+                fid.write('%.2f, %.2f\n' %(time, Q))
+                fid.close()
+
+            # Store value of time
+            self.last_time = time
+
+            
+    
+        # Execute flow term for each opening
+        # This is where Inflow objects are evaluated using the last rate that has been calculated
+        # 
+        # This will take place at every internal timestep and update the domain
+        for opening in self.openings:
+            opening(domain)
+
+
+
+        
+        
+
+class Culvert_flow_energy:
+    """Culvert flow - transfer water from one hole to another
+    
+    Using Momentum as Calculated by Culvert Flow !!
+    Could be Several Methods Investigated to do This !!!
+
+    2008_May_08
+    To Ole:
+    OK so here we need to get the Polygon Creating code to create 
+    polygons for the culvert Based on
+    the two input Points (X0,Y0) and (X1,Y1) - need to be passed 
+    to create polygon
+
+    The two centers are now passed on to create_polygon.
+    
+
+    Input: Two points, pipe_size (either diameter or width, height), 
+    mannings_rougness,
+    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
+    top-down_blockage_factor and bottom_up_blockage_factor
+    
+    
+    And the Delta H enquiry should be change from Openings in line 412 
+    to the enquiry Polygons infront of the culvert
+    At the moment this script uses only Depth, later we can change it to 
+    Energy...
+
+    Once we have Delta H can calculate a Flow Rate and from Flow Rate 
+    an Outlet Velocity
+    The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
+	
+    Invert levels are optional. If left out they will default to the 
+    elevation at the opening.
+        
+    """	
+
+    def __init__(self,
+                 domain,
+                 label=None, 
+                 description=None,
+                 end_point0=None, 
+                 end_point1=None,
+                 width=None,
+                 height=None,
+                 diameter=None,
+                 manning=None,          # Mannings Roughness for Culvert
+                 invert_level0=None,    # Invert level at opening 0
+                 invert_level1=None,    # Invert level at opening 1
+                 loss_exit=None,
+                 loss_entry=None,
+                 loss_bend=None,
+                 loss_special=None,
+                 blockage_topdwn=None,
+                 blockage_bottup=None,
+                 culvert_routine=None,
+                 number_of_barrels=1,
+                 enquiry_point0=None, 
+                 enquiry_point1=None,
+                 update_interval=None,
+                 verbose=False):
+        
+        # Input check
+        if diameter is not None:
+            self.culvert_type = 'circle'
+            self.diameter = diameter
+            if height is not None or width is not None:
+                msg = 'Either diameter or width&height must be specified, '
+                msg += 'but not both.'
+                raise Exception, msg
+        else:
+            if height is not None:
+                if width is None:
+                    self.culvert_type = 'square'                                
+                    width = height
+                else:
+                    self.culvert_type = 'rectangle'
+            elif width is not None:
+                if height is None:
+                    self.culvert_type = 'square'                                
+                    height = width
+            else:
+                msg = 'Either diameter or width&height must be specified.'
+                raise Exception, msg                
+                
+            if height == width:
+                self.culvert_type = 'square'                                                
+                
+            self.height = height
+            self.width = width
+
+            
+        assert self.culvert_type in ['circle', 'square', 'rectangle']
+        
+        assert number_of_barrels >= 1
+        self.number_of_barrels = number_of_barrels
+        
+        
+        # Set defaults
+        if manning is None: manning = 0.012   # Default roughness for pipe
+        if loss_exit is None: loss_exit = 1.00
+        if loss_entry is None: loss_entry = 0.50
+        if loss_bend is None: loss_bend=0.00
+        if loss_special is None: loss_special=0.00
+        if blockage_topdwn is None: blockage_topdwn=0.00
+        if blockage_bottup is None: blockage_bottup=0.00
+        if culvert_routine is None: 
+            culvert_routine=boyd_generalised_culvert_model
+            
+        if label is None: label = 'culvert_flow'
+        label += '_' + str(id(self)) 
+        self.label = label
+        
+        # File for storing culvert quantities
+        self.timeseries_filename = label + '_timeseries.csv'
+        fid = open(self.timeseries_filename, 'w')
+        fid.write('time, E0, E1, Velocity, Discharge\n')
+        fid.close()
+
+        # Log file for storing general textual output
+        self.log_filename = label + '.log'         
+        log_to_file(self.log_filename, self.label)        
+        log_to_file(self.log_filename, description)
+        log_to_file(self.log_filename, self.culvert_type)        
+
+
+        # Create the fundamental culvert polygons from POLYGON
+        if self.culvert_type == 'circle':
+            # Redefine width and height for use with create_culvert_polygons
+            width = height = diameter
+        
+        P = create_culvert_polygons(end_point0,
+                                    end_point1,
+                                    width=width,   
+                                    height=height,
+                                    number_of_barrels=number_of_barrels)
+        
+        # Select enquiry points
+        if enquiry_point0 is None:
+            enquiry_point0 = P['enquiry_point0']
+            
+        if enquiry_point1 is None:
+            enquiry_point1 = P['enquiry_point1']            
+            
+        if verbose is True:
+            pass
+            #plot_polygons([[end_point0, end_point1],
+            #               P['exchange_polygon0'],
+            #               P['exchange_polygon1'],
+            #               [enquiry_point0, 1.005*enquiry_point0],
+            #               [enquiry_point1, 1.005*enquiry_point1]],
+            #              figname='culvert_polygon_output')
+
+
+        self.enquiry_points = [enquiry_point0, enquiry_point1]                           
+        
+        
+        self.enquiry_indices = []                  
+        for point in self.enquiry_points:
+            # Find nearest centroid 
+            N = len(domain)    
+            points = domain.get_centroid_coordinates(absolute=True)
+
+            # Calculate indices in exchange area for this forcing term
+            
+            triangle_id = min_dist = sys.maxint
+            for k in range(N):
+                x, y = points[k,:] # Centroid
+
+                c = point                                
+                distance = (x-c[0])**2+(y-c[1])**2
+                if distance < min_dist:
+                    min_dist = distance
+                    triangle_id = k
+
+                    
+            if triangle_id < sys.maxint:
+                msg = 'found triangle with centroid (%f, %f)'\
+                    %tuple(points[triangle_id, :])
+                msg += ' for point (%f, %f)' %tuple(point)
+                
+                self.enquiry_indices.append(triangle_id)
+            else:
+                msg = 'Triangle not found for point (%f, %f)' %point
+                raise Exception, msg
+        
+                          
+
+            
+            
+
+        # Check that all polygons lie within the mesh.
+        bounding_polygon = domain.get_boundary_polygon()
+        for key in P.keys():
+            if key in ['exchange_polygon0', 
+                       'exchange_polygon1']:
+                for point in P[key]:
+                
+                    msg = 'Point %s in polygon %s for culvert %s did not'\
+                        %(str(point), key, self.label)
+                    msg += 'fall within the domain boundary.'
+                    assert is_inside_polygon(point, bounding_polygon), msg
+        
+
+        # Create inflow object at each end of the culvert. 
+        self.openings = []
+        self.openings.append(Inflow(domain,
+                                    polygon=P['exchange_polygon0']))
+
+        self.openings.append(Inflow(domain,
+                                    polygon=P['exchange_polygon1']))                                    
+
+
+        # Assume two openings for now: Referred to as 0 and 1
+        assert len(self.openings) == 2
+        
+        # Store basic geometry 
+        self.end_points = [end_point0, end_point1]
+        self.invert_levels = [invert_level0, invert_level1]                
+        #self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
+        #self.enquiry_polylines = [P['enquiry_polygon0'][:2], 
+        #                          P['enquiry_polygon1'][:2]]
+        self.vector = P['vector']
+        self.length = P['length']; assert self.length > 0.0
+        self.verbose = verbose
+        self.last_time = 0.0        
+        self.last_update = 0.0 # For use with update_interval        
+        self.update_interval = update_interval
+        
+
+        # Store hydraulic parameters 
+        self.manning = manning
+        self.loss_exit = loss_exit
+        self.loss_entry = loss_entry
+        self.loss_bend = loss_bend
+        self.loss_special = loss_special
+        self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
+        self.blockage_topdwn = blockage_topdwn
+        self.blockage_bottup = blockage_bottup
+
+        # Store culvert routine
+        self.culvert_routine = culvert_routine
+
+        
+        # Create objects to update momentum (a bit crude at this stage)
+        xmom0 = General_forcing(domain, 'xmomentum',
+                                polygon=P['exchange_polygon0'])
+
+        xmom1 = General_forcing(domain, 'xmomentum',
+                                polygon=P['exchange_polygon1'])
+
+        ymom0 = General_forcing(domain, 'ymomentum',
+                                polygon=P['exchange_polygon0'])
+
+        ymom1 = General_forcing(domain, 'ymomentum',
+                                polygon=P['exchange_polygon1'])
+
+        self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
+        
+
+        # Print something
+        s = 'Culvert Effective Length = %.2f m' %(self.length)
+        log_to_file(self.log_filename, s)
+   
+        s = 'Culvert Direction is %s\n' %str(self.vector)
+        log_to_file(self.log_filename, s)        
+        
+        
+    def __call__(self, domain):
+
+        log_filename = self.log_filename
+         
+        # Time stuff
+        time = domain.get_time()
+        
+        # Short hand
+        dq = domain.quantities
+                
+
+        update = False
+        if self.update_interval is None:
+            update = True
+            delta_t = domain.timestep # Next timestep has been computed in domain.py
+        else:    
+            if time - self.last_update > self.update_interval or time == 0.0:
+                update = True
+            delta_t = self.update_interval
+            
+        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
+        if hasattr(self, 'log_filename'):            
+            log_to_file(log_filename, s)
+                
+                                
+        if update is True:
+            self.last_update = time
+                        
+            msg = 'Time did not advance'
+            if time > 0.0: assert delta_t > 0.0, msg
+
+
+            # Get average water depths at each opening        
+            openings = self.openings   # There are two Opening [0] and [1]
+            for i, opening in enumerate(openings):
+                
+                # Compute mean values of selected quantitites in the 
+                # exchange area in front of the culvert
+                     
+                stage = opening.get_quantity_values(quantity_name='stage')
+                w = mean(stage) # Average stage
+
+                # Use invert level instead of elevation if specified
+                invert_level = self.invert_levels[i]
+                if invert_level is not None:
+                    z = invert_level
+                else:
+                    elevation = opening.get_quantity_values(quantity_name='elevation')
+                    z = mean(elevation) # Average elevation
+
+                # Estimated depth above the culvert inlet
+                d = w - z  # Used for calculations involving depth
+                if d < 0.0:
+                    # This is possible since w and z are taken at different locations
+                    #msg = 'D < 0.0: %f' %d
+                    #raise msg
+                    d = 0.0
+                
+
+                # Ratio of depth to culvert height.
+                # If ratio > 1 then culvert is running full
+                if self.culvert_type == 'circle':
+                    ratio = d/self.diameter
+                else:    
+                    ratio = d/self.height  
+                opening.ratio = ratio
+                    
+                    
+                # Average measures of energy in front of this opening
+                
+                id = [self.enquiry_indices[i]]
+                stage = dq['stage'].get_values(location='centroids',
+                                               indices=id)
+                elevation = dq['elevation'].get_values(location='centroids',
+                                                       indices=id)                                               
+                xmomentum = dq['xmomentum'].get_values(location='centroids',
+                                                       indices=id)                                               
+                ymomentum = dq['xmomentum'].get_values(location='centroids',
+                                                       indices=id)                                                                                              
+                depth = stage-elevation
+                if depth > 0.0:
+                    u = xmomentum/(depth + velocity_protection/depth)
+                    v = ymomentum/(depth + velocity_protection/depth)
+                else:
+                    u = v = 0.0
+
+                    
+                opening.total_energy = 0.5*(u*u + v*v)/g + stage
+
+                # Store current average stage and depth with each opening object
+                opening.depth = d
+                opening.depth_trigger = d # Use this for now
+                opening.stage = w
+                opening.elevation = z
+                
+
+            #################  End of the FOR loop ################################################
+
+            # We now need to deal with each opening individually
+                
+            # Determine flow direction based on total energy difference
+            delta_Et = openings[0].total_energy - openings[1].total_energy
+
+            if delta_Et > 0:
+                inlet = openings[0]
+                outlet = openings[1]
+
+                inlet.momentum = self.opening_momentum[0]
+                outlet.momentum = self.opening_momentum[1]
+
+            else:
+                inlet = openings[1]
+                outlet = openings[0]
+
+                inlet.momentum = self.opening_momentum[1]
+                outlet.momentum = self.opening_momentum[0]
+                
+                delta_Et = -delta_Et
+
+            self.inlet = inlet
+            self.outlet = outlet
+                
+            msg = 'Total energy difference is negative'
+            assert delta_Et > 0.0, msg
+
+            delta_h = inlet.stage - outlet.stage
+            delta_z = inlet.elevation - outlet.elevation
+            culvert_slope = (delta_z/self.length)
+
+            if culvert_slope < 0.0:
+                # Adverse gradient - flow is running uphill
+                # Flow will be purely controlled by uphill outlet face
+                if self.verbose is True:
+                    log.critical('WARNING: Flow is running uphill. Watch Out! '
+                                 'inlet.elevation=%s, outlet.elevation%s'
+                                 % (str(inlet.elevation), str(outlet.elevation)))
+
+
+            s = 'Delta total energy = %.3f' %(delta_Et)
+            log_to_file(log_filename, s)
+
+
+            # Calculate discharge for one barrel and set inlet.rate and outlet.rate
+            Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
+        
+            # Adjust discharge for multiple barrels
+            Q *= self.number_of_barrels
+
+            # Compute barrel momentum
+            barrel_momentum = barrel_velocity*culvert_outlet_depth
+                    
+            s = 'Barrel velocity = %f' %barrel_velocity
+            log_to_file(log_filename, s)
+
+            # Compute momentum vector at outlet
+            outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
+                
+            s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
+            log_to_file(log_filename, s)
+
+            # Log timeseries to file
+            fid = open(self.timeseries_filename, 'a')        
+            fid.write('%f, %f, %f, %f, %f\n'\
+                          %(time, 
+                            openings[0].total_energy,
+                            openings[1].total_energy,
+                            barrel_velocity,
+                            Q))
+            fid.close()
+
+            # Update momentum        
+            if delta_t > 0.0:
+                xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
+                xmomentum_rate /= delta_t
+                    
+                ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
+                ymomentum_rate /= delta_t
+                        
+                s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
+                log_to_file(log_filename, s)                    
+            else:
+                xmomentum_rate = ymomentum_rate = 0.0
+
+
+            # Set momentum rates for outlet jet
+            outlet.momentum[0].rate = xmomentum_rate
+            outlet.momentum[1].rate = ymomentum_rate
+
+            # Remember this value for next step (IMPORTANT)
+            outlet.momentum[0].value = outlet_mom_x
+            outlet.momentum[1].value = outlet_mom_y                    
+
+            if int(domain.time*100) % 100 == 0:
+                s = 'T=%.5f, Culvert Discharge = %.3f f'\
+                    %(time, Q)
+                s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
+                     %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
+                s += ' Momentum rate: (%.4f, %.4f)'\
+                     %(xmomentum_rate, ymomentum_rate)                    
+                s+='Outlet Vel= %.3f'\
+                    %(barrel_velocity)
+                log_to_file(log_filename, s)
+            
+            # Store value of time
+            self.last_time = time
+                
+
+
+        # Execute flow term for each opening
+        # This is where Inflow objects are evaluated and update the domain
+        for opening in self.openings:
+            opening(domain)
+
+        # Execute momentum terms
+        # This is where Inflow objects are evaluated and update the domain
+        self.outlet.momentum[0](domain)
+        self.outlet.momentum[1](domain)        
+            
+
+
+Culvert_flow = Culvert_flow_general        
Index: /trunk/anuga_core/source/anuga/structures/culvert_polygons.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/culvert_polygons.py	(revision 7939)
+++ /trunk/anuga_core/source/anuga/structures/culvert_polygons.py	(revision 7939)
@@ -0,0 +1,114 @@
+"""Functions for geometries related to culvert flows
+"""
+
+# Import necessary modules
+from math import sqrt
+from anuga.geometry.polygon import inside_polygon, polygon_area
+
+import numpy as num
+
+
+def create_culvert_polygons(end_point0,
+                            end_point1, 
+                            width, height=None,
+                            enquiry_gap_factor=0.2,
+                            number_of_barrels=1):
+    """Create polygons at the end of a culvert inlet and outlet.
+    At either end two polygons will be created; one for the actual flow to pass through and one a little further away
+    for enquiring the total energy at both ends of the culvert and transferring flow.
+
+    Input (mandatory):
+        end_point0 - one end of the culvert (x,y)
+        end_point1 - other end of the culvert (x,y)        
+        width - culvert width
+
+    Input (optional):        
+        height - culvert height, defaults to width making a square culvert
+        enquiry_gap_factor - sets the distance to the enquiry point as fraction of the height
+        number_of_barrels - number of identical pipes.
+        
+    Output:
+
+        Dictionary of four polygons. The dictionary keys are:
+            'exchange_polygon0' - polygon defining the flow area at end_point0
+            'exchange_polygon1' - polygon defining the flow area at end_point1
+            'enquiry_point0' - point beyond exchange_polygon0
+            'enquiry_point1' - point beyond exchange_polygon1            
+            'vector'
+            'length'
+            'normal'
+    """    
+
+
+    # Input check
+    if height is None:
+        height = width
+
+    # Dictionary for calculated polygons
+    culvert_polygons = {}
+    
+
+    # Calculate geometry
+    x0, y0 = end_point0
+    x1, y1 = end_point1
+
+    dx = x1-x0
+    dy = y1-y0
+
+    vector = num.array([dx, dy])
+    length = sqrt(num.sum(vector**2))
+
+    # Adjust polygon width to number of barrels in this culvert
+    width *= number_of_barrels
+    
+    
+    # Unit direction vector and normal 
+    vector /= length                 # Unit vector in culvert direction
+    normal = num.array([-dy, dx])/length # Normal vector
+    
+    culvert_polygons['vector'] = vector
+    culvert_polygons['length'] = length
+    culvert_polygons['normal'] = normal    
+
+    # Short hands
+    w = 0.5*width*normal # Perpendicular vector of 1/2 width 
+    h = height*vector    # Vector of length=height in the
+                         # direction of the culvert
+    gap = (1 + enquiry_gap_factor)*h 
+                         
+
+    # Build exchange polygon and enquiry point for opening 0
+    p0 = end_point0 + w
+    p1 = end_point0 - w
+    p2 = p1 - h
+    p3 = p0 - h
+    culvert_polygons['exchange_polygon0'] = num.array([p0,p1,p2,p3])
+    culvert_polygons['enquiry_point0'] = end_point0 - gap
+    
+
+    # Build exchange polygon and enquiry point for opening 1
+    p0 = end_point1 + w
+    p1 = end_point1 - w
+    p2 = p1 + h
+    p3 = p0 + h
+    culvert_polygons['exchange_polygon1'] = num.array([p0,p1,p2,p3])
+    culvert_polygons['enquiry_point1'] = end_point1 + gap  
+
+    # Check that enquiry polygons are outside exchange polygons
+    for key1 in ['exchange_polygon0',
+                 'exchange_polygon1']:
+        polygon = culvert_polygons[key1]
+        area = polygon_area(polygon)
+        
+        msg = 'Polygon %s ' %(polygon)
+        msg += ' has area = %f' % area
+        assert area > 0.0, msg
+
+        for key2 in ['enquiry_point0', 'enquiry_point1']:
+            point = culvert_polygons[key2]
+            msg = 'Enquiry point falls inside an enquiry point.'
+            msg += 'Email Ole.Nielsen@ga.gov.au'
+            assert not inside_polygon(point, polygon), msg
+
+    # Return results
+    return culvert_polygons
Index: /trunk/anuga_core/source/anuga/structures/culvert_routines.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/culvert_routines.py	(revision 7939)
+++ /trunk/anuga_core/source/anuga/structures/culvert_routines.py	(revision 7939)
@@ -0,0 +1,362 @@
+"""Collection of culvert routines for use with Culvert_flow in culvert_class
+
+This module holds various routines to determine FLOW through CULVERTS and SIMPLE BRIDGES
+
+
+
+
+Usage:
+
+   
+
+"""
+
+#NOTE:
+# Inlet control:  Delta_total_energy > inlet_specific_energy
+# Outlet control: Delta_total_energy < inlet_specific_energy
+# where total energy is (w + 0.5*v^2/g) and
+# specific energy is (h + 0.5*v^2/g)
+
+
+from math import pi, sqrt, sin, cos
+
+
+def boyd_generalised_culvert_model(inlet_depth, 
+                                   outlet_depth,
+                                   inlet_velocity,
+                                   outlet_velocity,
+                                   inlet_specific_energy, 
+                                   delta_total_energy, 
+                                   g,
+                                   culvert_length=0.0,
+                                   culvert_width=0.0,
+                                   culvert_height=0.0,
+                                   culvert_type='box',
+                                   manning=0.0,
+                                   sum_loss=0.0,
+                                   max_velocity=10.0,
+                                   log_filename=None):
+
+    """Boyd's generalisation of the US department of transportation culvert methods
+	
+	WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER
+	FULL PIPE OR CRITICAL DEPTH ONLY
+	For Supercritical flow this is UNDERESTIMATING the Outlet Velocity
+	The obtain the CORRECT velocity requires an iteration of Depth to Establish 
+	the Normal Depth of flow in the pipe.
+	
+	It is proposed to provide this in a seperate routine called
+    boyd_generalised_culvert_model_complex
+	
+    The Boyd Method is based on methods described by the following:
+	1.
+	US Dept. Transportation Federal Highway Administration (1965) 
+	Hydraulic Chart for Selection of Highway Culverts.
+	Hydraulic Engineering Circular No. 5 US Government Printing
+	2.
+	US Dept. Transportation Federal Highway Administration (1972) 
+	Capacity charts for the Hydraulic design of highway culverts.
+	Hydraulic Engineering Circular No. 10 US Government Printing
+    These documents provide around 60 charts for various configurations of culverts and inlets.
+	
+	Note these documents have been superceded by:
+	2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5), 
+	Which combines culvert design information previously contained in Hydraulic Engineering Circulars 
+	(HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information.
+	HEC-5 provides 20 Charts
+	HEC-10 Provides an additional 36 Charts
+	HEC-13 Discusses the Design of improved more efficient inlets
+	HDS-5 Provides 60 sets of Charts
+
+	In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in
+	1987 published "Generalised Head Discharge Equations for Culverts".
+	These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations.
+	
+	It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done.
+	The additional charts cover a range of culvert shapes and inlet configurations
+	
+   
+    """
+
+    from anuga.utilities.system_tools import log_to_file
+    from anuga.config import velocity_protection
+    from anuga.utilities.numerical_tools import safe_acos as acos
+
+    local_debug ='false'
+    if inlet_depth > 0.1: #this value was 0.01:
+        if local_debug =='true':
+            log.critical('Specific E & Deltat Tot E = %s, %s'
+                         % (str(inlet_specific_energy),
+                            str(delta_total_energy)))
+            log.critical('culvert type = %s' % str(culvert_type))
+        # Water has risen above inlet
+        
+        if log_filename is not None:                            
+            s = 'Specific energy  = %f m' % inlet_specific_energy
+            log_to_file(log_filename, s)
+        
+        msg = 'Specific energy at inlet is negative'
+        assert inlet_specific_energy >= 0.0, msg
+                      
+
+        if culvert_type =='circle':
+            # Round culvert (use height as diameter)
+            diameter = culvert_height   
+
+            """
+            For a circular pipe the Boyd method reviews 3 conditions
+            1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
+            2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
+            3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe 
+            
+            For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
+            """
+            
+            # Calculate flows for inlet control            
+            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged 
+            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
+            # Note for to SUBMERGED TO OCCUR inlet_specific_energy should be > 1.2 x diameter.... Should Check !!!
+
+            if log_filename is not None:                                        
+                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
+                log_to_file(log_filename, s)
+            Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
+
+            # THE LOWEST Value will Control Calcs From here
+            # Calculate Critical Depth Based on the Adopted Flow as an Estimate
+            dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
+            dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
+            # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
+            if dcrit1/diameter  > 0.85:
+                outlet_culvert_depth = dcrit2
+            else:
+                outlet_culvert_depth = dcrit1
+            #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
+            # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
+            if outlet_culvert_depth >= diameter:
+                outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
+                flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
+                perimeter = diameter * pi
+                flow_width= diameter
+                case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
+                if local_debug == 'true':
+                    log.critical('Inlet CTRL Outlet submerged Circular '
+                                 'PIPE FULL')
+            else:
+                #alpha = acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
+                alpha = acos(1-2*outlet_culvert_depth/diameter)*2
+                #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
+                flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
+                flow_width= diameter*sin(alpha/2.0)
+                perimeter = alpha*diameter/2.0
+                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
+                if local_debug =='true':
+                    log.critical('INLET CTRL Culvert is open channel flow '
+                                 'we will for now assume critical depth')
+                    log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
+                                 % (str(Q), str(outlet_culvert_depth),
+                                    str(alpha)))
+            if delta_total_energy < inlet_specific_energy:  #  OUTLET CONTROL !!!!
+                # Calculate flows for outlet control
+                
+                # Determine the depth at the outlet relative to the depth of flow in the Culvert
+                if outlet_depth > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
+                    outlet_culvert_depth=diameter
+                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
+                    perimeter = diameter * pi
+                    flow_width= diameter
+                    case = 'Outlet submerged'
+                    if local_debug =='true':
+                        log.critical('Outlet submerged')
+                else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
+                    # IF  outlet_depth < diameter
+                    dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
+                    dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
+                    if dcrit1/diameter >0.85:
+                        outlet_culvert_depth= dcrit2
+                    else:
+                        outlet_culvert_depth = dcrit1
+                    if outlet_culvert_depth > diameter:
+                        outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
+                        flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
+                        perimeter = diameter * pi
+                        flow_width= diameter
+                        case = 'Outlet unsubmerged PIPE FULL'
+                        if local_debug =='true':
+                            log.critical('Outlet unsubmerged PIPE FULL')
+                    else:
+                        alpha = acos(1-2*outlet_culvert_depth/diameter)*2
+                        flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
+                        flow_width= diameter*sin(alpha/2.0)
+                        perimeter = alpha*diameter/2.0
+                        case = 'Outlet is open channel flow we will for now assume critical depth'
+                        if local_debug == 'true':
+                            log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
+                                         % (str(Q), str(outlet_culvert_depth),
+                                            str(alpha)))
+                            log.critical('Outlet is open channel flow we '
+                                         'will for now assume critical depth')
+            if local_debug == 'true':
+                log.critical('FLOW AREA = %s' % str(flow_area))
+                log.critical('PERIMETER = %s' % str(perimeter))
+                log.critical('Q Interim = %s' % str(Q))
+            hyd_rad = flow_area/perimeter
+
+            if log_filename is not None:
+                s = 'hydraulic radius at outlet = %f' %hyd_rad
+                log_to_file(log_filename, s)
+
+            # Outlet control velocity using tail water
+            if local_debug =='true':
+                log.critical('GOT IT ALL CALCULATING Velocity')
+                log.critical('HydRad = %s' % str(hyd_rad))
+            culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
+            Q_outlet_tailwater = flow_area * culvert_velocity
+            if local_debug =='true':
+                log.critical('VELOCITY = %s' % str(culvert_velocity))
+                log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
+            if log_filename is not None:                
+                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
+                log_to_file(log_filename, s)
+            Q = min(Q, Q_outlet_tailwater)
+            if local_debug =='true':
+                log.critical('%s,%.3f,%.3f'
+                             % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
+                log.critical('%s,%.3f,%.3f,%.3f'
+                             % ('Q and Velocity and Depth=', Q,
+                                culvert_velocity, outlet_culvert_depth))
+
+        else:
+            pass
+                #FIXME(Ole): What about inlet control?
+        # ====  END OF CODE BLOCK FOR "IF" CIRCULAR PIPE 
+
+        #else....
+        if culvert_type == 'box':
+            if local_debug == 'true':
+                log.critical('BOX CULVERT')
+            # Box culvert (rectangle or square)   ========================================================================================================================
+
+            # Calculate flows for inlet control
+            height = culvert_height
+            width = culvert_width
+            flow_width=culvert_width            
+            
+            Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
+            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*inlet_specific_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
+
+            if log_filename is not None:                            
+                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
+                log_to_file(log_filename, s)
+
+            # FIXME(Ole): Are these functions really for inlet control?    
+            if Q_inlet_unsubmerged < Q_inlet_submerged:
+                Q = Q_inlet_unsubmerged
+                dcrit = (Q**2/g/width**2)**0.333333
+                if dcrit > height:
+                    dcrit = height
+                flow_area = width*dcrit
+                outlet_culvert_depth = dcrit
+                case = 'Inlet unsubmerged Box Acts as Weir'
+            else:    
+                Q = Q_inlet_submerged
+                flow_area = width*height
+                outlet_culvert_depth = height
+                case = 'Inlet submerged Box Acts as Orifice'                    
+
+            dcrit = (Q**2/g/width**2)**0.333333
+
+            outlet_culvert_depth = dcrit
+            if outlet_culvert_depth > height:
+                outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
+                flow_area = width*height  # Cross sectional area of flow in the culvert
+                perimeter = 2*(width+height)
+                case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
+            else:
+                flow_area = width * outlet_culvert_depth
+                perimeter = width+2*outlet_culvert_depth
+                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
+	
+            if delta_total_energy < inlet_specific_energy:
+                # Calculate flows for outlet control
+                
+                # Determine the depth at the outlet relative to the depth of flow in the Culvert
+                if outlet_depth > height:        # The Outlet is Submerged
+                    outlet_culvert_depth=height
+                    flow_area=width*height       # Cross sectional area of flow in the culvert
+                    perimeter=2.0*(width+height)
+                    case = 'Outlet submerged'
+                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
+                    dcrit = (Q**2/g/width**2)**0.333333
+                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
+                    if outlet_culvert_depth > height:
+                        outlet_culvert_depth=height
+                        flow_area=width*height
+                        perimeter=2.0*(width+height)
+                        case = 'Outlet is Flowing Full'
+                    else:
+                        flow_area=width*outlet_culvert_depth
+                        perimeter=(width+2.0*outlet_culvert_depth)
+                        case = 'Outlet is open channel flow'
+
+                hyd_rad = flow_area/perimeter
+
+                if log_filename is not None:                                
+                    s = 'hydraulic radius at outlet = %f' % hyd_rad
+                    log_to_file(log_filename, s)
+
+                # Outlet control velocity using tail water
+                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
+                Q_outlet_tailwater = flow_area * culvert_velocity
+
+                if log_filename is not None:                            
+                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
+                    log_to_file(log_filename, s)
+                Q = min(Q, Q_outlet_tailwater)
+            else:
+                pass
+                #FIXME(Ole): What about inlet control?
+        # ====  END OF CODE BLOCK FOR "IF" BOX 
+
+                
+        # Common code for circle and box geometries ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+        if log_filename is not None:
+            log_to_file(log_filename, 'Case: "%s"' % case)
+            
+        if log_filename is not None:                        
+            s = 'Flow Rate Control = %f' % Q
+            log_to_file(log_filename, s)
+
+        culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
+        if local_debug =='true':
+            log.critical('FLOW AREA = %s' % str(flow_area))
+            log.critical('PERIMETER = %s' % str(perimeter))
+            log.critical('Q final = %s' % str(Q))
+            log.critical('FROUDE = %s' % str(culv_froude))
+        if log_filename is not None:                            
+            s = 'Froude in Culvert = %f' % culv_froude
+            log_to_file(log_filename, s)
+
+        # Determine momentum at the outlet 
+        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
+
+    # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow  
+
+    else: # inlet_depth < 0.01:
+        Q = barrel_velocity = outlet_culvert_depth = 0.0
+
+    # Temporary flow limit
+    if barrel_velocity > max_velocity:
+        if log_filename is not None:                            
+            s = 'Barrel velocity was reduced from = %f m/s to %f m/s' % (barrel_velocity, max_velocity)
+            log_to_file(log_filename, s)
+        
+        barrel_velocity = max_velocity
+        Q = flow_area * barrel_velocity
+        
+        
+
+        
+    return Q, barrel_velocity, outlet_culvert_depth 
+
+
Index: /trunk/anuga_core/source/anuga/structures/test_all.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/test_all.py	(revision 7939)
+++ /trunk/anuga_core/source/anuga/structures/test_all.py	(revision 7939)
@@ -0,0 +1,83 @@
+"""Regression testing framework
+This module will search for scripts in the same directory named
+test_*.py.  Each such script should be a test suite that tests a
+module through PyUnit. This script will aggregate all
+found test suites into one big test suite and run them all at once.
+"""
+
+# Author: Mark Pilgrim
+# Modified by Ole Nielsen
+
+import unittest
+import os
+
+
+#List files that should be excluded from the testing process.
+#E.g. if they are known to fail and under development
+exclude = [] 
+
+
+def get_test_files(path):
+
+    import sys
+
+    files = os.listdir(path)
+
+    #Check sub directories
+    test_files = []
+    for file in files:
+        if os.path.isdir(file):
+            sys.path.append(file)
+            #print 'Recursing into', file
+            test_files += get_test_files(path + os.sep + file)
+        elif file[:5] == 'test_' and file[-2:] == 'py':
+            #print 'Appending', file
+            test_files.append(file)
+        else:
+            pass
+    return test_files
+
+
+
+def regressionTest():
+    import sys, os, re, unittest
+    path = os.path.split(sys.argv[0])[0] or os.getcwd()
+
+
+    files = get_test_files(path)
+
+
+
+    #test = re.compile('^test_[\w]*.py$', re.IGNORECASE)
+    #files = filter(test.search, files)
+
+
+    try:
+        files.remove(__file__)  #Remove self from list (Ver 2.3. or later)
+    except:
+        files.remove('test_all.py')
+
+    print 'Testing:'
+    for file in files:
+        print '  ' + file
+
+    if globals().has_key('exclude'):
+        for file in exclude:
+            files.remove(file)
+            print 'WARNING: File '+ file + ' excluded from testing'
+
+
+    filenameToModuleName = lambda f: os.path.splitext(f)[0]
+    #print "files",files 
+    moduleNames = map(filenameToModuleName, files)
+    #print "moduleNames",moduleNames
+    modules = map(__import__, moduleNames)
+    load = unittest.defaultTestLoader.loadTestsFromModule
+    return unittest.TestSuite(map(load, modules))
+
+################################################################################
+
+if __name__ == '__main__':   
+    suite = regressionTest()
+    runner = unittest.TextTestRunner() #verbosity=2)
+    runner.run(suite)
Index: /trunk/anuga_core/source/anuga/structures/test_culvert_class.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/test_culvert_class.py	(revision 7939)
+++ /trunk/anuga_core/source/anuga/structures/test_culvert_class.py	(revision 7939)
@@ -0,0 +1,810 @@
+#!/usr/bin/env python
+
+
+import unittest
+import os.path
+import sys
+
+from anuga.utilities.system_tools import get_pathname_from_package
+from anuga.geometry.polygon_function import Polygon_function
+        
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+
+import anuga
+
+from anuga.culvert_flows.culvert_class import Culvert_flow, \
+                            Culvert_flow_rating, Culvert_flow_energy
+from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
+     
+from math import pi, pow, sqrt
+
+import numpy as num
+
+
+# Helper functions
+def run_culvert_flow_problem(depth):
+    """Run flow with culvert given depth
+    """
+
+    length = 40.
+    width = 5.
+
+    dx = dy = 1           # Resolution: Length of subdivisions on both axes
+
+    points, vertices, boundary = rectangular_cross(int(length/dx),
+                                                   int(width/dy),
+                                                   len1=length, 
+                                                   len2=width)
+    domain = anuga.Domain(points, vertices, boundary)   
+    domain.set_name('Test_culvert_shallow') # Output name
+    domain.set_default_order(2)
+
+
+    #----------------------------------------------------------------------
+    # Setup initial conditions
+    #----------------------------------------------------------------------
+
+    def topography(x, y):
+        """Set up a weir
+        
+        A culvert will connect either side
+        """
+        # General Slope of Topography
+        z=-x/1000
+        
+        N = len(x)
+        for i in range(N):
+
+           # Sloping Embankment Across Channel
+            if 5.0 < x[i] < 10.1:
+                # Cut Out Segment for Culvert face                
+                if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+                   z[i]=z[i]
+                else:
+                   z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+            if 10.0 < x[i] < 12.1:
+               z[i] +=  2.5                    # Flat Crest of Embankment
+            if 12.0 < x[i] < 14.5:
+                # Cut Out Segment for Culvert face                
+                if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+                   z[i]=z[i]
+                else:
+                   z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+            		   
+    		
+        return z
+
+
+    domain.set_quantity('elevation', topography) 
+    domain.set_quantity('friction', 0.01)         # Constant friction 
+    domain.set_quantity('stage',
+                        expression='elevation + %f' % depth) # Shallow initial condition
+
+    # Boyd culvert
+    culvert = Culvert_flow(domain,
+                           label='Culvert No. 1',
+                           description='This culvert is a test unit 1.2m Wide by 0.75m High',   
+                           end_point0=[9.0, 2.5], 
+                           end_point1=[13.0, 2.5],
+                           width=1.20, height=0.75,
+                           culvert_routine=boyd_generalised_culvert_model,
+                           number_of_barrels=1,
+                           update_interval=2,
+                           verbose=False)
+    
+
+    domain.forcing_terms.append(culvert)
+    
+
+    #-----------------------------------------------------------------------
+    # Setup boundary conditions
+    #-----------------------------------------------------------------------
+
+    # Inflow based on Flow Depth and Approaching Momentum
+
+    Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+    domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+    
+    #-----------------------------------------------------------------------
+    # Evolve system through time
+    #-----------------------------------------------------------------------
+
+    #print 'depth', depth
+    ref_volume = domain.get_quantity('stage').get_integral()
+    for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
+        new_volume = domain.get_quantity('stage').get_integral()
+        
+        msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
+               % (new_volume, ref_volume))
+        assert num.allclose(new_volume, ref_volume), msg        
+        
+
+    os.remove('Test_culvert_shallow.sww')
+
+class Test_Culvert(unittest.TestCase):
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+
+    def test_that_culvert_runs_rating(self):
+        """test_that_culvert_runs_rating
+        
+        This test exercises the culvert and checks values outside rating curve
+        are dealt with       
+        """
+
+        path = get_pathname_from_package('anuga.culvert_flows')    
+        
+        length = 40.
+        width = 5.
+
+        dx = dy = 1           # Resolution: Length of subdivisions on both axes
+
+        points, vertices, boundary = rectangular_cross(int(length/dx),
+                                                       int(width/dy),
+                                                       len1=length, 
+                                                       len2=width)
+        domain = anuga.Domain(points, vertices, boundary)   
+        domain.set_name('Test_culvert')                 # Output name
+        domain.set_default_order(2)
+
+
+        #----------------------------------------------------------------------
+        # Setup initial conditions
+        #----------------------------------------------------------------------
+
+        def topography(x, y):
+            """Set up a weir
+            
+            A culvert will connect either side
+            """
+            # General Slope of Topography
+            z=-x/1000
+            
+            N = len(x)
+            for i in range(N):
+
+               # Sloping Embankment Across Channel
+                if 5.0 < x[i] < 10.1:
+                    # Cut Out Segment for Culvert face                
+                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+                       z[i]=z[i]
+                    else:
+                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+                if 10.0 < x[i] < 12.1:
+                   z[i] +=  2.5                    # Flat Crest of Embankment
+                if 12.0 < x[i] < 14.5:
+                    # Cut Out Segment for Culvert face                
+                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+                       z[i]=z[i]
+                    else:
+                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                		   
+        		
+            return z
+
+
+        domain.set_quantity('elevation', topography) 
+        domain.set_quantity('friction', 0.01)         # Constant friction 
+        domain.set_quantity('stage',
+                            expression='elevation')   # Dry initial condition
+
+        filename=os.path.join(path, 'example_rating_curve.csv')
+        culvert = Culvert_flow(domain,
+                               culvert_description_filename=filename,        
+                               end_point0=[9.0, 2.5], 
+                               end_point1=[13.0, 2.5],
+                               width=1.00,
+                               use_velocity_head=True,
+                               verbose=False)
+
+        domain.forcing_terms.append(culvert)
+        
+
+        #-----------------------------------------------------------------------
+        # Setup boundary conditions
+        #-----------------------------------------------------------------------
+
+        # Inflow based on Flow Depth and Approaching Momentum
+        Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])
+        Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+        Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # Outflow
+        
+        # Upstream and downstream conditions that will exceed the rating curve
+        # I.e produce delta_h outside the range [0, 10] specified in the the 
+        # file example_rating_curve.csv
+        Btus = anuga.Time_boundary(domain, \
+                    lambda t: [100*num.sin(2*pi*(t-4)/10), 0.0, 0.0])
+        Btds = anuga.Time_boundary(domain, \
+                    lambda t: [-5*(num.cos(2*pi*(t-4)/20)), 0.0, 0.0])
+        domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
+
+
+        #-----------------------------------------------------------------------
+        # Evolve system through time
+        #-----------------------------------------------------------------------
+
+        min_delta_w = sys.maxint 
+        max_delta_w = -min_delta_w
+        for t in domain.evolve(yieldstep = 1, finaltime = 25):
+            delta_w = culvert.inlet.stage - culvert.outlet.stage
+            
+            if delta_w > max_delta_w: max_delta_w = delta_w
+            if delta_w < min_delta_w: min_delta_w = delta_w            
+            
+            pass
+
+        # Check that extreme values in rating curve have been exceeded
+        # so that we know that condition has been exercised
+        assert min_delta_w < 0
+        assert max_delta_w > 10        
+
+
+        os.remove('Test_culvert.sww')
+
+    def test_that_culvert_dry_bed_rating_does_not_produce_flow(self):
+        """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
+        
+        Test that culvert on a sloping dry bed doesn't produce flows
+        although there will be a 'pressure' head due to delta_w > 0
+
+        This one is using the rating curve variant
+        """
+
+        path = get_pathname_from_package('anuga.culvert_flows')    
+        
+        length = 40.
+        width = 5.
+
+        dx = dy = 1           # Resolution: Length of subdivisions on both axes
+
+        points, vertices, boundary = rectangular_cross(int(length/dx),
+                                                       int(width/dy),
+                                                       len1=length, 
+                                                       len2=width)
+        domain = anuga.Domain(points, vertices, boundary)   
+        domain.set_name('Test_culvert_dry') # Output name
+        domain.set_default_order(2)
+
+
+        #----------------------------------------------------------------------
+        # Setup initial conditions
+        #----------------------------------------------------------------------
+
+        def topography(x, y):
+            """Set up a weir
+            
+            A culvert will connect either side
+            """
+            # General Slope of Topography
+            z=-x/1000
+            
+            N = len(x)
+            for i in range(N):
+
+               # Sloping Embankment Across Channel
+                if 5.0 < x[i] < 10.1:
+                    # Cut Out Segment for Culvert face                
+                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+                       z[i]=z[i]
+                    else:
+                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+                if 10.0 < x[i] < 12.1:
+                   z[i] +=  2.5                    # Flat Crest of Embankment
+                if 12.0 < x[i] < 14.5:
+                    # Cut Out Segment for Culvert face                
+                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+                       z[i]=z[i]
+                    else:
+                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                		   
+        		
+            return z
+
+
+        domain.set_quantity('elevation', topography) 
+        domain.set_quantity('friction', 0.01)         # Constant friction 
+        domain.set_quantity('stage',
+                            expression='elevation')   # Dry initial condition
+
+
+        filename = os.path.join(path, 'example_rating_curve.csv')
+        culvert = Culvert_flow(domain,
+                               culvert_description_filename=filename,
+                               end_point0=[9.0, 2.5], 
+                               end_point1=[13.0, 2.5],
+                               height=0.75,
+                               verbose=False)
+
+        domain.forcing_terms.append(culvert)
+        
+
+        #-----------------------------------------------------------------------
+        # Setup boundary conditions
+        #-----------------------------------------------------------------------
+
+        # Inflow based on Flow Depth and Approaching Momentum
+
+        Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+        #-----------------------------------------------------------------------
+        # Evolve system through time
+        #-----------------------------------------------------------------------
+
+        ref_volume = domain.get_quantity('stage').get_integral()
+        for t in domain.evolve(yieldstep = 1, finaltime = 25):
+            new_volume = domain.get_quantity('stage').get_integral()
+            
+            msg = 'Total volume has changed'
+            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
+            pass
+    
+
+    
+        os.remove('Test_culvert_dry.sww')
+        
+    def test_that_culvert_flows_conserves_volume(self):
+        """test_that_culvert_flows_conserves_volume
+
+        Test that culvert on a sloping dry bed limits flows when very little water
+        is present at inlet.
+
+        Uses helper function: run_culvert_flow_problem(depth):
+
+        """
+
+        # Try this for a range of depths
+        for depth in [0.1, 0.2, 0.5, 1.0]:
+            run_culvert_flow_problem(depth)
+
+
+    def OBSOLETE_XXXtest_that_culvert_rating_limits_flow_in_shallow_inlet_condition(self):
+        """test_that_culvert_rating_limits_flow_in_shallow_inlet_condition
+        
+        Test that culvert on a sloping dry bed limits flows when very little water
+        is present at inlet
+
+        This one is using the rating curve variant
+        """
+
+        
+
+        path = get_pathname_from_package('anuga.culvert_flows')    
+        
+        length = 40.
+        width = 5.
+
+        dx = dy = 1           # Resolution: Length of subdivisions on both axes
+
+        points, vertices, boundary = rectangular_cross(int(length/dx),
+                                                       int(width/dy),
+                                                       len1=length, 
+                                                       len2=width)
+        domain = anuga.Domain(points, vertices, boundary)   
+        domain.set_name('Test_culvert_shallow') # Output name
+        domain.set_default_order(2)
+
+
+        #----------------------------------------------------------------------
+        # Setup initial conditions
+        #----------------------------------------------------------------------
+
+        def topography(x, y):
+            """Set up a weir
+            
+            A culvert will connect either side
+            """
+            # General Slope of Topography
+            z=-x/1000
+            
+            N = len(x)
+            for i in range(N):
+
+               # Sloping Embankment Across Channel
+                if 5.0 < x[i] < 10.1:
+                    # Cut Out Segment for Culvert face                
+                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+                       z[i]=z[i]
+                    else:
+                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+                if 10.0 < x[i] < 12.1:
+                   z[i] +=  2.5                    # Flat Crest of Embankment
+                if 12.0 < x[i] < 14.5:
+                    # Cut Out Segment for Culvert face                
+                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+                       z[i]=z[i]
+                    else:
+                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                		   
+        		
+            return z
+
+
+        domain.set_quantity('elevation', topography) 
+        domain.set_quantity('friction', 0.01)         # Constant friction 
+        domain.set_quantity('stage',
+                            expression='elevation + 0.1') # Shallow initial condition
+
+        # Boyd culvert
+        culvert = Culvert_flow(domain,
+                               label='Culvert No. 1',
+                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
+                               end_point0=[9.0, 2.5], 
+                               end_point1=[13.0, 2.5],
+                               width=1.20, height=0.75,
+                               culvert_routine=boyd_generalised_culvert_model,
+                               number_of_barrels=1,
+                               update_interval=2,
+                               verbose=False)
+        
+        # Rating curve
+        #filename = os.path.join(path, 'example_rating_curve.csv')
+        #culvert = Culvert_flow(domain,
+        #                       culvert_description_filename=filename,
+        #                       end_point0=[9.0, 2.5], 
+        #                       end_point1=[13.0, 2.5],
+        #                       trigger_depth=0.01,
+        #                       verbose=False)
+
+        domain.forcing_terms.append(culvert)
+        
+
+        #-----------------------------------------------------------------------
+        # Setup boundary conditions
+        #-----------------------------------------------------------------------
+
+        # Inflow based on Flow Depth and Approaching Momentum
+
+        Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+        
+        #-----------------------------------------------------------------------
+        # Evolve system through time
+        #-----------------------------------------------------------------------
+
+        print 'depth', 0.1
+        ref_volume = domain.get_quantity('stage').get_integral()
+        for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
+            new_volume = domain.get_quantity('stage').get_integral()
+
+            msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
+                   % (new_volume, ref_volume))
+            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg        
+        
+        
+        return
+        # Now try this again for a depth of 10 cm and for a range of other depths
+        for depth in [0.1, 0.2, 0.5, 1.0]:
+            print 'depth', depth
+            domain.set_time(0.0)
+            
+            domain.set_quantity('elevation', topography) 
+            domain.set_quantity('friction', 0.01)         # Constant friction 
+            domain.set_quantity('stage',
+                                expression='elevation + %f' % depth)
+            
+        
+            ref_volume = domain.get_quantity('stage').get_integral()
+            for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
+                new_volume = domain.get_quantity('stage').get_integral()
+            
+                msg = 'Total volume has changed: Is %.8f m^3 should have been %.8f m^3'\
+                    % (new_volume, ref_volume)
+
+                assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
+    
+    
+    
+    def test_that_culvert_dry_bed_boyd_does_not_produce_flow(self):
+        """test_that_culvert_in_dry_bed_boyd_does_not_produce_flow(self):
+        
+        Test that culvert on a sloping dry bed doesn't produce flows
+        although there will be a 'pressure' head due to delta_w > 0
+
+        This one is using the 'Boyd' variant        
+        """
+
+        path = get_pathname_from_package('anuga.culvert_flows')    
+        
+        length = 40.
+        width = 5.
+
+        dx = dy = 1           # Resolution: Length of subdivisions on both axes
+
+        points, vertices, boundary = rectangular_cross(int(length/dx),
+                                                       int(width/dy),
+                                                       len1=length, 
+                                                       len2=width)
+        domain = anuga.Domain(points, vertices, boundary)   
+        domain.set_name('Test_culvert_dry') # Output name
+        domain.set_default_order(2)
+
+
+        #----------------------------------------------------------------------
+        # Setup initial conditions
+        #----------------------------------------------------------------------
+
+        def topography(x, y):
+            """Set up a weir
+            
+            A culvert will connect either side
+            """
+            # General Slope of Topography
+            z=-x/1000
+            
+            N = len(x)
+            for i in range(N):
+
+               # Sloping Embankment Across Channel
+                if 5.0 < x[i] < 10.1:
+                    # Cut Out Segment for Culvert face                
+                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+                       z[i]=z[i]
+                    else:
+                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+                if 10.0 < x[i] < 12.1:
+                   z[i] +=  2.5                    # Flat Crest of Embankment
+                if 12.0 < x[i] < 14.5:
+                    # Cut Out Segment for Culvert face                
+                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+                       z[i]=z[i]
+                    else:
+                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                		   
+        		
+            return z
+
+
+        domain.set_quantity('elevation', topography) 
+        domain.set_quantity('friction', 0.01)         # Constant friction 
+        domain.set_quantity('stage',
+                            expression='elevation')   # Dry initial condition
+
+
+        filename = os.path.join(path, 'example_rating_curve.csv')
+
+
+        culvert = Culvert_flow(domain,
+                               label='Culvert No. 1',
+                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
+                               end_point0=[9.0, 2.5], 
+                               end_point1=[13.0, 2.5],
+                               width=1.20, height=0.75,
+                               culvert_routine=boyd_generalised_culvert_model,
+                               number_of_barrels=1,
+                               update_interval=2,
+                               verbose=False)
+        
+        domain.forcing_terms.append(culvert)
+        
+
+        #-----------------------------------------------------------------------
+        # Setup boundary conditions
+        #-----------------------------------------------------------------------
+
+        # Inflow based on Flow Depth and Approaching Momentum
+
+        Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+        #-----------------------------------------------------------------------
+        # Evolve system through time
+        #-----------------------------------------------------------------------
+
+        ref_volume = domain.get_quantity('stage').get_integral()
+        for t in domain.evolve(yieldstep = 1, finaltime = 25):
+            
+            new_volume = domain.get_quantity('stage').get_integral()
+
+            msg = 'Total volume has changed'
+            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
+            pass
+    
+
+    
+    
+
+    def test_predicted_boyd_flow(self):
+        """test_predicted_boyd_flow
+        
+        Test that flows predicted by the boyd method are consistent with what what
+        is calculated in engineering codes.
+        The data was supplied by Petar Milevski
+        """
+
+        # FIXME(Ole) this is nowhere near finished
+        path = get_pathname_from_package('anuga.culvert_flows')    
+        
+        length = 12.
+        width = 5.
+
+        dx = dy = 0.5           # Resolution: Length of subdivisions on both axes
+
+        points, vertices, boundary = rectangular_cross(int(length/dx),
+                                                       int(width/dy),
+                                                       len1=length, 
+                                                       len2=width)
+        domain = anuga.Domain(points, vertices, boundary)   
+        
+        domain.set_name('test_culvert')                 # Output name
+        domain.set_default_order(2)
+
+
+        #----------------------------------------------------------------------
+        # Setup initial conditions
+        #----------------------------------------------------------------------
+
+        def topography(x, y):
+            # General Slope of Topography
+            z=-x/10
+            
+            return z
+
+
+        domain.set_quantity('elevation', topography) 
+        domain.set_quantity('friction', 0.01)         # Constant friction 
+        domain.set_quantity('stage', expression='elevation')
+
+
+        Q0 = domain.get_quantity('stage')
+        Q1 = Quantity(domain)
+        
+        # Add depths to stage        
+        head_water_depth = 0.169
+        tail_water_depth = 0.089
+        
+        inlet_poly = [[0,0], [6,0], [6,5], [0,5]]
+        outlet_poly = [[6,0], [12,0], [12,5], [6,5]]        
+        
+        Q1.set_values(Polygon_function([(inlet_poly, head_water_depth),
+                                        (outlet_poly, tail_water_depth)]))
+        
+        domain.set_quantity('stage', Q0 + Q1)
+
+
+
+        culvert = Culvert_flow(domain,
+                               label='Test culvert',
+                               description='4 m test culvert',   
+                               end_point0=[4.0, 2.5], 
+                               end_point1=[8.0, 2.5],
+                               width=1.20, 
+                               height=0.75,
+                               culvert_routine=boyd_generalised_culvert_model,        
+                               number_of_barrels=1,
+                               verbose=False)
+                               
+        
+        domain.forcing_terms.append(culvert)
+        
+        # Call
+        culvert(domain)
+    
+
+    
+
+    def test_momentum_jet(self):
+        """test_momentum_jet
+        
+        Test that culvert_class can accept keyword use_momentum_jet
+        This does not yet imply that the values have been tested. FIXME
+        """
+
+
+        length = 40.
+        width = 5.
+
+        dx = dy = 1           # Resolution: Length of subdivisions on both axes
+
+        points, vertices, boundary = rectangular_cross(int(length/dx),
+                                                       int(width/dy),
+                                                       len1=length, 
+                                                       len2=width)
+        domain = anuga.Domain(points, vertices, boundary)   
+        domain.set_name('Test_culvert_shallow') # Output name
+        domain.set_default_order(2)
+
+
+        #----------------------------------------------------------------------
+        # Setup initial conditions
+        #----------------------------------------------------------------------
+
+        def topography(x, y):
+            """Set up a weir
+            
+            A culvert will connect either side
+            """
+            # General Slope of Topography
+            z=-x/1000
+            
+            N = len(x)
+            for i in range(N):
+
+               # Sloping Embankment Across Channel
+                if 5.0 < x[i] < 10.1:
+                    # Cut Out Segment for Culvert face                
+                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+                       z[i]=z[i]
+                    else:
+                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+                if 10.0 < x[i] < 12.1:
+                   z[i] +=  2.5                    # Flat Crest of Embankment
+                if 12.0 < x[i] < 14.5:
+                    # Cut Out Segment for Culvert face                
+                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+                       z[i]=z[i]
+                    else:
+                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                		   
+        		
+            return z
+
+
+        domain.set_quantity('elevation', topography) 
+        domain.set_quantity('friction', 0.01)         # Constant friction 
+        domain.set_quantity('stage',
+                            expression='elevation + 1.0') # Shallow initial condition
+
+        # Boyd culvert
+        culvert = Culvert_flow(domain,
+                               label='Culvert No. 1',
+                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
+                               end_point0=[9.0, 2.5], 
+                               end_point1=[13.0, 2.5],
+                               width=1.20, height=0.75,
+                               culvert_routine=boyd_generalised_culvert_model,
+                               number_of_barrels=1,
+                               use_momentum_jet=True,
+                               update_interval=2,
+                               verbose=False)
+    
+
+        domain.forcing_terms.append(culvert)
+
+        
+        # Call
+        culvert(domain)
+    
+
+        #-----------------------------------------------------------------------
+        # Setup boundary conditions
+        #-----------------------------------------------------------------------
+
+
+        Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+        #-----------------------------------------------------------------------
+        # Evolve system through time
+        #-----------------------------------------------------------------------
+
+        ref_volume = domain.get_quantity('stage').get_integral()
+        for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
+            new_volume = domain.get_quantity('stage').get_integral()
+            
+            msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
+                   % (new_volume, ref_volume))
+            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg        
+        
+
+        
+
+               
+#-------------------------------------------------------------
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Culvert, 'test')
+    runner = unittest.TextTestRunner() #verbosity=2)
+    runner.run(suite)
+        
Index: /trunk/anuga_core/source/anuga/structures/test_culvert_polygons.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/test_culvert_polygons.py	(revision 7939)
+++ /trunk/anuga_core/source/anuga/structures/test_culvert_polygons.py	(revision 7939)
@@ -0,0 +1,84 @@
+#!/usr/bin/env python
+
+from culvert_polygons import *
+
+import unittest
+import os.path
+
+from anuga.geometry.polygon import inside_polygon, polygon_area
+
+
+class Test_poly(unittest.TestCase):
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+
+    def test_1(self):
+        
+        end_point0=[307138.813,6193474]
+        end_point1=[307150.563,6193469]
+        width=3 
+        height=3
+        number_of_barrels=1
+
+        P = create_culvert_polygons(end_point0,
+                                    end_point1,
+                                    width=width,   
+                                    height=height,
+                                    number_of_barrels=number_of_barrels)
+        
+        # Compute area and check that it is greater than 0
+        for key1 in ['exchange_polygon0',
+                     'exchange_polygon1']:
+            polygon = P[key1]
+            area = polygon_area(polygon)
+            
+            msg = 'Polygon %s ' %(polygon)
+            msg += ' has area = %f' % area
+            assert area > 0.0, msg
+
+            for key2 in ['enquiry_point0', 'enquiry_point1']:
+                point = P[key2]
+                assert not inside_polygon(point, polygon)
+        
+
+    def test_2(self):
+        #end_point0=[307138.813,6193474]
+        #end_point1=[307150.563,6193469]
+        end_point0=[10., 5.]
+        end_point1=[10., 10.]     
+        width = 1
+        height = 3.5 
+        number_of_barrels=1
+
+        P = create_culvert_polygons(end_point0,
+                                    end_point1,
+                                    width=width,   
+                                    height=height,
+                                    number_of_barrels=number_of_barrels)
+        
+        # Compute area and check that it is greater than 0
+        for key1 in ['exchange_polygon0',
+                    'exchange_polygon1']:
+            polygon = P[key1]
+            area = polygon_area(polygon)
+            
+            msg = 'Polygon %s ' % (polygon)
+            msg += ' has area = %f' % area
+            assert area > 0.0, msg
+
+            for key2 in ['enquiry_point0', 'enquiry_point1']:
+                point = P[key2]
+                assert not inside_polygon(point, polygon)                        
+
+    
+               
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_poly, 'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
+        
Index: /trunk/anuga_core/source/anuga/structures/test_culvert_routines.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/test_culvert_routines.py	(revision 7939)
+++ /trunk/anuga_core/source/anuga/structures/test_culvert_routines.py	(revision 7939)
@@ -0,0 +1,586 @@
+#!/usr/bin/env python
+
+
+import unittest
+import os.path
+import sys
+
+from anuga.utilities.system_tools import get_pathname_from_package
+from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
+import numpy as num
+
+
+class Test_culvert_routines(unittest.TestCase):
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+
+    def test_boyd_0(self):
+        """test_boyd_0
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski
+        This test is the only one that passed in late February 2009
+        """
+      
+        g=9.81
+        culvert_slope=0.1  # Downward
+
+        inlet_depth=2.0
+        outlet_depth=0.0
+        
+        inlet_velocity=0.0,
+        outlet_velocity=0.0,        
+
+        culvert_length=4.0
+        culvert_width=1.2
+        culvert_height=0.75
+
+        culvert_type='box'
+        manning=0.013
+        sum_loss=0.0
+
+        inlet_specific_energy=inlet_depth #+0.5*v**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth # + 
+        E_out = z_out+outlet_depth # +
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth, 
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        #print Q, v, d
+        assert num.allclose(Q, 3.118, rtol=1.0e-3)
+        
+
+        #assert num.allclose(v, 0.93)
+        #assert num.allclose(d, 0.0)
+        
+
+    def Xtest_boyd_00(self):
+        """test_boyd_00
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=0.1  # Downward
+
+        inlet_depth=0.2
+        outlet_depth=0.0
+
+        inlet_velocity=0.0,
+        outlet_velocity=0.0,                
+        
+        culvert_length=4.0
+        culvert_width=1.2
+        culvert_height=0.75
+
+        culvert_type='box'
+        manning=0.013
+        sum_loss=0.0
+
+        inlet_specific_energy=inlet_depth #+0.5*v**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth # + 
+        E_out = z_out+outlet_depth # +
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth, 
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        #print Q, v, d
+        assert num.allclose(Q, 0.185, rtol=1.0e-3)
+        #assert num.allclose(v, 0.93)
+        #assert num.allclose(d, 0.0)
+        
+    def Xtest_boyd_1(self):
+        """test_boyd_1
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=0.01  # Downward
+
+        inlet_depth=0.263
+        outlet_depth=0.0
+
+        culvert_length=4.0
+        culvert_width=0.75
+        culvert_height=0.75
+        
+        culvert_type='pipe'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth #+0.5*v**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
+        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth, 
+                                                 outlet_depth,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        print Q, v, d
+        assert num.allclose(Q, 0.10, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 1.13, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.15, rtol=1.0e-2) #depth at outlet used to calc v 
+        
+    def Xtest_boyd_2(self):
+        """test_boyd_2
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=0.01  # Downward
+
+        inlet_depth=1.135
+        outlet_depth=0.0
+
+        culvert_length=4.0
+        culvert_width=0.75
+        culvert_height=0.75
+        
+        culvert_type='pipe'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth #+0.5*v**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
+        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth, 
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        print Q, v, d
+        assert num.allclose(Q, 1.00, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 2.59, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.563, rtol=1.0e-2) #depth at outlet used to calc v  
+
+    def Xtest_boyd_3(self):
+        """test_boyd_3
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=0.01  # Downward
+
+        inlet_depth=12.747
+        outlet_depth=0.0
+
+        culvert_length=4.0
+        culvert_width=0.75
+        culvert_height=0.75
+        
+        culvert_type='pipe'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth #+0.5*v**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
+        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth, 
+                                                 outlet_depth,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        print Q, v, d
+        assert num.allclose(Q, 5.00, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 11.022, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.72, rtol=1.0e-2) #depth at outlet used to calc v
+
+    def Xtest_boyd_4(self):
+        """test_boyd_4
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=0.01  # Downward
+
+        inlet_depth=1.004
+        outlet_depth=1.00
+
+        culvert_length=4.0
+        culvert_width=0.75
+        culvert_height=0.75
+        
+        culvert_type='pipe'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth #+0.5*v**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
+        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth, 
+                                                 outlet_depth,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        print Q, v, d
+        assert num.allclose(Q, 0.10, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 0.22, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.76, rtol=1.0e-2) #depth at outlet used to calc v
+
+    def Xtest_boyd_5(self):
+        """test_boyd_5
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=0.01  # Downward
+
+        inlet_depth=1.401
+        outlet_depth=1.00
+
+        culvert_length=4.0
+        culvert_width=0.75
+        culvert_height=0.75
+        
+        culvert_type='pipe'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth #+0.5*v**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
+        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth, 
+                                                 outlet_depth,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        print Q, v, d
+        assert num.allclose(Q, 1.00, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 2.204, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.76, rtol=1.0e-2) #depth at outlet used to calc v
+
+
+    def Xtest_boyd_6(self):
+        """test_boyd_5
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=0.01  # Downward
+
+        inlet_depth=12.747
+        outlet_depth=1.00
+
+        culvert_length=4.0
+        culvert_width=0.75
+        culvert_height=0.75
+        
+        culvert_type='pipe'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth #+0.5*v**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
+        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth, 
+                                                 outlet_depth,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        print Q, v, d
+        assert num.allclose(Q, 5.00, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 11.022, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.76, rtol=1.0e-2) #depth at outlet used to calc v
+
+
+    def Xtest_boyd_7(self):
+        """test_boyd_7
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=0.1  # Downward
+
+        inlet_depth=0.303
+        outlet_depth=0.00
+
+        culvert_length=4.0
+        culvert_width=0.75
+        culvert_height=0.75
+        
+        culvert_type='pipe'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth #+0.5*v**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
+        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth, 
+                                                 outlet_depth,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        print Q, v, d
+        assert num.allclose(Q, 0.10, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 1.13, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.19, rtol=1.0e-2) #depth at outlet used to calc v
+
+
+    def Xtest_boyd_8(self):
+        """test_boyd_8
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=0.1  # Downward
+
+        inlet_depth=1.135
+        outlet_depth=0.00
+
+        culvert_length=4.0
+        culvert_width=0.75
+        culvert_height=0.75
+        
+        culvert_type='pipe'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth #+0.5*v**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
+        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth, 
+                                                 outlet_depth,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        print Q, v, d
+        assert num.allclose(Q, 1.00, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 2.204, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.76, rtol=1.0e-2) #depth at outlet used to calc v
+
+    def Xtest_boyd_9(self):
+        """test_boyd_9
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=0.1  # Downward
+
+        inlet_depth=1.1504
+        outlet_depth=1.5
+
+        culvert_length=4.0
+        culvert_width=0.75
+        culvert_height=0.75
+        
+        culvert_type='pipe'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth #+0.5*v**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
+        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth, 
+                                                 outlet_depth,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        print Q, v, d
+        assert num.allclose(Q, 0.10, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 0.22, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.76, rtol=1.0e-2) #depth at outlet used to calc v
+
+
+    def Xtest_boyd_10(self):
+        """test_boyd_9
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=0.1  # Downward
+
+        inlet_depth=1.901
+        outlet_depth=1.5
+
+        culvert_length=4.0
+        culvert_width=0.75
+        culvert_height=0.75
+        
+        culvert_type='pipe'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth #+0.5*v**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
+        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth, 
+                                                 outlet_depth,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        print Q, v, d
+        assert num.allclose(Q, 1.00, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 2.204, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.76, rtol=1.0e-2) #depth at outlet used to calc v
+    
+               
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_culvert_routines, 'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
+
Index: /trunk/anuga_core/source/anuga/structures/test_culvert_routines_box_10pct.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/test_culvert_routines_box_10pct.py	(revision 7939)
+++ /trunk/anuga_core/source/anuga/structures/test_culvert_routines_box_10pct.py	(revision 7939)
@@ -0,0 +1,342 @@
+#!/usr/bin/env python
+
+
+import unittest
+import os.path
+import sys
+
+from anuga.utilities.system_tools import get_pathname_from_package
+from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
+import numpy as num
+
+
+class Test_culvert_routines_box_10pct(unittest.TestCase):
+    """
+	This unit test sets up 6 tests for various culvert conditions for a Box Culvert on a 10% Slope
+    """
+
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+
+    def test_boyd_1(self):
+        """test_boyd_1
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+
+
+        inlet_depth=0.150
+        outlet_depth=0.15
+        inlet_velocity=1.00
+        outlet_velocity=0.5
+        
+        culvert_length=10.0
+        culvert_width=3.6
+        culvert_height=1.20
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        culvert_slope=10.0  # % Downward
+        z_in = 10.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST01 Q-v-d',Q,v,d))
+        #print('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 0.5526, 1.146, 0.1339))
+        assert num.allclose(Q, 0.5526, rtol=1.0e-1) #inflow
+        assert num.allclose(v, 1.146, rtol=1.0e-1) #outflow velocity
+        assert num.allclose(d, 0.1339, rtol=1.0e-1) #depth at outlet used to calc v 
+        
+    def test_boyd_2(self):
+        """test_boyd_2
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=10  # Downward
+
+        inlet_depth=0.500
+        outlet_depth=0.700
+        inlet_velocity=1.0
+        outlet_velocity=0.50
+        
+        culvert_length=10.0
+        culvert_width=3.60
+        culvert_height=1.20
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST02 Q-v-d',Q,v,d))
+        #print ('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 2.508, 1.897, 0.367))
+        assert num.allclose(Q, 2.508, rtol=1.0e-1) #inflow
+        assert num.allclose(v, 1.897, rtol=1.0e-1) #outflow velocity
+        assert num.allclose(d, 0.367, rtol=1.0e-1) #depth at outlet used to calc v  
+
+    def test_boyd_3(self):
+        """test_boyd_3
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=10  # Downward
+
+        inlet_depth=1.800
+        outlet_depth=0.80
+        inlet_velocity=1.0
+        outlet_velocity=0.5
+        
+        culvert_length=10.0
+        culvert_width=3.60
+        culvert_height=1.20
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        #print ('%s,%.2f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.2f'%('Delta E = ',delta_total_energy))
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST03 Q-v-d',Q,v,d))
+        #print ('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 13.554, 3.329, 1.131))
+        assert num.allclose(Q, 13.554, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 3.329, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 1.131, rtol=1.0e-2) #depth at outlet used to calc v
+
+#NOTE FROM HERE DOWN THE UNITS TEST HAVE NOT BEEN AMENDED TO ALLOW VELOCITY COMPONENT TO BE USED. ONLY ABOVE 3 TESTS WORK. PM WILL FIX THE ONES BELOW WHEN THE ABOVE 3 ARE WORKING
+    def test_boyd_4(self):
+        """test_boyd_4
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=10  # Downward
+
+        inlet_depth=1.00
+        outlet_depth=0.8
+        inlet_velocity=1.0
+        outlet_velocity=0.5 
+        culvert_length=10.0
+        culvert_width=3.60
+        culvert_height=1.20
+       
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)        
+        #print ('%s,%.2f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.2f'%('Delta E = ',delta_total_energy))
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST04 Q-v-d',Q,v,d))
+        #print ('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 6.609, 2.621, 0.70))
+        assert num.allclose(Q, 6.609, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 2.621, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.70, rtol=1.0e-2) #depth at outlet used to calc v
+
+    def test_boyd_5(self):
+        """test_boyd_5
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=10  # Downward
+
+        inlet_depth=1.50
+        inlet_velocity= 1.0
+        outlet_depth=2.5
+        outlet_velocity=0.5
+        culvert_length=10.0
+        culvert_width=3.60
+        culvert_height=1.20
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)        
+        #print ('%s,%.3f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.3f'%('Delta E = ',delta_total_energy))
+        
+        #print ('%s,%.3f,%.3f,%.3f' %('ANUGAcalcsTEST05 Q-v-d',Q,v,d))
+        #print ('%s,%.3f,%.3f,%.3f' %('Spreadsheet_Boydcalcs',2.961, 0.685, 1.20))
+        assert num.allclose(Q, 2.961, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 0.685, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 1.20, rtol=1.0e-2) #depth at outlet used to calc v
+
+    def test_boyd_6(self):
+        """test_boyd_6
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=10  # Downward
+
+        inlet_depth=1.50
+        inlet_velocity= 4.0
+        outlet_depth=0.80
+        outlet_velocity=4.0
+        culvert_length=10.0
+        culvert_width=3.60
+        culvert_height=1.20
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)        
+        #print ('%s,%.3f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.3f'%('Delta E = ',delta_total_energy))
+        
+        #print ('%s,%.3f,%.3f,%.3f' %('ANUGAcalcsTEST06 Q-v-d',Q,v,d))
+        #print ('%s,%.3f,%.3f,%.3f' %('Spreadsheet_Boydcalcs',15.537, 3.597, 1.20))
+        assert num.allclose(Q, 15.537, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 3.597, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 1.20, rtol=1.0e-2) #depth at outlet used to calc v
+# =========================================================================
+# =========================================================================
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_culvert_routines_box_10pct, 'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: /trunk/anuga_core/source/anuga/structures/test_culvert_routines_box_1pct.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/test_culvert_routines_box_1pct.py	(revision 7939)
+++ /trunk/anuga_core/source/anuga/structures/test_culvert_routines_box_1pct.py	(revision 7939)
@@ -0,0 +1,342 @@
+#!/usr/bin/env python
+
+
+import unittest
+import os.path
+import sys
+
+from anuga.utilities.system_tools import get_pathname_from_package
+from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
+import numpy as num
+
+
+class Test_culvert_routines_box_1pct(unittest.TestCase):
+    """
+	This unit test sets up 6 tests for various culvert conditions for a Box Culvert on a 1% Slope
+    """
+
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+
+    def test_boyd_1(self):
+        """test_boyd_1
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+
+
+        inlet_depth=0.150
+        outlet_depth=0.15
+        inlet_velocity=1.00
+        outlet_velocity=0.5
+        
+        culvert_length=10.0
+        culvert_width=3.6
+        culvert_height=1.20
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        culvert_slope=1  # % Downward
+        z_in = 10.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST01 Q-v-d',Q,v,d))
+        #print('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 0.5526, 1.146, 0.1339))
+        assert num.allclose(Q, 0.5526, rtol=1.0e-1) #inflow
+        assert num.allclose(v, 1.146, rtol=1.0e-1) #outflow velocity
+        assert num.allclose(d, 0.1339, rtol=1.0e-1) #depth at outlet used to calc v 
+        
+    def test_boyd_2(self):
+        """test_boyd_2
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=1  # Downward
+
+        inlet_depth=0.500
+        outlet_depth=0.700
+        inlet_velocity=1.50
+        outlet_velocity=0.50
+        
+        culvert_length=10.0
+        culvert_width=3.60
+        culvert_height=1.20
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST02 Q-v-d',Q,v,d))
+        #print ('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 0.224, 0.152, 0.409))
+        assert num.allclose(Q, 0.224, rtol=1.0e-1) #inflow
+        assert num.allclose(v, 0.152, rtol=1.0e-1) #outflow velocity
+        assert num.allclose(d, 0.409, rtol=1.0e-1) #depth at outlet used to calc v  
+
+    def test_boyd_3(self):
+        """test_boyd_3
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=1  # Downward
+
+        inlet_depth=1.800
+        outlet_depth=0.80
+        inlet_velocity=1.0
+        outlet_velocity=0.5
+        
+        culvert_length=10.0
+        culvert_width=3.60
+        culvert_height=1.20
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        #print ('%s,%.2f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.2f'%('Delta E = ',delta_total_energy))
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST03 Q-v-d',Q,v,d))
+        #print ('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 13.554, 3.329, 1.131))
+        assert num.allclose(Q, 13.554, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 3.329, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 1.131, rtol=1.0e-2) #depth at outlet used to calc v
+
+#NOTE FROM HERE DOWN THE UNITS TEST HAVE NOT BEEN AMENDED TO ALLOW VELOCITY COMPONENT TO BE USED. ONLY ABOVE 3 TESTS WORK. PM WILL FIX THE ONES BELOW WHEN THE ABOVE 3 ARE WORKING
+    def test_boyd_4(self):
+        """test_boyd_4
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=1  # Downward
+
+        inlet_depth=1.00
+        outlet_depth=0.8
+        inlet_velocity=1.0
+        outlet_velocity=0.5 
+        culvert_length=10.0
+        culvert_width=3.60
+        culvert_height=1.20
+       
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)        
+        #print ('%s,%.2f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.2f'%('Delta E = ',delta_total_energy))
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST04 Q-v-d',Q,v,d))
+        #print ('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 5.164, 2.047, 0.70))
+        assert num.allclose(Q, 5.164, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 2.047, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.70, rtol=1.0e-2) #depth at outlet used to calc v
+
+    def test_boyd_5(self):
+        """test_boyd_5
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=1  # Downward
+
+        inlet_depth=1.50
+        inlet_velocity= 1.0
+        outlet_depth=1.3
+        outlet_velocity=0.5
+        culvert_length=10.0
+        culvert_width=3.60
+        culvert_height=1.20
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)        
+        #print ('%s,%.3f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.3f'%('Delta E = ',delta_total_energy))
+        
+        #print ('%s,%.3f,%.3f,%.3f' %('ANUGAcalcsTEST05Q-v-d',Q,v,d))
+        #print ('%s,%.3f,%.3f,%.3f' %('Spreadsheet_Boydcalcs',8.808, 2.039, 1.20))
+        assert num.allclose(Q, 8.808, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 2.039, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 1.20, rtol=1.0e-2) #depth at outlet used to calc v
+
+
+    def test_boyd_6(self):
+        """test_boyd_6
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=1  # Downward
+
+        inlet_depth=1.50
+        inlet_velocity= 4.0
+        outlet_depth=0.8
+        outlet_velocity=4.0
+        culvert_length=10.0
+        culvert_width=3.60
+        culvert_height=1.20
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)        
+        #print ('%s,%.3f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.3f'%('Delta E = ',delta_total_energy))
+        
+        #print ('%s,%.3f,%.3f,%.3f' %('ANUGAcalcsTEST06 Q-v-d',Q,v,d))
+        #print ('%s,%.3f,%.3f,%.3f' %('Spreadsheet_Boydcalcs',13.546, 3.136, 1.20))
+        assert num.allclose(Q, 13.546, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 3.136, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 1.20, rtol=1.0e-2) #depth at outlet used to calc v
+
+# =========================================================================
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_culvert_routines_box_1pct, 'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: /trunk/anuga_core/source/anuga/structures/test_culvert_routines_pipe_10pct.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/test_culvert_routines_pipe_10pct.py	(revision 7939)
+++ /trunk/anuga_core/source/anuga/structures/test_culvert_routines_pipe_10pct.py	(revision 7939)
@@ -0,0 +1,338 @@
+#!/usr/bin/env python
+
+
+import unittest
+import os.path
+import sys
+
+from anuga.utilities.system_tools import get_pathname_from_package
+from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
+import numpy as num
+
+
+class Test_culvert_routines_pipe_10pct(unittest.TestCase):
+    """
+	This unit test sets up 6 tests for various culvert conditions for a Circular Pipe Culvert on a 10% Slope
+    """
+
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+
+    def test_boyd_1(self):
+        """test_boyd_1
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+
+
+        inlet_depth=0.150
+        outlet_depth=0.15
+        inlet_velocity=1.00
+        outlet_velocity=0.5
+        
+        culvert_length=10.0
+        culvert_width=0.0
+        culvert_height=1.20
+        
+        culvert_type='circle'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        culvert_slope=10.0  # % Downward
+        z_in = 10.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST01 Q-v-d',Q,v,d))
+        #print('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 0.113, 0.297, 0.443))
+        assert num.allclose(Q, 0.113, rtol=1.0e-1) #inflow
+        assert num.allclose(v, 0.297, rtol=1.0e-1) #outflow velocity
+        assert num.allclose(d, 0.443, rtol=1.0e-1) #depth at outlet used to calc v 
+        
+    def test_boyd_2(self):
+        """test_boyd_2
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=10  # Downward
+
+        inlet_depth=0.500
+        outlet_depth=0.700
+        inlet_velocity=1.0
+        outlet_velocity=0.50
+        
+        culvert_length=10.0
+        culvert_width=0.0
+        culvert_height=1.20
+        culvert_width=0.0
+        culvert_type='circle'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST02 Q-v-d',Q,v,d))
+        #print ('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 0.585, 0.873, 0.687))
+        assert num.allclose(Q, 0.585, rtol=1.0e-1) #inflow
+        assert num.allclose(v, 0.873, rtol=1.0e-1) #outflow velocity
+        assert num.allclose(d, 0.687, rtol=1.0e-1) #depth at outlet used to calc v  
+
+    def test_boyd_3(self):
+        """test_boyd_3
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=10  # Downward
+
+        inlet_depth=1.800
+        outlet_depth=0.80
+        inlet_velocity=1.0
+        outlet_velocity=0.5
+        
+        culvert_length=10.0
+        
+        culvert_height=1.20
+        culvert_width=0.0
+        culvert_type='circle'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        #print ('%s,%.2f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.2f'%('Delta E = ',delta_total_energy))
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST03 Q-v-d',Q,v,d))
+        #print ('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 3.441, 3.042, 1.2))
+        assert num.allclose(Q, 3.441, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 3.042, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 1.2, rtol=1.0e-2) #depth at outlet used to calc v
+
+#NOTE FROM HERE DOWN THE UNITS TEST HAVE NOT BEEN AMENDED TO ALLOW VELOCITY COMPONENT TO BE USED. ONLY ABOVE 3 TESTS WORK. PM WILL FIX THE ONES BELOW WHEN THE ABOVE 3 ARE WORKING
+    def test_boyd_4(self):
+        """test_boyd_4
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=10  # Downward
+
+        inlet_depth=1.00
+        outlet_depth=0.8
+        inlet_velocity=1.0
+        outlet_velocity=0.5 
+        culvert_length=10.0
+        culvert_width=0.0
+        culvert_height=1.20
+       
+        culvert_type='circle'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)        
+        #print ('%s,%.2f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.2f'%('Delta E = ',delta_total_energy))
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST04 Q-v-d',Q,v,d))
+        #print ('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 1.676, 1.821, 0.91))
+        assert num.allclose(Q, 1.676, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 1.821, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.91, rtol=1.0e-2) #depth at outlet used to calc v
+
+    def test_boyd_5(self):
+        """test_boyd_5
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=10  # Downward
+
+        inlet_depth=1.50
+        inlet_velocity= 1.0
+        outlet_depth=2.5
+        outlet_velocity=0.5
+        culvert_length=10.0
+        culvert_width=0.0
+        culvert_height=1.20
+        
+        culvert_type='circle'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)        
+        #print ('%s,%.3f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.3f'%('Delta E = ',delta_total_energy))
+        
+        #print ('%s,%.3f,%.3f,%.3f' %('ANUGAcalcsTEST05 Q-v-d',Q,v,d))
+        #print ('%s,%.3f,%.3f,%.3f' %('Spreadsheet_Boydcalcs',0.759, 0.671, 1.20))
+        assert num.allclose(Q, 0.759, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 0.671, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 1.20, rtol=1.0e-2) #depth at outlet used to calc v
+
+    def test_boyd_6(self):
+        """test_boyd_6
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=10  # Downward
+
+        inlet_depth=1.50
+        inlet_velocity= 4.0
+        outlet_depth=0.80
+        outlet_velocity=4.0
+        culvert_length=10.0
+        culvert_height=1.20
+        culvert_width=0.0
+        culvert_type='circle'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)        
+        #print ('%s,%.3f,%.3f,%.3f' %('ANUGAcalcsTEST06 Q-v-d',Q,v,d))
+        #print ('%s,%.3f,%.3f,%.3f' %('Spreadsheet_Boydcalcs',3.962, 3.503, 1.20))
+        assert num.allclose(Q, 3.962, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 3.503, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 1.20, rtol=1.0e-2) #depth at outlet used to calc v
+# =========================================================================
+# =========================================================================
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_culvert_routines_pipe_10pct, 'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: /trunk/anuga_core/source/anuga/structures/test_culvert_routines_pipe_1pct.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/test_culvert_routines_pipe_1pct.py	(revision 7939)
+++ /trunk/anuga_core/source/anuga/structures/test_culvert_routines_pipe_1pct.py	(revision 7939)
@@ -0,0 +1,338 @@
+#!/usr/bin/env python
+
+
+import unittest
+import os.path
+import sys
+
+from anuga.utilities.system_tools import get_pathname_from_package
+from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
+import numpy as num
+
+
+class Test_culvert_routines_pipe_1pct(unittest.TestCase):
+    """
+	This unit test sets up 6 tests for various culvert conditions for a Circular Pipe Culvert on a 1.0% Slope
+    """
+
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+
+    def test_boyd_1(self):
+        """test_boyd_1
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+
+
+        inlet_depth=0.150
+        outlet_depth=0.15
+        inlet_velocity=1.00
+        outlet_velocity=0.5
+        
+        culvert_length=10.0
+        culvert_width=0.0
+        culvert_height=1.20
+        
+        culvert_type='circle'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        culvert_slope=1.0  # % Downward
+        z_in = 10.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST01 Q-v-d',Q,v,d))
+        #print('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 0.113, 0.297, 0.443))
+        assert num.allclose(Q, 0.113, rtol=1.0e-1) #inflow
+        assert num.allclose(v, 0.297, rtol=1.0e-1) #outflow velocity
+        assert num.allclose(d, 0.443, rtol=1.0e-1) #depth at outlet used to calc v 
+        
+    def test_boyd_2(self):
+        """test_boyd_2
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=1  # Downward
+
+        inlet_depth=0.500
+        outlet_depth=0.700
+        inlet_velocity=1.5
+        outlet_velocity=0.50
+        
+        culvert_length=10.0
+        culvert_width=0.0
+        culvert_height=1.20
+        culvert_width=0.0
+        culvert_type='circle'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        #print ('%s,%.3f,%.3f,%.3f' %('ANUGAcalcsTEST02 Q-v-d',Q,v,d))
+        #print ('%s,%.3f,%.3f,%.3f' %('Spreadsheet_Boydcalcs', 0.108, 0.152, 0.721))
+        assert num.allclose(Q, 0.108, rtol=1.0e-1) #inflow
+        assert num.allclose(v, 0.152, rtol=1.0e-1) #outflow velocity
+        assert num.allclose(d, 0.721, rtol=1.0e-1) #depth at outlet used to calc v  
+
+    def test_boyd_3(self):
+        """test_boyd_3
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=1  # Downward
+
+        inlet_depth=1.800
+        outlet_depth=0.80
+        inlet_velocity=1.0
+        outlet_velocity=0.5
+        
+        culvert_length=10.0
+        
+        culvert_height=1.20
+        culvert_width=0.0
+        culvert_type='circle'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        #print ('%s,%.2f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.2f'%('Delta E = ',delta_total_energy))
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST03 Q-v-d',Q,v,d))
+        #print ('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 3.441, 3.042, 1.2))
+        assert num.allclose(Q, 3.441, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 3.042, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 1.2, rtol=1.0e-2) #depth at outlet used to calc v
+
+#NOTE FROM HERE DOWN THE UNITS TEST HAVE NOT BEEN AMENDED TO ALLOW VELOCITY COMPONENT TO BE USED. ONLY ABOVE 3 TESTS WORK. PM WILL FIX THE ONES BELOW WHEN THE ABOVE 3 ARE WORKING
+    def test_boyd_4(self):
+        """test_boyd_4
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=1  # Downward
+
+        inlet_depth=1.00
+        outlet_depth=0.9
+        inlet_velocity=1.0
+        outlet_velocity=0.5 
+        culvert_length=10.0
+        culvert_width=0.0
+        culvert_height=1.20
+       
+        culvert_type='circle'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)        
+        #print ('%s,%.2f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.2f'%('Delta E = ',delta_total_energy))
+        #print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST04 Q-v-d',Q,v,d))
+        #print ('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', 1.559, 1.694, 0.91))
+        assert num.allclose(Q, 1.559, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 1.694, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 0.91, rtol=1.0e-2) #depth at outlet used to calc v
+
+    def test_boyd_5(self):
+        """test_boyd_5
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=1  # Downward
+
+        inlet_depth=1.50
+        inlet_velocity= 1.0
+        outlet_depth=1.3
+        outlet_velocity=0.5
+        culvert_length=10.0
+        culvert_width=0.0
+        culvert_height=1.20
+        
+        culvert_type='circle'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)        
+        #print ('%s,%.3f'%('SPEC_E = ',inlet_specific_energy))
+        #print ('%s,%.3f'%('Delta E = ',delta_total_energy))
+        
+        #print ('%s,%.3f,%.3f,%.3f' %('ANUGAcalcsTEST05 Q-v-d',Q,v,d))
+        #print ('%s,%.3f,%.3f,%.3f' %('Spreadsheet_Boydcalcs',2.258, 1.996, 1.20))
+        assert num.allclose(Q, 2.258, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 1.996, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 1.20, rtol=1.0e-2) #depth at outlet used to calc v
+
+    def test_boyd_6(self):
+        """test_boyd_6
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=1  # Downward
+
+        inlet_depth=1.50
+        inlet_velocity= 4.0
+        outlet_depth=0.80
+        outlet_velocity=4.0
+        culvert_length=10.0
+        culvert_height=1.20
+        culvert_width=0.0
+        culvert_type='circle'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+                                                 outlet_depth,
+                                                 inlet_velocity,
+                                                 outlet_velocity,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)        
+        #print ('%s,%.3f,%.3f,%.3f' %('ANUGAcalcsTEST06 Q-v-d',Q,v,d))
+        #print ('%s,%.3f,%.3f,%.3f' %('Spreadsheet_Boydcalcs',3.472, 3.070, 1.20))
+        assert num.allclose(Q, 3.472, rtol=1.0e-2) #inflow
+        assert num.allclose(v, 3.070, rtol=1.0e-2) #outflow velocity
+        assert num.allclose(d, 1.20, rtol=1.0e-2) #depth at outlet used to calc v
+# =========================================================================
+# =========================================================================
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_culvert_routines_pipe_1pct, 'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
