Index: /anuga_work/development/culvert_flow/culvert_boyd_channel.py
===================================================================
--- /anuga_work/development/culvert_flow/culvert_boyd_channel.py	(revision 5427)
+++ /anuga_work/development/culvert_flow/culvert_boyd_channel.py	(revision 5427)
@@ -0,0 +1,740 @@
+""" 
+
+Ole Check Culvert Routine from Line 258
+
+Although it is Setup as a Culvert with the Opening presented vertically, 
+for now the calculation of flow rate is assuming a horizontal hole in the
+ground (Fix this Later)
+
+MOST importantly 2 things...
+1. How to use the Create Polygon Routine to enquire Depth ( or later energy)
+infront of the Culvert
+
+Done (Ole)
+
+2. How to apply the Culvert velocity and thereby Momentum to the Outlet
+Ject presented at the Outlet
+
+
+
+Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction
+
+This Version CALCULATES the Culvert Velocity and Uses it to establish 
+The Culvert Outlet Momentum
+
+The Aim is to define a Flow Transfer function that Simulates a Culvert
+by using the Total Available Energy to Drive the Culvert
+as Derived by determining the Difference in Total Energy
+between 2 Points, Just Up stream and Just Down Stream of the Culvert
+away from the influence of the Flow Abstraction etc..
+
+This example includes a Model Topography that shows a
+TYPICAL Headwall Configuration
+
+The aim is to change the Culvert Routine to Model more precisely the
+abstraction
+from a vertical face.
+
+The inflow must include the impact of Approach velocity.
+Similarly the Outflow has MOMENTUM Not just Up welling as in the
+Horizontal Style
+abstraction
+
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.shallow_water import Domain
+from anuga.shallow_water.shallow_water_domain import Reflective_boundary
+from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
+from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
+from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
+from anuga.utilities.polygon import plot_polygons
+
+from math import pi,sqrt
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+
+
+# Open file for storing some specific results...
+fid = open('Culvert_Headwall_VarM.txt', 'w')
+
+length = 40.
+width = 5.
+
+#dx = dy = 1           # Resolution: Length of subdivisions on both axes
+#dx = dy = .5           # Resolution: Length of subdivisions on both axes
+dx = dy = .25           # Resolution: Length of subdivisions on both axes
+#dx = dy = .1           # Resolution: Length of subdivisions on both axes
+
+# OLE.... How do I refine the resolution... in the area where I have the Culvert Opening ???......
+#  Can I refine in a X & Y Range ???
+points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
+                                               len1=length, len2=width)
+domain = Domain(points, vertices, boundary)   
+domain.set_name('culv_dev_HW_Var_Mom')                 # Output name
+domain.set_default_order(2)
+domain.H0 = 0.01
+domain.tight_slope_limiters = True
+domain.set_minimum_storable_height(0.02)
+
+s='Size'+str(len(domain))
+fid.write(s)
+print s
+
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+
+# Define the topography (land scape)
+def topography(x, y):
+    """Set up a weir
+    
+    A culvert will connect either side of an Embankment with a Headwall type structure
+    The aim is for the Model to REALISTICALY model flow through the Culvert
+    """
+    # General Slope of Topography
+    z=-x/100
+    floorhgt = 5
+    embank_hgt=10
+    embank_upslope=embank_hgt/5
+    embank_dnslope=embank_hgt/2.5
+    # Add bits and Pieces to topography
+    N = len(x)
+    for i in range(N):
+
+        # Sloping Embankment Across Channel
+        
+        if 0.0 < x[i] < 7.51:
+            z[i]=z[i]+5.0
+        if 7.5 < x[i] < 10.1:
+            if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: # Cut Out Segment for Culvert FACE
+               z[i]=z[i]+5.0
+            else:
+               z[i] +=  embank_upslope*(x[i] -5.0)    # Sloping Segment  U/S Face
+        if 10.0 < x[i] < 12.1:
+            if  2.0 < y[i]  < 3.0: # Cut Out Segment for Culvert (open Channel)
+                #z[i] +=  z[i]+5-(x[i]-10)*2              # Sloping Channel in Embankment
+                z[i] +=  embank_hgt                 # Flat Crest of Embankment
+            else:
+                z[i] +=  embank_hgt                 # Flat Crest of Embankment
+        if 12.0 < x[i] < 14.5:
+            if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5: # Cut Out Segment for Culvert FACE
+               z[i]=z[i]
+            else:
+               z[i] += embank_hgt-embank_dnslope*(x[i] -12.0)       # Sloping D/S Face
+      
+        
+        # Constriction
+        #if 27 < x[i] < 29 and y[i] > 3:
+        #    z[i] += 2        
+        
+        # Pole
+        #if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2:
+        #    z[i] += 2
+
+        # HOLE For Pit at Opening[0]
+        #if (x[i]-4)**2 + (y[i]-1.5)**2<0.75**2:
+        #  z[i]-=1
+
+        # HOLE For Pit at Opening[1]
+        #if (x[i]-20)**2 + (y[i]-3.5)**2<0.5**2:
+        #  z[i]-=1
+
+    return z
+
+
+domain.set_quantity('elevation', topography)  # Use function for elevation
+domain.set_quantity('friction', 0.01)         # Constant friction 
+domain.set_quantity('stage',
+                    expression='elevation')   # Dry initial condition
+
+
+
+
+#------------------------------------------------------------------------------
+# Setup specialised forcing terms
+#------------------------------------------------------------------------------
+
+ #   NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet)
+ # 
+ # The First Attempt has a Simple Horizontal Circle as a Hole on the Bed
+ # Flow Is Removed at a Rate of INFLOW 
+ # Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges
+ #
+ # This SHould be changed to a Vertical Opening Both BOX and Circular
+ # There will be several Culvert Routines such as:
+ # CULVERT_Simple_FLOOR
+ # CULVERT_Simple_WALL
+ # CULVERT_Eqn_Floor
+ # CULVERT_Eqn_Wall
+ # CULVERT_Tab_Floor
+ # CULVERT_Tab_Wall
+ # BRIDGES.....
+ # NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!!
+ 
+ #  COULD USE EPA SWMM Model !!!!
+ 
+ 
+ 
+class Culvert_flow:
+    """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...
+	
+    """	
+
+    def __init__(self,
+                 domain,
+                 label=None, 
+                 description=None,
+                 end_point0=None, 
+                 end_point1=None,
+                 width=None,
+                 height=None,
+                 manning=None,          # Mannings Roughness for Culvert
+                 inv_point0=None,       # Invert level if not the same as the Elevation on the Domain
+                 inv_point1=None,       # Invert level if not the same as the Elevation on the Domain
+                 loss_exit=None,
+                 loss_entry=None,
+                 loss_bend=None,
+                 loss_special=None,
+                 blockage_topdwn=None,
+                 blockage_bottup=None,
+                 verbose=False):
+        
+        from Numeric import sqrt, sum
+
+        if height is None: height = width
+        if manning is None: manning = 0.012   # Set a Default Mannings 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
+        #sum_loss=loss_exit+loss_entry+loss_bend+loss_special    # Note if losses are a function can't do this !!!!
+
+        # Create the fundamental culvert polygons from POLYGON 
+        P = create_culvert_polygons(end_point0,
+                                    end_point1,
+                                    width=width,   
+                                    height=height)
+        
+        if verbose is True:
+            pass
+            #plot_polygons([[end_point0, end_point1],
+            #               P['exchange_polygon0'],
+            #               P['exchange_polygon1'],
+            #               P['enquiry_polygon0'],
+            #               P['enquiry_polygon1']],
+            #              figname='culvert_polygon_output')
+        
+        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.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
+        self.vector = P['vector']
+        self.distance = P['length']
+        self.verbose = verbose
+        self.width = width
+        self.height = height
+        self.last_time = 0.0        
+        self.temp_keep_delta_t=0.0
+        # Create objects to update momentum (a bit crude at this stage)
+        self.xmom_forcing0 = General_forcing(domain,
+                                             'xmomentum',
+                                             polygon=P['exchange_polygon0'])
+
+        self.xmom_forcing1 = General_forcing(domain,
+                                             'xmomentum',
+                                             polygon=P['exchange_polygon1'])
+
+        self.ymom_forcing0 = General_forcing(domain,
+                                             'ymomentum',
+                                             polygon=P['exchange_polygon0'])
+
+        self.ymom_forcing1 = General_forcing(domain,
+                                             'ymomentum',
+                                             polygon=P['exchange_polygon1'])                
+
+        # Print something
+        
+        s= 'Culvert Effective Length ='+str( self.distance)
+        print s
+        fid.write(s)
+   
+        s= 'Culvert Direction is '+str( self.vector)
+        print s
+        fid.write(s)
+        s= 'Point 1m Up Stream is X,Y ='
+        print s
+        fid.write(s)
+        s= 'Point 1m Down Stream is X,Y ='
+        print s
+        fid.write(s)
+       
+
+        
+    def __call__(self, domain):
+        from anuga.utilities.numerical_tools import mean
+        from anuga.utilities.polygon import inside_polygon
+        from anuga.config import g, epsilon
+        from Numeric import take, sqrt
+        
+        #  OLE    I think this routine gets Called twice every time Step, so that each time the second time around Delta_t =0 and therefore Momentum is Set to ZERO
+		# Why does this NEED to be calle d twice ???
+
+
+        # Get average water depths at each opening
+
+        time = domain.get_time()
+        openings = self.openings   # There are two Opening  [0] and [1]   so this happens TWICE Right ???????????
+        for i, opening in enumerate(openings):
+            stage = domain.quantities['stage'].get_values(location='centroids',
+                                                          indices=opening.indices)
+            elevation = domain.quantities['elevation'].get_values(location='centroids',
+                                                                  indices=opening.indices)
+
+            # Indices corresponding to energy enquiry field for this opening
+            coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y)
+            idx = inside_polygon(coordinates, self.enquiry_polygons[i]) 
+            
+            if self.verbose:
+                pass
+                #print 'Opening %d: Got %d points in enquiry polygon:\n%s'\
+                #      %(i, len(idx), self.enquiry_polygons[i])
+                
+            
+            # Get average model values for points in ?????? WHAT IF THE DOMAIN IS DRY ???????
+            # enquiry polygon for this opening
+            dq = domain.quantities
+            stage = dq['stage'].get_values(location='centroids', indices=idx)
+            xmomentum = dq['xmomentum'].get_values(location='centroids', indices=idx)
+            ymomentum = dq['ymomentum'].get_values(location='centroids', indices=idx)                        
+            elevation = dq['elevation'].get_values(location='centroids', indices=idx)
+            depth = stage - elevation            
+
+            # Compute mean velocity in the area in front of the culvert (taking zero depths into account)
+            ux = xmomentum/(depth+epsilon)   # Does epsilon.... handle a Dry Cell Case ????????
+            uy = ymomentum/(depth+epsilon)
+            v = mean(sqrt(ux**2+uy**2))
+            d = mean(depth)
+            w = mean(stage)
+            z = mean(elevation)
+            
+            # Ratio of depth to culvert height
+            #ratio = d/(2*self.height)  # Is Height Now REAL Height or Radius   if Radius then Need 2 else NOT !!
+            ratio = d/self.height  # Use this if Height is NOT Radius!
+            if ratio > 1.0:    # Assume culvert is running full & 
+                ratio = 1.0    # under pressure. Note this is usually ~ 1.35
+            opening.ratio = ratio
+                
+
+            # Average measure of total energy (D + V^2/2g) in enquiry field in front of this opening
+            Es = d + 0.5*v**2/g  #  SPECIFIC ENERGY
+            Et = w + 0.5*v**2/g  #  TOTAL ENERGY
+            opening.energy = Et
+            opening.vel=v
+            
+
+            # Store current average stage and depth at enquiry field with each opening object 
+            opening.depth = d
+            opening.stage = w
+            opening.elevation = z
+
+        #  End of the FOR  LOOP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+        # we now need to deal with each opening
+
+        # Now we are done calculating energy etc for each opening.
+        # At this point we can work on the transfer functions etc.
+        
+        # Handy values (all calculated at enquiry polygon -
+        # if you need them at exchange polygons we can easily do that.)
+
+
+        #Then we need delta_h or better delta_E to be used to drive water through culvert      IS THIS FROM EXCHANGE or ENQUIRY POLYGON ?????? Need to Know !!!!
+        #delta_h = openings[0].energy - openings[1].energy
+        delta_h = openings[0].stage - openings[1].stage
+        delta_Et = openings[0].energy - openings[1].energy
+        delta_z = openings[1].elevation - openings[0].elevation
+        culvert_slope=(delta_z/self.distance)      # COULD USE CULVERT SLOPE To Determine Internal Flow Velocity That is Not Inlet Or Outlet Velocity
+        s= 'Delta_Energy = %.3f'%(delta_Et)
+        print s
+        fid.write(s)
+        
+        # Ideas.....  GET RID OF THIS  !!!!!!!!!!!!!!!!
+        if openings[0].depth > 0 and openings[1].depth > 0:
+            flow_rate = delta_h
+        else:
+            flow_rate = 0.0
+
+        if openings[0].depth > self.height:
+            # This could be usefull mayby
+            #print 'Inlet has been overflowed'
+            pass
+
+        if openings[1].depth > self.height:
+            # This could be usefull mayby
+            #print 'Outlet has been overflowed'            
+            pass
+        #    ====================== NOW ENTER INTO THE CULVERT EQUATIONS AS DEFINED BY BOYD GENERALISED METHOD
+        #  ==    The quantity of flow passing through a culvert is controlled by many factors
+        # ==     It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation
+        # ==    Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation
+        # ==    Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled 
+        
+        if delta_Et > 0:
+            #print 'Flow U/S ---> D/S'
+            inlet=openings[0]
+            outlet=openings[1]
+        else:
+            #print 'Flow D/S ---> U/S'
+            inlet=openings[1]
+            outlet=openings[0]
+
+        if inlet.depth < 0.01:   #epsilon:      10mm Minium should be OK....for NOW
+            # Do nothing because No water over Opening	That is Set Flow to Zero!!
+            #print 'NO water over Inlet yet!'
+            self.openings[0].rate = 0
+            self.openings[1].rate = 0
+            Qb1=0.0
+            Qb2=0.0
+            Qb3tw=0.0
+            Qb3dc=0.0
+
+
+        else: #For PIPES
+            print 'Water over Inlet NOW!'
+            if self.width==self.height:  # Need something better to Flag Diamater !!!!!!! Can't we just have DIAMETER as well as Width & Height ???
+                pass
+                #Q1[t]= 0.421*g**0.5*Diam**0.87*HW**1.63    # Inlet Ctrl Inlet Unsubmerged
+                #Q2[t]= 0.530*g**0.5*Diam**1.87*HW**0.63    # Inlet Ctrl Inlet Submerged
+                #PipeDcrit=
+                #Delta_E=HW-TW
+            else:     #FOR BOX CULVERT
+			# Why does it not KNOW these Losses from __INIT__  ????
+                #sum_loss=loss_exit+loss_entry+loss_bend+loss_special    # Note if losses are a function can't do this !!!!
+                sum_loss=1.5
+                manning=0.012     # THis should come from Culvert data !!!
+                print 'Depth over inlet =',inlet.depth
+                Qb1=0.540*g**0.5*self.width*inlet.depth**1.50                    # Flow based on Inlet Ctrl Inlet Unsubmerged
+                Qb2=0.702*g**0.5*self.width*self.height**0.89*inlet.depth**0.61  # Flow based on Inlet Ctrl Inlet Submerged
+                print'Qb1= %.6f, Qb2 = %.6f'%(Qb1,Qb2)
+                Qb_Inlet=min(Qb1,Qb2)
+                # NOW FOR OUTLET CONTROL 
+                # Need to determine the depth at the outlet relative to the depth of flow in the Culvert
+#   THIS IS WHERE THE TROUBLE STARTS
+                # USUALLY YOU NEED TO ITERATE THECALCULATION OF Q as itis Reliant on Knowing Dcrit which is reliant on Q
+                # Instead we will calculate Q assuming TW and then Q again using Dcrit assuming the Qtw
+                #If TW >H then Assume Dcrit =H else Assume Dcrit= TW
+                if outlet.depth > self.height:  # The Outlet is Submerged
+                    outlet_culv_depth=self.height
+                    flow_area=self.width*self.height
+                    perimeter=2.0*(self.width+self.height)
+                    hyd_rad=flow_area/perimeter
+                elif outlet.depth==0.0: 
+                    outlet_culv_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
+                    flow_area=self.width*inlet.depth
+                    perimeter=(self.width+2.0*inlet.depth)
+                    hyd_rad=flow_area/perimeter
+                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
+                    outlet_culv_depth=outlet.depth
+                    flow_area=self.width*outlet.depth
+                    perimeter=(self.width+2.0*outlet.depth)
+                    hyd_rad=flow_area/perimeter
+                print 'hydrad at outlet= ',hyd_rad
+                Vel3=sqrt(delta_Et/((sum_loss/2*g)+(manning**2*self.distance)/hyd_rad**1.33333))  # Outlet Ctrl usingTW as Control Depth Either TW or H
+                Qb3tw=flow_area*Vel3
+                Qb3dc=Qhold=Qb_Inlet
+                Qcheck=10.0
+                while Qcheck > 0.001:
+                    Dcrit = 0.4672*(Qb3dc/self.width)**0.66666   # This assumes that Flow is at Critical Depth  Not Alsways True...
+					# Better to Use Mannings Equation to Estimate Q and iterate Depth in the Culvert to Match the Energy Loss
+                    print 'Dcrit =',Dcrit
+                    if Dcrit > self.height:  # 
+                        flow_area=self.width*self.height
+                        perimeter=2.0*(self.width+self.height)
+                        hyd_rad=flow_area/perimeter
+                    else:
+                        flow_area=self.width*Dcrit
+                        perimeter=(self.width+2.0*Dcrit)
+                        hyd_rad=flow_area/perimeter
+                    print 'Adjusted hydrad at outlet= ',hyd_rad
+                
+                    Qb3dc=flow_area*sqrt(delta_Et/((sum_loss/2*g)+(manning**2*self.distance)/hyd_rad**1.33333))  # Outlet Ctrl using Dcrit as Control Depth Either Dcrit or H
+                    Qcheck =abs(Qhold-Qb3dc)
+                    Qhold=Qb3dc
+                    print 'Qcheck Qb3dc  ',Qcheck,Qb3dc
+                flow_rate_control=min(Qb1,Qb2) #
+                print 'OutletCtrl Qbtw and Qbdc= ',Qb3tw,Qb3dc
+                print 'Flow Rate Control = ',flow_rate_control
+				
+				# NOW UPDATE Flows on the DOMAIN once again have to check direction of FLOW 
+                if delta_Et > 0:
+                    print 'Flow U/S ---> D/S'
+                    self.openings[0].rate=-flow_rate_control
+                    self.openings[1].rate=flow_rate_control
+                else:
+                    print 'Flow D/S ---> U/S'
+                    self.openings[1].rate=-flow_rate_control
+                    self.openings[0].rate=flow_rate_control     # Flow Rates have been set  NOW  Add jet in the form of absolute momentum to outlet
+					
+					
+                culv_froude=sqrt(flow_rate_control**2*self.width/(g*flow_area**3))
+                print 'Froude in Culvert = ',culv_froude
+
+                #+++++++++++ NOW DETERMINE MOMENTUM AT THE OUTLET !!!! +++++++++++++++
+
+
+                # This Should be Based on the VELOCITY in the Culvert 
+                # Based on the Flow Depth in the Culvert for part full  So need to determine actual flow depth in the culvert based on slope !!!
+                # or Flow Jet of the Pressurised Culvert if FUll
+                # If Part Full Flow Calculate Part Full Depth
+                # Based on Depth Calculate Area... the Vel = flow_rate_control / Area
+                # Need to Break Velocity into X & Y Components
+                #self.openings[1].set_quantity_values(delta_h*speed, 'xmomentum')
+                # Previous Calculated Depth/ Culvert Heigh Ratio Use it to Determine Velocity !!!! FIX LAter
+                
+                #if (ratio*self.area)== 0: # Don't Really need this already established water depth here
+                #   outlet_vel=0.0
+                #else:
+
+
+                # FIXME (Ole): Shouldn't this be openings[1]??
+                # This attempts to resolve the Velocity in the Culvert at the Entrance and Outlet.... Again if flow reversal Occurs  then need to swap END !!!
+                #inlet_vel =(flow_rate_control/(openings[0].ratio*openings[0].area))
+                #inlet_dep = openings[0].depth*openings[0].ratio #?????
+                #inlet_mom = outlet_vel*outlet_dep
+                #inlet_E_Loss= 0.5*0.5*inlet_vel**2/g
+                #  NEED to Calculate Critical depth inthe Culvert for Box & Pipe Sections.. If Greater than Depth on the Down Stream Flood Plain use Dcrit
+                outlet_vel =(flow_rate_control/flow_area)
+                outlet_dep = flow_area/self.width #?????
+                outlet_mom = outlet_vel*outlet_dep
+                outlet_E_Loss= 0.5*0.5*outlet_vel**2/g   # Ke v^2/2g
+                print 'Outlet Velocity =',outlet_vel
+
+                #barrel_vel=(inlet_vel+outlet_vel)/2
+                # use this to determine Barrel Loss
+                
+                # Eventually will Need Momentum in X & Y Components based on the orientation of
+                # the culvert from X0,Y0, X1, Y1 from Create Polygon Routine
+                # YES - use self.vector which is a unit vector in the direction of the culvert.
+				# Will also need to Check Normal Depth For a Steep Culvert using the Mannings Equation to determine the Flow Depth in the Culvert this requires Iteration and a Simpsons Rule Approach to Refine Depth
+				#    So that Qcalculated = Qestimated....
+                    
+                outlet_mom_x, outlet_mom_y = self.vector * outlet_mom
+                #print 'Directional momentum', outlet_mom_x, outlet_mom_y
+                xmomentum_rate = outlet_mom_x
+                ymomentum_rate = outlet_mom_y
+                print 'X Y MOM ',xmomentum_rate,ymomentum_rate
+                # Update the momentum forcing terms with directional momentum at the outlet
+				# >>>>>>>>>>>>>>OLE Passing through this the 2nd time ALWAYS sets delta_t to 0 so that Momentum is also set to 0.0    This is NOT RIGHT !!!!  THIS SEGMENT OF CODE NEEDS ONLY TO BE RUn FOR THE OUTLET !!!
+				# ONLY NEEDS TO BE CALLED ONCE for the CULVERT not TWICE for Each Opening as a check of each opening is DONE HERe any WAY !!!
+				
+				# I have tried to trick it by introducing this term temp_keep_Delta_t...... but this is really bad !!!      How do we avoid this Situation ????
+                delta_t = time - self.last_time
+                if delta_t==0: 
+                    delta_t=self.temp_keep_delta_t
+                    xmomentum_rate=self.xmom_forcing1.rate
+                    ymomentum_rate=self.ymom_forcing1.rate
+                self.temp_keep_delta_t=delta_t
+                print 'time and delta t = ',time,delta_t
+                if delta_t > 0.0:
+                    print 'delta_t Calc MOM'
+                    xmomentum_rate = outlet_mom_x - self.xmom_forcing1.value
+                    xmomentum_rate /= delta_t
+                    
+                    ymomentum_rate = outlet_mom_y - self.ymom_forcing1.value
+                    ymomentum_rate /= delta_t
+                    print 'X Y MOM_RATE ',xmomentum_rate,ymomentum_rate
+                    #ymomentum_rate = 3
+                #else:
+                #    xmomentum_rate = ymomentum_rate = 0.0
+
+                # Set momentum rates for outlet jet
+                self.xmom_forcing1.rate = xmomentum_rate
+                self.ymom_forcing1.rate = ymomentum_rate
+
+                # Remember this value for next step (IMPORTANT)    
+                self.xmom_forcing1.value = outlet_mom_x
+                self.ymom_forcing1.value = outlet_mom_y                    
+
+
+                if int(domain.time*100) % 100 == 0:
+                    s = 'T=%.5f, Culvert Discharge = %.3f  Culv. Vel. %0.3f'\
+                        %(time, flow_rate_control, outlet_vel)
+                    s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
+                         %(outlet_dep, outlet_mom_x,outlet_mom_y)
+                    s += ' Momentum rate: (%.4f, %.4f)'\
+                         %(xmomentum_rate, ymomentum_rate)                    
+                    s+='Outlet Vel= %.3f Outlet Dep = %.3f'\
+                         %(outlet_vel,outlet_dep)
+                    fid.write(s)
+                    print s
+         
+       
+        # 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.xmom_forcing0(domain)
+        self.ymom_forcing0(domain)                
+        self.xmom_forcing1(domain)
+        self.ymom_forcing1(domain)        
+
+            
+
+        # Print out flows every 1 seconds
+        if int(time*100) % 100 == 0:
+            s = 'Time %.6f\n' %time
+            s += '    Opening[0]: d=%.3f, v=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\
+                 %(openings[0].depth,
+                   openings[0].vel,
+                   openings[0].area,
+                   openings[0].energy,
+                   openings[0].ratio)                   
+            s += '    Opening[1]: d=%.3f, v=%.3f,  A=%.3f, E=%.3f, r=%.3f\n'\
+                 %(openings[1].depth,
+                   openings[1].vel,
+                   openings[1].area,
+                   openings[1].energy,
+                   openings[1].ratio)                  
+            s += '    Distance=%.2f, Qb1=%.3f, Qb2=%.3f, Qb3tw=%.3f,Qb3dc=%.3f\n'\
+                 %(self.distance,
+                   Qb1,
+                   Qb2,
+                   Qb3tw,Qb3dc)
+            
+            print s
+            fid.write(s)
+
+        # Store value of time    
+        self.last_time = time
+        #self.temp_keep_delta_t=delta_t
+#------------------------------------------------------------------------------
+# Setup culvert inlets and outlets in current topography
+#------------------------------------------------------------------------------
+
+# Define culvert inlet and outlets
+# NEED TO ADD Mannings N as Fixed Value or Function
+# Energy Loss Coefficients as Fixed or Function
+# Also Set the Shape & Gap Factors fo rthe Enquiry PolyGons
+# ALSO Allow the Invert Level to be provided by the USER
+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,                       
+                       verbose=True)
+                       
+domain.forcing_terms.append(culvert)
+
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+#Bi = Dirichlet_boundary([0.5, 0.0, 0.0])          # Inflow based on Flow Depth (0.5m) and Approaching Momentum !!!
+Bi = Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
+Br = Reflective_boundary(domain)              # Solid reflective wall
+Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Setup Application of  specialised forcing terms
+#------------------------------------------------------------------------------
+
+# This is the new element implemented by Ole to allow direct input of Inflow in m^3/s
+fixed_flow = Inflow(domain,
+                    rate=20.00,
+                    center=(2.1, 2.1),
+                    radius=1.261566)   #   Fixed Flow Value Over Area of 5m2 at 1m/s = 5m^3/s
+
+#	     flow=file_function('Q/QPMF_Rot_Sub13.tms'))        # Read Time Series in  from File
+#             flow=lambda t: min(0.01*t, 0.01942))                             # Time Varying Function   Tap turning up  	
+
+domain.forcing_terms.append(fixed_flow)
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+temp_keep_delta_t=0.0
+
+for t in domain.evolve(yieldstep = 0.1, finaltime = 20):
+    pass
+    #if int(domain.time*100) % 100 == 0:
+    #    domain.write_time()
+    
+    #if domain.get_time() >= 4 and tap.flow != 0.0:
+    #    print 'Turning tap off'
+    #    tap.flow = 0.0
+	
+    #if domain.get_time() >= 3 and sink.flow < 0.0:
+    #    print 'Turning drain on'
+    #    sink.flow = -0.8	
+    # Close
+
+fid.close()
+
+
+#------------------------------------------------------------------------------
+# Query output
+#------------------------------------------------------------------------------
+
+from anuga.shallow_water.data_manager import get_flow_through_cross_section
+
+swwfilename = domain.get_name()+'.sww'  # Output name from script
+print swwfilename
+
+polyline = [[17., 0.], [17., 5.]]
+
+time, Q = get_flow_through_cross_section(swwfilename, polyline, verbose=True)
+
+from pylab import ion, plot
+ion()
+plot(time, Q)
+raw_input('done')
