Index: /anuga_work/development/Hinwood_2008/run_dam.py
===================================================================
--- /anuga_work/development/Hinwood_2008/run_dam.py	(revision 5502)
+++ /anuga_work/development/Hinwood_2008/run_dam.py	(revision 5503)
@@ -61,9 +61,10 @@
          friction=0.012,  # planed wood. http://www.lmnoeng.com/manningn.htm
          outputdir_name=None,
-         run_type=0):
+         run_type=0,
+         end_tag = '_limiterD'):
 
     
     basename = 'zz_' + metadata_dic['scenario_id']
-    end_tag = '_limiterD'
+    
     if run_type == 1:
         outputdir_name += '_test' + end_tag
@@ -261,5 +262,5 @@
     # 4 is 0.01
     # 5 is 0.001
-    run_type = 4
+    run_type = 5
     #for run_data in [scenarios[5]]:
     #scenarios = scenarios[2:]
@@ -269,4 +270,5 @@
                              run_data,
                              run_type = run_type,
-                             outputdir_name=run_data['scenario_id'])
+                             outputdir_name=run_data['scenario_id'],
+                             end_tag='_limiterD')
         #gauges_for_slope(pro_instance.outputdir,[run_data])
Index: /anuga_work/development/Hinwood_2008/slope.py
===================================================================
--- /anuga_work/development/Hinwood_2008/slope.py	(revision 5502)
+++ /anuga_work/development/Hinwood_2008/slope.py	(revision 5503)
@@ -10,5 +10,5 @@
 
 from Numeric import arange, array, zeros, Float, where, greater, less, \
-     compress, argmin, choose
+     compress, argmin, choose, searchsorted
 
 from anuga.fit_interpolate.interpolate import interpolate_sww2csv
@@ -69,6 +69,6 @@
     
     for i in arange(n_slope_locations):
-        slope_locations[i] = (locations[i+1] + locations[i+1])/2.
         delta_locations[i] = (locations[i+1] - locations[i])
+        slope_locations[i] = locations[i] + 0.5*delta_locations[i]
     
     for j in arange(n_time):
@@ -118,5 +118,5 @@
 
     # Add the contour line levels to the colorbar
-    cbar.ax.set_ylabel('verbosity coefficient')
+    cbar.ax.set_ylabel('stage slope')
     #cbar.add_lines(CS2)
          
@@ -146,7 +146,11 @@
         
     # Can't seem to reshape this info once it is in the function
-    CS = contourf(x_data, times, z_data, [-1,0.6,0.8,1,2,4],
-                  colors = ('black', 'r', 'g', 'b','r'),
-                  #cmap=cm.bone,
+    #CS = contourf(x_data, times, z_data, [-1,0.6,0.8,1,2,4],
+     #             colors = ('black', 'r', 'g', 'b','r'),
+      #            #cmap=cm.bone,
+       #           origin=origin)
+    CS = contourf(x_data, times, z_data, 10,
+                  #colors = ('black', 'r', 'g', 'b','r'),
+                  cmap=cm.bone,
                   origin=origin)
     
@@ -168,5 +172,5 @@
 
     # Add the contour line levels to the colorbar
-    cbar.ax.set_ylabel('verbosity coefficient')
+    cbar.ax.set_ylabel('Froude Number')
     #cbar.add_lines(CS2)
          
@@ -190,14 +194,16 @@
             anuga_break_times.append( \
                 break_time -  run_data['ANUGA_start_time'])
-        plot_title = "Stage slope " + id
         stage_file = pro_instance.outputdir + "fslope_stage_" + end
+        plot_title = "Stage slope " + id + "\n file:" + \
+                     outputdir_name + "_slope_stage" + plot_type
+        print "Creating ", stage_file
         save_as = pro_instance.plots_dir + sep + \
                             outputdir_name + "_slope_stage" + plot_type
         times, locations, slopes = load_slopes(stage_file)
-        times, slopes = get_band(anuga_break_times[0]-TIME_BORDER,
-                                   100, times, slopes, 0)
-        locations, slopes = get_band(
-            run_data['break_xs'][0]- 2*LOCATION_BORDER,
-            100, locations, slopes, -1)
+        #times, slopes = get_band(anuga_break_times[0]-TIME_BORDER,
+         #                          100, times, slopes, 0)
+        #locations, slopes = get_band(
+         #   min(run_data['break_xs'])- 2*LOCATION_BORDER,
+          #  100, locations, slopes, -1)
         graph_contours(times, locations, slopes,
                        plot_title=plot_title,
@@ -221,7 +227,8 @@
             anuga_break_times.append( \
                 break_time -  run_data['ANUGA_start_time'])
-        plot_title = "Froude Number" + id
-            
+        plot_title = "Froude Number" + id + "\n file:" + \
+                     outputdir_name + "_froude" + plot_type
         froude_file = pro_instance.outputdir + "fslope_froude_" + end
+        print "Creating ", froude_file
         save_as = pro_instance.plots_dir + sep + \
                             outputdir_name + "_froude" + plot_type
@@ -229,6 +236,7 @@
         dtimes, sensors = get_band(anuga_break_times[0]-TIME_BORDER,
                                    100, dtimes, sensors, 0)
-        locations, sensors = get_band(run_data['break_xs'][0]-LOCATION_BORDER,
-                                      100, locations, sensors, -1)
+        locations, sensors = get_band(
+            min(run_data['break_xs'])-LOCATION_BORDER,
+            100, locations, sensors, -1)
         #print "dtimes", dtimes
         #print "sensors", sensors
@@ -241,6 +249,35 @@
                      is_interactive=is_interactive)
         
- 
-def auto_find_min_slopes(outputdir_tag, scenarios):
+def find_froude(times_froude, locations_froude, froudes_array,
+                times, locations):
+    if len(times) == 0:
+        return []
+    time_indexes = searchsorted(times_froude, times)
+    location_indexes = searchsorted(locations_froude, locations)
+
+    
+    assert len(time_indexes) == len(location_indexes)
+
+    froudes = []
+    for time_i, loc_i, time, location in map(None, time_indexes,
+                                             location_indexes,
+                                             times, locations):
+        # the time values should be the same
+        assert times_froude[time_i] == time
+
+        # The distance value should be half way between the froude locations
+        midpoint = locations_froude[loc_i-1] + \
+                   (locations_froude[loc_i]-locations_froude[loc_i-1])*0.5
+        #print "location", location
+        #print "midpoint", midpoint
+        assert location == midpoint
+        froude = froudes_array[time_i, loc_i-1] + \
+                   (froudes_array[time_i, loc_i]- \
+                    froudes_array[time_i, loc_i-1])*0.5
+        froudes.append(froude)
+                
+    return froudes
+    
+def auto_find_min_slopes(slope_tag, outputdir_tag, scenarios):
     """
     
@@ -258,32 +295,88 @@
                 break_time -  run_data['ANUGA_start_time'])
             
-        stage_file = pro_instance.outputdir + "fslope_stage_" + end
+        stage_file = pro_instance.outputdir + slope_tag + "slope_stage_" + end
+        froude_file = pro_instance.outputdir + slope_tag + "slope_froude_" + \
+                      end
         
         times, slope_locations, slopes = load_slopes(stage_file)
+        #print "slope_locations", slope_locations
+        times_froude, locations_froude, froudes_a = load_sensors(froude_file)
+        #print "locations_froude", locations_froude
         waves = find_min_slopes(times, slope_locations, slopes,
-                                run_data['break_times'],
                                 anuga_break_times,
                                 run_data['band_offset'])
+        
         # write the wave info here
-        # return the keys actually.
-        for wave_key in waves.iterkeys():
-            wave_file = stage_file[:-3] + '_'+ wave_key + ".csv"
+        # and find the froude values
+        for i, wave in enumerate(waves):
+            
+            id = "wave_" + str(i)
+            wave_file = stage_file[:-4] + '_'+ id + ".csv"
             print "wave_file", wave_file
+            froudes = find_froude(times_froude, locations_froude,
+                             froudes_a, wave[TIME_STR],
+                                  slope_locations)
             wave_writer = writer(file(wave_file, "wb"))
-            wave_writer.writerow(["x location", "min slope", "Time"])
+            wave_writer.writerow(["x location", "min slope", "Time", "Froude"])
             wave_writer.writerows(map(None,
                                       slope_locations,
-                                      waves[wave_key][SLOPE_STR],
-                                      waves[wave_key][TIME_STR]))
+                                      wave[SLOPE_STR],
+                                      wave[TIME_STR],
+                                      froudes))
+
             
-            #wave_writer.close()
+def auto_plot_froude_slopes(slope_tag, outputdir_tag, scenarios):
+    """
+    
+    """
+    
+    plot_type = ".pdf"
+    
+    for run_data in scenarios:
+        id = run_data['scenario_id']
+        outputdir_name = id + outputdir_tag
+        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
+                                       outputdir_name=outputdir_name)
+
+        assert len(run_data['break_times']) == len(run_data['break_xs'])
+        assert len(run_data['break_times']) == len(run_data['break_type'])
+        
+        end = id + ".csv"
+        
+        anuga_break_times = []
+        for break_time in run_data['break_times']:
+            anuga_break_times.append( \
+                break_time -  run_data['ANUGA_start_time'])
             
+        #run_data['break_type'] = (run_data['break_type'][0])
+        for i in range(len(run_data['break_type'])):
             
-        
-    
-def gauges_for_slope(outputdir_tag, scenarios):
+            wave_tag = "wave_" + str(i)
+            stage_file = pro_instance.outputdir + slope_tag + \
+                         "slope_stage_" + end
+            wave_file = stage_file[:-4] + '_'+ wave_tag + ".csv"
+            save_as = pro_instance.plots_dir + sep + \
+                            outputdir_name + "_" + wave_tag + plot_type
+            print "wave_file", wave_file
+            break_type = run_data['break_type'][i]
+            plot_title = id  + ' Wave: ' + str(i) + \
+                        " Break Type: " + break_type + '\n' + \
+                        "File: " + wave_file[34:] # not good!
+            plot_foude_slope_stage(wave_file,
+                                   anuga_break_times[i],
+                                   run_data['break_xs'][i],
+                                   plot_title=plot_title,
+                                   break_type=break_type,
+                                   save_as=save_as,
+                                   is_interactive=False)
+        
+        
+    
+def gauges_for_slope(slope_tag, outputdir_tag, scenarios):
     """
     This is used to create a stage file, using gauges relivent to
     finding a slope.
+
+    It also create's a frounde file.
     """
     dx = 0.05
@@ -306,22 +399,20 @@
             pro_instance.outputdir + basename +".sww",
             points,
-            pro_instance.outputdir + "fslope_depth_" + end,
-            pro_instance.outputdir + "fslope_velocity_x_" + end,
-            pro_instance.outputdir + "fslope_velocity_y_" + end,
-            pro_instance.outputdir + "fslope_stage_" + end,
-            pro_instance.outputdir + "fslope_froude_" + end,
+            pro_instance.outputdir + slope_tag + "slope_depth_" + end,
+            pro_instance.outputdir + slope_tag + "slope_velocity_x_" + end,
+            pro_instance.outputdir + slope_tag + "slope_velocity_y_" + end,
+            pro_instance.outputdir + slope_tag + "slope_stage_" + end,
+            pro_instance.outputdir + slope_tag + "slope_froude_" + end,
             time_thinning=1)
 
 def find_min_slopes(times, slope_locations, slopes,
-                    break_times, anuga_break_times, band_offset):
-    bands = break_times2bands(break_times, band_offset)
-
-    waves = {}
+                    anuga_break_times, band_offset):
+    bands = break_times2bands(anuga_break_times, band_offset)
+
+    waves = []
     for i,_ in enumerate(bands[0:-1]):
-        id = "wave " + str(i)
         max_q, max_q_times = get_min_in_band(bands[i], bands[i+1],
                                              times, slopes)
-        waves[id] = {SLOPE_STR:max_q,
-                     TIME_STR:max_q_times}
+        waves.append({SLOPE_STR:max_q, TIME_STR:max_q_times})
     return waves
 
@@ -366,5 +457,12 @@
                                          time_vector, quantity_array, 0)
     #print "band_quantity",band_quantity
-    max_quantity_indices = argmin(band_quantity, axis=0)
+    try:
+        max_quantity_indices = argmin(band_quantity, axis=0)
+    except:
+        #print "time_vector", time_vector
+        print "min_time",min_time 
+        print "max_time", max_time
+        return [],[]
+        
     #print "max_quantity_indices", max_quantity_indices
     max_quantity_times = choose(max_quantity_indices, band_time)
@@ -386,4 +484,5 @@
     bands.append(break_times[0]-0.5*(break_times[1]-break_times[0]))
 
+
     for i,break_x in enumerate(break_times[0:-1]):
         bands.append(break_times[i]+0.5*(break_times[i+1]-break_times[i]))
@@ -394,4 +493,92 @@
     return bands
 
+def plot_foude_slope_stage(wave_file,
+                           break_time,
+                           break_x,
+                           save_as=None,
+                           plot_title="",
+                           is_interactive=False,
+                           break_type="",
+                           use_axis=None):
+    """ 
+    """
+    from pylab import ion, plot, xlabel, ylabel, close, legend, \
+         savefig, title, axis, setp, subplot, grid, axvspan
+    from anuga.shallow_water.data_manager import csv2dict
+
+
+
+    # Load in the csv files and convert info from strings to floats
+    simulation, _ = csv2dict(wave_file)
+    location = [float(x) for x in simulation['x location']]
+    slope = [float(x) for x in simulation['min slope']]
+    time = [float(x) for x in simulation['Time']]
+    froude = [float(x) for x in simulation['Froude']]
+
+    min_location = min(location)
+    max_location = max(location)
+    
+    if is_interactive:
+        ion()
+    # The upper subplot
+    subplot(311)
+    l_froude = plot(location, froude)
+    #setp(l_froude, color='r')
+
+    # Add axis stuff
+    title(plot_title)
+    y_label = "Froude Number"
+    ylabel(y_label)
+    grid(True)
+    axvspan(break_x-0.001,break_x+0.001, facecolor='g')
+    axis(ymin=0, ymax=1.8)
+    
+    # The slope subplot
+    subplot(312)
+    l_slope = plot(location, slope)
+    setp(l_slope, color='r')
+
+    # Add axis stuff and legend
+    x_label = "X location, m"
+    y_label = "Stage slope"
+    #xlabel(x_label)
+    ylabel(y_label)
+    grid(True)
+    axvspan(break_x-0.001,break_x+0.001, facecolor='g')
+    axis(ymin=-0.5, ymax=0)
+
+    # The time, x location subplot
+    subplot(313)
+    l_time = plot(location, time)
+    setp(l_time, color='g')
+    #print "break_x", break_x
+    #print "break_time", break_time
+    plot([break_x], [break_time], 'yo')
+    #plot([break_x-1], [], 'yo')
+
+    # Add axis stuff and legend
+    x_label = "X location, m"
+    y_label = "time, sec"
+    xlabel(x_label)
+    ylabel(y_label)
+    grid(True)
+
+    
+    # The order defines the label
+    #legend((legend_exp, legend_sim),'upper left')
+    #legend(('Wave front'),'upper left')
+    if use_axis is not None:
+        axis(use_axis)
+    
+    if is_interactive:
+        # Wait for enter pressed
+        raw_input()
+
+    if save_as is not None:
+        savefig(save_as)
+    
+    #Need to close this plot
+    close()
+    
 #-------------------------------------------------------------
 if __name__ == "__main__":
@@ -400,9 +587,15 @@
     from scenarios import scenarios
     #scenarios = [scenarios[0]]
-    outputdir_tag = "_good_tri_area_0.01_limiterC"
+    outputdir_tag = "_good_tri_area_0.01_limiterD"
+    slope_tag = ""
+    #outputdir_tag = "_good_tri_area_0.01_limiterC"
+    #slope_tag = "f"
     #outputdir_tag = "_test_limiterC"
     #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
-    #gauges_for_slope(outputdir_tag, scenarios)
-    auto_graph_slopes(outputdir_tag, scenarios) #, is_interactive=True) 
-    #auto_find_min_slopes(outputdir_tag, scenarios)
-    auto_graph_froudes(outputdir_tag, scenarios) 
+    #scenarios = scenarios[4:] # !!!!!!!!!!!!!!!!!!!!!!
+    
+    #gauges_for_slope(slope_tag, outputdir_tag, scenarios)
+    #auto_graph_slopes(outputdir_tag, scenarios) #, is_interactive=True) 
+    #auto_find_min_slopes(slope_tag, outputdir_tag, scenarios)
+    #auto_graph_froudes(outputdir_tag, scenarios) 
+    auto_plot_froude_slopes(slope_tag, outputdir_tag, scenarios)
Index: /anuga_work/development/Hinwood_2008/test_slope.py
===================================================================
--- /anuga_work/development/Hinwood_2008/test_slope.py	(revision 5502)
+++ /anuga_work/development/Hinwood_2008/test_slope.py	(revision 5503)
@@ -38,5 +38,5 @@
         actual = array([0,0.01,0.02,0.03])
         self.assert_ (allclose(dtimes, actual, 0.001))
-        actual = array([0.5, 0.501, 0.502])
+        actual = array([0.5, 0.5005, 0.5015])
         self.assert_ (allclose(depths, actual, 0.001))
 
@@ -111,5 +111,5 @@
         self.assert_ (allclose(ensure_numeric(bands), actual, 0.001))
         
-    def test_get_time_band(self):
+    def not_used_test_get_time_band(self):
         time = array([0,1,2,3,4,5])
         slopes = array([[10,12],
@@ -158,9 +158,23 @@
                         [10,20,30,40]])
         self.assert_ (allclose(band_quantity, actual, 0.001))      
+
+    def test_find_froude(self):
+        times_froude = array([0,1])
+        locations_froude = array([0,1,2,3,4,5])
+        froudes = array([[0,1,2,3,4,5],
+                        [0,10,20,30,40,50]])
+        times = array([0,0,1,1])
+        locations = array([0.5, 1.5, 2.5, 3.5])
+        froudes = find_froude(times_froude, locations_froude,
+                             froudes, times, locations)
+        actual = array([0.5, 1.5, 25, 35])
+        self.assert_ (allclose(ensure_numeric(froudes), actual, 0.001))
+        
+        
         
 if __name__ == "__main__":
 
     suite = unittest.makeSuite(slopeTestCase,'test')
-    suite = unittest.makeSuite(slopeTestCase,'test_get_band')
+    #suite = unittest.makeSuite(slopeTestCase,'test_find_froude')
 
     runner = unittest.TextTestRunner()  #verbosity=2)
