Index: /documentation/experimentation/boundary_ANUGA_MOST/report/50100MOSTcomparison_onslow.tex
===================================================================
--- /documentation/experimentation/boundary_ANUGA_MOST/report/50100MOSTcomparison_onslow.tex	(revision 3190)
+++ /documentation/experimentation/boundary_ANUGA_MOST/report/50100MOSTcomparison_onslow.tex	(revision 3190)
@@ -0,0 +1,27 @@
+\begin{table} \label{table:locationsonslow} 
+\caption{Defined point locations for comparison study for Onslow region.}
+\begin{center}
+\begin{tabular}{|l|l|l|l|}\hline
+\bf{Point Name} & \bf{Easting} & \bf{Northing} & \bf{Elevation (m)}\\ \hline
+Region polygon1 & 305000.00 & 7601000.00 & 6.55 \\ \hline 
+Onshore Centre & 305000.00 & 7605000.00 & 1.42 \\ \hline 
+Ocean polygon 2 & 305000.00 & 7646000.00 & -30.83 \\ \hline 
+Ocean polygon 1 & 305000.00 & 7624000.00 & -3.18 \\ \hline 
+Boundary 1 & 305000.00 & 7670000.00 & -49.04 \\ \hline 
+Western region & 292000.00 & 7611000.00 & -9.03 \\ \hline 
+Boundary 2 & 280000.00 & 7624000.00 & -11.34 \\ \hline 
+Western corner & 265000.00 & 7620000.00 & -31.18 \\ \hline 
+Centre & 292000.00 & 7630000.00 & -9.48 \\ \hline 
+Region polygon2 & 293000.00 & 7602500.00 & -9.69 \\ \hline 
+\end{tabular} 
+  \end{center} 
+ \end{table} 
+ 
+\begin{figure}[hbt] 
+ \centerline{ \includegraphics[width=150mm, height=100mm]{../report_figures/onslow_boundary_gauges.png}}
+\caption{Point locations used for comparison between boundary at 50m and 100m contour for Onslow region.}  
+\label{fig:points}
+\end{figure}
+\input{latexoutput_compare50_100_MOST} 
+ \clearpage 
+ 
Index: /documentation/experimentation/boundary_ANUGA_MOST/report/MOST_ANUGA.tex
===================================================================
--- /documentation/experimentation/boundary_ANUGA_MOST/report/MOST_ANUGA.tex	(revision 3189)
+++ /documentation/experimentation/boundary_ANUGA_MOST/report/MOST_ANUGA.tex	(revision 3190)
@@ -226,4 +226,5 @@
 see for example the output for the Ocean polygon 1 and 2 locations.
 
+\input{50100MOSTcomparison_onslow}
 
 It is more instructive in this case to compare differences in
Index: /documentation/experimentation/boundary_ANUGA_MOST/report/latexoutput_compare50_100_MOST.tex
===================================================================
--- /documentation/experimentation/boundary_ANUGA_MOST/report/latexoutput_compare50_100_MOST.tex	(revision 3190)
+++ /documentation/experimentation/boundary_ANUGA_MOST/report/latexoutput_compare50_100_MOST.tex	(revision 3190)
@@ -0,0 +1,100 @@
+\begin{figure}[hbt] 
+ \centering 
+ \begin{tabular}{cc} 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeRegionpolygon1stage.png}& 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeRegionpolygon1speed.png}\\ 
+\end{tabular} 
+ \caption{Time series for stage and speed at Region polygon1 gauge location (elevation 6.55m)} 
+ \label{fig:gaugeRegionpolygon1} 
+ \end{figure} 
+ 
+\begin{figure}[hbt] 
+ \centering 
+ \begin{tabular}{cc} 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeOnshoreCentrestage.png}& 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeOnshoreCentrespeed.png}\\ 
+\end{tabular} 
+ \caption{Time series for stage and speed at Onshore Centre gauge location (elevation 1.42m)} 
+ \label{fig:gaugeOnshoreCentre} 
+ \end{figure} 
+ 
+\begin{figure}[hbt] 
+ \centering 
+ \begin{tabular}{cc} 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeOceanpolygon2stage.png}& 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeOceanpolygon2speed.png}\\ 
+\end{tabular} 
+ \caption{Time series for stage and speed at Ocean polygon 2 gauge location (elevation -30.83m)} 
+ \label{fig:gaugeOceanpolygon2} 
+ \end{figure} 
+ 
+\begin{figure}[hbt] 
+ \centering 
+ \begin{tabular}{cc} 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeOceanpolygon1stage.png}& 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeOceanpolygon1speed.png}\\ 
+\end{tabular} 
+ \caption{Time series for stage and speed at Ocean polygon 1 gauge location (elevation -3.18m)} 
+ \label{fig:gaugeOceanpolygon1} 
+ \end{figure} 
+ 
+\begin{figure}[hbt] 
+ \centering 
+ \begin{tabular}{cc} 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeBoundary1stage.png}& 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeBoundary1speed.png}\\ 
+\end{tabular} 
+ \caption{Time series for stage and speed at Boundary 1 gauge location (elevation -49.04m)} 
+ \label{fig:gaugeBoundary1} 
+ \end{figure} 
+ 
+\begin{figure}[hbt] 
+ \centering 
+ \begin{tabular}{cc} 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeWesternregionstage.png}& 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeWesternregionspeed.png}\\ 
+\end{tabular} 
+ \caption{Time series for stage and speed at Western region gauge location (elevation -9.03m)} 
+ \label{fig:gaugeWesternregion} 
+ \end{figure} 
+ 
+\begin{figure}[hbt] 
+ \centering 
+ \begin{tabular}{cc} 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeBoundary2stage.png}& 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeBoundary2speed.png}\\ 
+\end{tabular} 
+ \caption{Time series for stage and speed at Boundary 2 gauge location (elevation -11.34m)} 
+ \label{fig:gaugeBoundary2} 
+ \end{figure} 
+ 
+\begin{figure}[hbt] 
+ \centering 
+ \begin{tabular}{cc} 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeWesterncornerstage.png}& 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeWesterncornerspeed.png}\\ 
+\end{tabular} 
+ \caption{Time series for stage and speed at Western corner gauge location (elevation -31.18m)} 
+ \label{fig:gaugeWesterncorner} 
+ \end{figure} 
+ 
+\begin{figure}[hbt] 
+ \centering 
+ \begin{tabular}{cc} 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeCentrestage.png}& 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeCentrespeed.png}\\ 
+\end{tabular} 
+ \caption{Time series for stage and speed at Centre gauge location (elevation -9.48m)} 
+ \label{fig:gaugeCentre} 
+ \end{figure} 
+ 
+\begin{figure}[hbt] 
+ \centering 
+ \begin{tabular}{cc} 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeRegionpolygon2stage.png}& 
+\includegraphics[width=0.49\linewidth, height=50mm]{../report_figures/gaugeRegionpolygon2speed.png}\\ 
+\end{tabular} 
+ \caption{Time series for stage and speed at Region polygon2 gauge location (elevation -9.69m)} 
+ \label{fig:gaugeRegionpolygon2} 
+ \end{figure} 
+ 
Index: /documentation/user_manual/examples/project.py
===================================================================
--- /documentation/user_manual/examples/project.py	(revision 3189)
+++ /documentation/user_manual/examples/project.py	(revision 3190)
@@ -83,7 +83,8 @@
 # demo poly which isn't rectangular
 j0 = [385000, 6280000]
-j1 = [350000, 6270000]
-j2 = [325000, 6263000]
-j3 = [325000, 6260000]
+j1 = [360000, 6272500]
+j2 = [335000, 6272500]
+j3 = [330000, 6265000]
+j31 = [325000, 6260000]
 j4 = [316000, 6260000]
 j5 = [316000, 6247000]
@@ -91,12 +92,16 @@
 j7 = [385000, 6238000]
 
-demopoly = [j0, j2, j3, j4, j5, j6, j7]
+demopoly = [j0, j1, j2, j3, j31, j4, j5, j6, j7]
 
+polygonptsfile4 = polygondir + 'poly1'
+polygonptsfile0 = polygondir + 'poly2'
 polygonptsfile1 = polygondir + 'poly3'
 polygonptsfile2 = polygondir + 'poly4'
 polygonptsfile3 = polygondir + 'poly5'
-northern_polygon = read_polygon(polygonptsfile1 + '.csv')
+northern_polygon = read_polygon(polygonptsfile0 + '.csv')
+manly_polygon = read_polygon(polygonptsfile1 + '.csv')
 harbour_polygon = read_polygon(polygonptsfile2 + '.csv')
 southern_polygon = read_polygon(polygonptsfile3 + '.csv')
+top_polygon = read_polygon(polygonptsfile4 + '.csv')
 
 #Interior regions - the Harbour
Index: /documentation/user_manual/examples/runsydney.py
===================================================================
--- /documentation/user_manual/examples/runsydney.py	(revision 3189)
+++ /documentation/user_manual/examples/runsydney.py	(revision 3190)
@@ -57,15 +57,19 @@
 interior_regions = [[project.northern_polygon, interior_res],
                     [project.harbour_polygon, interior_res],
-                    [project.southern_polygon, interior_res]]
+                    [project.southern_polygon, interior_res],
+                    [project.manly_polygon, interior_res],
+                    [project.top_polygon, interior_res]]
 
 
 create_mesh_from_regions(project.demopoly,
                          boundary_tags= {'oceannorth': [0], 
-                                         'onshorewest': [1],
-                                         'harbournorth': [2],
-                                         'harbourwest': [3],
-                                         'harboursouth': [4],
-                                         'oceansouth': [5],
-                                         'oceaneast': [6]},
+                                         'onshorenorth1': [1],
+                                         'onshorenorthwest1': [2],
+                                         'onshorenorthwest2': [3],
+                                         'onshorenorth2': [4],
+                                         'onshorewest': [5],
+                                         'onshoresouth': [6],
+                                         'oceansouth': [7],
+                                         'oceaneast': [8]},
                          maximum_triangle_area=100000,
                          filename=meshname,           
@@ -115,5 +119,4 @@
                     verbose = True)
 
-
 #------------------------------------------------------------------------------
 # Setup boundary conditions (all Dirichlet)
@@ -123,9 +126,13 @@
 
 Bd = Dirichlet_boundary([0.0,0.0,0.0])
-domain.set_boundary( {'oceannorth': Bd, 'onshorewest': Bd,
-                      'harbournorth': Bd, 'harbourwest': Bd,
-                      'harboursouth': Bd, 'oceansouth': Bd,
-                      'oceaneast': Bd})
-
+domain.set_boundary( {'oceannorth': Bd,
+                      'onshorenorth1': Bd,
+                      'onshorenorthwest1': Bd,
+                      'onshorenorthwest2': Bd,
+                      'onshorenorth2': Bd,
+                      'onshorewest': Bd,
+                      'onshoresouth': Bd,
+                      'oceansouth': Bd,
+                      'oceaneast': Bd} )
 
 #------------------------------------------------------------------------------
@@ -135,8 +142,6 @@
 import time
 t0 = time.time()
-tend = 7200.0
-print 'Simulate to time = %d secs' %tend
 
-for t in domain.evolve(yieldstep = 60, duration = tend):
+for t in domain.evolve(yieldstep = 60, finaltime = 7200):
     print domain.timestepping_statistics()
     print domain.boundary_statistics(tags = 'oceaneast')    
Index: /inundation/pyvolution/smf.py
===================================================================
--- /inundation/pyvolution/smf.py	(revision 3189)
+++ /inundation/pyvolution/smf.py	(revision 3190)
@@ -346,6 +346,7 @@
         kappa = self.kappa
         kappad = self.kappad
-        #amin = self.find_min(x0,wa,kappad,kappa,dx,am)
-        amin = 1.0
+        amin = self.find_min(x0,wa,kappad,kappa,dx,am)
+        print 'hello amin', amin
+        #amin = 1.0
 
         #double Gaussian calculation assumes water displacement is oriented
@@ -369,7 +370,11 @@
                 if z[i] > maxz: maxz = z[i]
                 if z[i] < minz: minz = z[i]
+                
             except OverflowError:
                 pass
-        
+
+        print 'max z', maxz
+        print 'min z', minz
+                
         return z
 
@@ -407,4 +412,5 @@
         count_max = 1000000
         c = 0
+        am = 1.0
 
         while c < count_max and deriv > 0:
Index: /inundation/pyvolution/util.py
===================================================================
--- /inundation/pyvolution/util.py	(revision 3189)
+++ /inundation/pyvolution/util.py	(revision 3190)
@@ -548,4 +548,5 @@
                    reportname = None,
                    plot_quantity = None,
+                   surface = None,
                    time_min = None,
                    time_max = None,
@@ -596,4 +597,7 @@
                     - default will be ['stage', 'speed', 'bearing']
 
+    surface          - if True, then generate solution surface with 3d plot
+                       and save to current working directory
+    
     time_min         - beginning of user defined time range for plotting purposes
                         - default will be first available time found in swwfile
@@ -631,4 +635,5 @@
                         reportname,
                         plot_quantity,
+                        surface,
                         time_min,
                         time_max,
@@ -644,4 +649,5 @@
                     reportname = None,
                     plot_quantity = None,
+                    surface = None,
                     time_min = None,
                     time_max = None,
@@ -670,4 +676,7 @@
         check_list(plot_quantity)
 
+    if surface is None:
+        surface = False
+        
     if title_on is None:
         title_on = True
@@ -727,6 +736,6 @@
     if verbose: print 'Inputs OK - going to generate figures'
 
-    return generate_figures(plot_quantity, file_loc, report, reportname, leg_label,
-                            f_list, gauges, locations, elev, production_dirs,
+    return generate_figures(plot_quantity, file_loc, report, reportname, surface,
+                            leg_label, f_list, gauges, locations, elev, production_dirs,
                             time_min, time_max, title_on, label_id, verbose)
                          
@@ -789,13 +798,17 @@
     return bearing
 
-def generate_figures(plot_quantity, file_loc, report, reportname, leg_label, f_list, gauges,
-                     locations, elev, production_dirs, time_min, time_max,
+def generate_figures(plot_quantity, file_loc, report, reportname, surface,
+                     leg_label, f_list,
+                     gauges, locations, elev, production_dirs, time_min, time_max,
                      title_on, label_id, verbose):
 
     from math import sqrt, atan, degrees
-    from Numeric import ones, allclose, zeros, Float
+    from Numeric import ones, allclose, zeros, Float, ravel
     from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
     from pylab import ion, hold, plot, axis, figure, legend, savefig, \
          xlabel, ylabel, title, close, subplot
+
+    import pylab as p1
+    import mpl3d.mplot3d as p3
 
     if report == True:    
@@ -837,4 +850,5 @@
     bearings = zeros((n0,m,p), Float)
     depths = zeros((n0,m,p), Float)
+    eastings = zeros((n0,m,p), Float)
     min_stages = []
     max_stages = []
@@ -842,4 +856,7 @@
     max_speeds = []
     c = 0
+    model_time_plot3d = zeros((n0,m), Float)
+    stages_plot3d = zeros((n0,m), Float)
+    eastings_plot3d = zeros((n0,m),Float)
     ##### loop over each swwfile #####
     for j, f in enumerate(f_list):
@@ -858,5 +875,4 @@
             fid_out.write(s)
             #### generate quantities #######
-            
             for i, t in enumerate(f.get_time()):
                 if time_min <= t <= time_max:
@@ -872,6 +888,6 @@
                         vel = m / (depth + 1.e-30) 
                     bearing = calc_bearing(uh, vh)
-                    model_time[i,k,j] = t/60.0         
-                    stages[i,k,j] = w 
+                    model_time[i,k,j] = t/60.0
+                    stages[i,k,j] = w
                     elevations[i,k,j] = z 
                     xmom[i,k,j] = uh 
@@ -880,5 +896,7 @@
                     speed[i,k,j] = vel 
                     bearings[i,k,j] = bearing 
-                    depths[i,k,j] = depth 
+                    depths[i,k,j] = depth
+                    thisgauge = gauges[k]
+                    eastings[i,k,j] = thisgauge[0]
                     s = '%.2f, %.2f, %.2f, %.2f\n' %(t, w, m, vel)
                     fid_out.write(s)
@@ -886,15 +904,39 @@
                     if w < min_stage: min_stage = w
                     if m > max_momentum: max_momentum = m
-                    if vel > max_speed: max_speed = vel
-
+                    if vel > max_speed: max_speed = vel            
+            
             max_stages.append(max_stage)
             min_stages.append(min_stage)
             max_momentums.append(max_momentum)
             max_speeds.append(max_speed)        
+            #### finished generating quantities for each swwfile #####
+
+        model_time_plot3d[:,:] = model_time[:,:,j]
+        stages_plot3d[:,:] = stages[:,:,j]
+        eastings_plot3d[:,] = eastings[:,:,j]
             
-            #### finished generating quantities for each swwfile #####
-            
+        if surface ==  True:
+            print 'Printing surface figure'
+            for i in range(2):
+                fig = p1.figure(10)
+                ax = p3.Axes3D(fig)
+                if len(gauges) > 80:
+                    ax.plot_surface(eastings,model_time_plot3d,stages_plot3d)
+                    #ax.plot_surface(eastings[:,:,j],model_time[:,:,j],stages[:,:,j])
+                else:
+                    #ax.plot_wireframe(eastings[:,:,j],model_time[:,:,j],stages[:,:,j])
+                    ax.plot_wireframe(eastings,model_time_plot3d,stages_plot3d)
+                    #ax.plot3D(ravel(eastings[:,:,j]),ravel(model_time[:,:,j]),ravel(stages[:,:,j]))
+                ax.set_xlabel('x')
+                ax.set_ylabel('time')
+                ax.set_zlabel('stage')
+                fig.add_axes(ax)
+                p1.show()
+                surfacefig = 'solution_surface%s' %leg_label[j]
+                p1.savefig(surfacefig)
+                p1.close()
+        
     #### finished generating quantities for all swwfiles #####
-            
+
     stage_axis = axis([time_min/60.0, time_max/60.0, min(min_stages), max(max_stages)*1.1])
     vel_axis = axis([time_min/60.0, time_max/60.0, min(max_speeds), max(max_speeds)*1.1])
@@ -975,8 +1017,6 @@
                         mkdir (figdir)
                     latex_file_loc = figdir.replace(sep,altsep) 
-                    # storing files in production directory
-                    graphname_latex = '%sgauge%s%s' %(latex_file_loc, gaugeloc2, which_quantity)     
-                    # giving location in latex output file
-                    graphname_report_input = '%sgauge%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity)
+                    graphname_latex = '%sgauge%s%s' %(latex_file_loc, gaugeloc2, which_quantity) # storing files in production directory    
+                    graphname_report_input = '%sgauge%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity) # giving location in latex output file
                     graphname_report.append(graphname_report_input)
                     
@@ -991,8 +1031,6 @@
 
                     if len(label_id) == 1: 
-                        # storing files in production directory
-                        graphname_latex = '%sgauge%s%s%s' %(latex_file_loc, gaugeloc2, which_quantity, label_id2)     
-                        # giving location in latex output file
-                        graphname_report = '%sgauge%s%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity, label_id2)
+                        graphname_latex = '%sgauge%s%s%s' %(latex_file_loc, gaugeloc2, which_quantity, label_id2) # storing files in production directory  
+                        graphname_report = '%sgauge%s%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity, label_id2) # giving location in latex output file
                         s = '\includegraphics[width=0.49\linewidth, height=50mm]{%s%s}' %(graphname_report, '.png')
                         fid.write(s)
Index: /production/onslow_2006/MOST_timeseries.py
===================================================================
--- /production/onslow_2006/MOST_timeseries.py	(revision 3189)
+++ /production/onslow_2006/MOST_timeseries.py	(revision 3190)
@@ -61,5 +61,5 @@
 
 # Start report generation
-input_name = project.comparereportdir + sep + 'comparison_onslow.tex'
+input_name = project.comparereportdir + sep + 'MOST_ANUGA_comparison_onslow.tex'
 fid = open(input_name, 'w')
 
Index: /production/onslow_2006/compare_timeseries.py
===================================================================
--- /production/onslow_2006/compare_timeseries.py	(revision 3189)
+++ /production/onslow_2006/compare_timeseries.py	(revision 3190)
@@ -23,5 +23,6 @@
 # User defined inputs
 production_dirs = {'20060515_001733': '100m boundary',
-                   '20060530_102753': '50m boundary'}
+                   '20060530_102753': '50m boundary',
+                   'MOST': 'MOST'}
 
 gauge_map = 'onslow_boundary_gauges.png'
@@ -36,4 +37,6 @@
     file_loc = project.outputdir + label_id + sep
     swwfile = file_loc + project.basename + '.sww'
+    if label_id == 'MOST':
+        swwfile = project.boundarydir + project.boundary_basename + '.sww'
     swwfiles[swwfile] = label_id
 
@@ -42,5 +45,5 @@
                               production_dirs,
                               report = True,
-                              reportname = 'latexoutput_compare',
+                              reportname = 'latexoutput_compare50_100_MOST',
                               surface = False,
                               plot_quantity = plot_quantity,
@@ -57,5 +60,5 @@
 
 # Start report generation
-input_name = project.comparereportdir + sep + 'comparison_onslow.tex'
+input_name = project.comparereportdir + sep + '50100MOSTcomparison_onslow.tex'
 fid = open(input_name, 'w')
 
Index: /production/onslow_2006/make_report.py
===================================================================
--- /production/onslow_2006/make_report.py	(revision 3189)
+++ /production/onslow_2006/make_report.py	(revision 3190)
@@ -20,4 +20,5 @@
 The report structure is
 
+* Executive Summary
 * Introduction
 * Modelling Methodology
@@ -159,10 +160,10 @@
 \\begin{document}
   \maketitle
-  
-  \\begin{abstract}
-    \input{abstract}
-  \end{abstract}
-  
+   
   \\tableofcontents
+
+    \section{Executive Summary}
+    \label{sec:execsum}
+  \input{execsum}
   
   \section{Introduction}
Index: /production/pt_hedland_2006/compare_timeseries.py
===================================================================
--- /production/pt_hedland_2006/compare_timeseries.py	(revision 3190)
+++ /production/pt_hedland_2006/compare_timeseries.py	(revision 3190)
@@ -0,0 +1,97 @@
+"""
+Generate latex output comparing time series of given sww files
+
+Inputs:
+
+output dirs: dictionary of output directories with a
+             association to simulation, eg 50m boundary, 100m boundary
+                    
+Outputs:
+
+* latex output to generate figures for report 
+* figures used for report stored in the report_figure directory
+NOTE, this directory will need to be committed, as well as
+the latex files.
+
+June 2006                    
+"""
+
+from os import sep
+import project
+from pyvolution.util import sww2timeseries, get_gauges_from_file
+
+# User defined inputs
+production_dirs = {'': '100m boundary',
+                   '': '50m boundary'}
+
+gauge_map = 'pt_hedland_boundary_gauges.png'
+
+plot_quantity = ['stage', 'speed']
+
+# Create sections and graphs for each designated production directory
+latex_output = []
+swwfiles = {}
+for label_id in production_dirs.keys():
+    
+    file_loc = project.outputdir + label_id + sep
+    swwfile = file_loc + project.basename + '.sww'
+    swwfiles[swwfile] = label_id
+
+texname, vec = sww2timeseries(swwfiles,
+                              project.gauges50,
+                              production_dirs,
+                              report = True,
+                              reportname = 'latexoutput_compare_pt_hedland',
+                              plot_quantity = plot_quantity,
+                              time_min = None,
+                              time_max = None,
+                              title_on = False,
+                              verbose = True)
+
+latex_output.append(texname)
+
+from shutil import copy, move
+copy ('report' + sep + texname + '.tex', project.comparereportdir + sep + texname + '.tex')
+
+
+# Start report generation
+input_name = project.comparereportdir + sep + 'comparison_pt_hedland.tex'
+fid = open(input_name, 'w')
+
+# Generate latex output for location points
+s = '\\begin{table} \label{table:locationspthedland} \n'
+fid.write(s)
+s = '\caption{Defined point locations for comparison study for Pt Hedland region.}' 
+fid.write(s)
+s = """
+\\begin{center}
+\\begin{tabular}{|l|l|l|l|}\hline
+\\bf{Point Name} & \\bf{Easting} & \\bf{Northing} & \\bf{Elevation (m)}\\\\ \hline
+"""
+fid.write(s)
+
+for i, thisvec in enumerate(vec):
+    name = thisvec[0]
+    east = thisvec[1]
+    north = thisvec[2]
+    elev = thisvec[3]
+    s = '%s & %.2f & %.2f & %.2f \\\\ \hline \n' %(name, east, north, elev)
+    fid.write(s)
+    for which_one in plot_quantity:
+        move ('report_figures' + sep + 'gauge'+name.replace(' ','')+which_one+'.png', project.comparereportfigdir)
+
+s = '\\end{tabular} \n  \end{center} \n \end{table} \n \n'
+fid.write(s)
+
+s = '\\begin{figure}[hbt] \n \centerline{ \includegraphics[width=150mm, height=100mm]{../report_figures/%s}}' %gauge_map
+fid.write(s)
+
+s  = """
+\caption{Point locations used for comparison between boundary at 50m and 100m contour for Onslow region.}  
+\label{fig:points}
+\end{figure}
+"""
+fid.write(s)
+    
+s = '\input{%s} \n \clearpage \n \n' %latex_output[0]
+fid.write(s)
Index: /production/pt_hedland_2006/make_report.py
===================================================================
--- /production/pt_hedland_2006/make_report.py	(revision 3189)
+++ /production/pt_hedland_2006/make_report.py	(revision 3190)
@@ -20,5 +20,7 @@
 The report structure is
 
+* Executive Summary
 * Introduction
+* Modelling Methodology
 * Tsunami scenario
 * Inundation model
@@ -68,15 +70,11 @@
 report_title = 'Tsunami impact modelling for the North West shelf: %s' %scenario_name.title()
 
-production_dirs = {'': '1.5 AHD',
-                   '': '-1.5 AHD',
+production_dirs = {'': '3.6 AHD',
+                   '': '-3.9 AHD',
                    '': '0 AHD'}
 
-max_maps = {'1.5 AHD': 'HAT_map',
-            '-1.5 AHD': 'LAT_map',
+max_maps = {'3.6 AHD': 'HAT_map',
+            '-3.9 AHD': 'LAT_map',
             '0 AHD': 'MSL_map'}
-
-damage_maps = {'1.5 AHD': 'HAT_damage',
-               '-1.5 AHD': 'LAT_damage',
-               '0 AHD': 'MSL_damage'}
 
 gauge_map = 'pt_hedland_gauge_map.jpg'
@@ -91,14 +89,15 @@
     swwfiles[swwfile] = label_id
 
-texname = sww2timeseries(swwfiles,
-                         project.gauge_filename,
-                         #label_id,
-                         production_dirs,
-                         report = True,
-                         plot_quantity = ['stage', 'speed'],
-                         time_min = None,
-                         time_max = None,
-                         title_on = False,
-                         verbose = True)
+texname, elev_output = sww2timeseries(swwfiles,
+                                      project.gauge_filename,
+                                      production_dirs,
+                                      report = True,
+                                      reportname = 'latexoutput',
+                                      plot_quantity = ['stage', 'speed'],
+                                      surface = False,
+                                      time_min = None,
+                                      time_max = None,
+                                      title_on = False,
+                                      verbose = True)
 
 latex_output.append(texname)
@@ -117,4 +116,5 @@
 % * an introduction must be written in introduction.tex; a basic outline and
 %   some of the core inputs are already in place
+% * outline of the modelling methodology provided in modelling_methodology.tex
 % * the tsunami-genic event should be discussed in tsunami_scenario.tex
 % * an computational_setup.tex file needs to be written for the particular scenario
@@ -157,19 +157,23 @@
   \maketitle
   
-  \\begin{abstract}
-    \input{abstract}
-  \end{abstract}
-  
   \\tableofcontents
+
+   \section{Executive Summary}
+    \label{sec:execsum}
+  \input{execsum}
   
   \section{Introduction}
     \label{sec:intro}
   \input{introduction}
-   
+
+   \section{Modelling methodology}
+    \label{sec:methodology}
+    \input{modelling_methodology}
+    
   \section{Tsunami scenarios}
     \label{sec:tsunamiscenario}
     \input{tsunami_scenario}
 
-   \section{Inundation Model}
+   \section{Inundation model}
     \label{sec:anuga}
     \input{anuga}
@@ -213,5 +217,5 @@
 
 s  = """
-\caption{Point locations used for Pt Hedland study.}  
+\caption{Point locations used for Onslow study.}  
 \label{fig:points}
 \end{figure}
@@ -235,14 +239,14 @@
 
 s = """
-   \section{Damage modelling}
-    \label{sec:damage}
+   \section{Impact modelling}
+    \label{sec:impact}
      \input{damage}
 """
 fid.write(s)
 
-for i, name in enumerate(production_dirs.keys()):
-
-    s = '\input{%s} \n \clearpage \n \n' %damage_maps[production_dirs[name]]
-    fid.write(s)
+#for i, name in enumerate(production_dirs.keys()):
+
+#    s = '\input{%s} \n \clearpage \n \n' %damage_maps[production_dirs[name]]
+#    fid.write(s)
 
 s = """
Index: /production/pt_hedland_2006/project.py
===================================================================
--- /production/pt_hedland_2006/project.py	(revision 3189)
+++ /production/pt_hedland_2006/project.py	(revision 3190)
@@ -108,8 +108,8 @@
 # region to export (used from export_results.py)
 
-e_min_area = 659000#633000
-e_max_area = 678000#690000
-n_min_area = 7746000#7740000
-n_max_area = 7757000#7761000
+e_min_area = 633000
+e_max_area = 690000
+n_min_area = 7740000
+n_max_area = 7761000
 
 refzone = 50 # confirm with Hamish
Index: /production/sydney_2006/export_results.py
===================================================================
--- /production/sydney_2006/export_results.py	(revision 3189)
+++ /production/sydney_2006/export_results.py	(revision 3190)
@@ -4,5 +4,5 @@
 from pyvolution.ermapper_grids import read_ermapper_grid
 
-name = project.outputname2
+name = project.outputname
 
 #print 'Which variable do you want to export?'
Index: /production/sydney_2006/find_max.py
===================================================================
--- /production/sydney_2006/find_max.py	(revision 3189)
+++ /production/sydney_2006/find_max.py	(revision 3190)
@@ -35,29 +35,34 @@
 c = 0
 direction = 'negative'
+xstar = x0
+print 'start of loop', deriv
 while c < 100000 and deriv > 0:
     deriv = dfuncdx(x,0.83)
-    if deriv < 0: xstar = x
+    if deriv < 0:
+        xstar = x
+        print 'hello', deriv, xstar/lam0, c
     if direction == 'positive': x += step
     if direction == 'negative': x -= step
     c += 1
+    
 
 print 'location of maximum of surface elevation function: xstar = %f' % (xstar/lam0)
 const = 1.0 #a3D = ? #sydney = 86? and kappad=1
 
-x = arange(-20,25,0.001)
-y = arange(-10,10,0.1)
-X,Y = meshgrid(x,y)
+#x = arange(-20,25,0.001)
+#y = arange(-10,10,0.1)
+#X,Y = meshgrid(x,y)
 
-test = 0.0
-for xi in x:
-    testi = func(xi,0.0,0.83)
-    if direction == 'positive':
-        if testi > test:
-            test = testi
-            xstar = xi
-    else:
-        if testi < test:
-            test = testi
-            xstar = xi
+#test = 0.0
+#for xi in x:
+#    testi = func(xi,0.0,0.83)
+#    if direction == 'positive':
+#        if testi > test:
+#            test = testi
+#            xstar = xi
+#    else:
+#        if testi < test:
+#            test = testi
+#            xstar = xi
     
-print 'check: xstar = %f and eta(xstar) = %f' %(xstar/lam0, test)
+#print 'check: xstar = %f and eta(xstar) = %f' %(xstar/lam0, test)
Index: /production/sydney_2006/get_timeseries.py
===================================================================
--- /production/sydney_2006/get_timeseries.py	(revision 3189)
+++ /production/sydney_2006/get_timeseries.py	(revision 3190)
@@ -5,15 +5,17 @@
 import project
 from pyvolution.util import file_function
-from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
+#from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
 from pylab import *
 from matplotlib.ticker import MultipleLocator, FormatStrFormatter
 
-swwfile = project.outputname + '.sww' 
-
-
+swwfile = project.outputname + '.sww'
+#swwfile = project.outputdir + timestampdir + sep + project.outputname + '.sww'
+# place to store figures
+#graphloc = project.outputdir + timestampdir + sep
+graphloc = project.outputdir
 
 #Time interval to plot
-tmin = 13000
-tmax = 21000
+tmin = 2*60
+tmax = 10*60
 
 def get_gauges_from_file(filename):
@@ -34,5 +36,4 @@
         gaugelocation.append(location)
         
-    #Return gauges and raw data for subsequent storage
     return gauges, lines, gaugelocation
 
@@ -40,5 +41,5 @@
 gauges, lines, locations = get_gauges_from_file(project.gauge_filename)
 
-print 'number of gauges for Benfield:   ', len(gauges)
+print 'number of gauges:   ', len(gauges)
 
 #Read model output
@@ -50,6 +51,6 @@
                   use_cache = True)
 
-
-               
+print 'size f', size(f.quantities['stage'],axis=0), size(f.quantities['stage'],axis=1)
+
 from math import sqrt, atan, degrees
 from Numeric import ones
@@ -85,7 +86,7 @@
         uh = f(t, point_id = k)[2]
         vh = f(t, point_id = k)[3]
-        myloc = locations[k]
+        gaugeloc = locations[k]
         depth = w-z
-
+        
         m = sqrt(uh*uh + vh*vh)   #Absolute momentum
         vel = sqrt(uh*uh + vh*vh) / (w-z + 1.e-30) #Absolute velocity
@@ -135,5 +136,5 @@
              model_time, elevations, '-k')
     #name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])
-    name = 'Gauge_%d: (%.1f, %.1f) Location: %s' %(k, g[0], g[1], myloc)
+    name = 'Gauge_%d: (%.1f, %.1f) Location: %s' %(k, g[0], g[1], gaugeloc)
     title(name)
 
@@ -144,8 +145,7 @@
            shadow=True,
            loc='upper right')
-    #savefig('Gauge_%d_stage' %k)
-    savefig('Gauge_%s_stage' %myloc)
-
-    # raw_input('Next')
+    #('Gauge_%d_stage' %k)
+    savefig('%sGauge_%s_stage' %(graphloc, gaugeloc))
+    #savefig('Gauge_%s_stage.eps' %gaugeloc)
     
     #Momentum plot
@@ -159,8 +159,6 @@
     ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')    
     #savefig('Gauge_%d_momentum' %k)
-    savefig('Gauge_%s_momentum' %myloc)
+    savefig('%sGauge_%s_momentum' %(graphloc, gaugeloc))
     
-    # raw_input('Next')
-
     #Bearing plot
     ion()
@@ -181,7 +179,5 @@
     ylabel(' atan(vh/uh) [degrees from North]')    
     #savefig('Gauge_%d_bearing' %k)
-    savefig('Gauge_%s_bearing' %myloc)
-    
-    # raw_input('Next')
+    savefig('%sGauge_%s_bearing' %(graphloc, gaugeloc))
 
     #Speed plot
@@ -195,9 +191,7 @@
     ylabel('sqrt( uh^2 + vh^2 ) / depth [m/s]')    
     #savefig('Gauge_%d_speed' %k)
-    savefig('Gauge_%s_speed' %myloc)
-    
-    # raw_input('Next')
-
-    whichone = '_%s' %myloc
+    savefig('%sGauge_%s_speed' %(graphloc, gaugeloc))
+
+    whichone = '_%s' %gaugeloc
     thisfile = project.gaugetimeseries+whichone+'.csv'
     fid = open(thisfile, 'w')
Index: /production/sydney_2006/plot_integral.py
===================================================================
--- /production/sydney_2006/plot_integral.py	(revision 3189)
+++ /production/sydney_2006/plot_integral.py	(revision 3190)
@@ -30,4 +30,5 @@
 hold(False)
 plot(t, integral, '.-b')
+axis([0, max(t), min(integral)-0.1, min(integral) + 1])
 title('Checking for water conservation')
 xlabel('Time (sec)')
Index: /production/sydney_2006/process_gauges.py
===================================================================
--- /production/sydney_2006/process_gauges.py	(revision 3189)
+++ /production/sydney_2006/process_gauges.py	(revision 3190)
@@ -5,5 +5,5 @@
 import project
 from pyvolution.util import file_function
-from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
+#from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
 from pylab import *
 
@@ -11,5 +11,5 @@
 plot = False
     
-swwfile = project.outputname + '.sww'
+swwfile = project.outputname4 + '.sww'
 
 def get_gauges_from_file(filename):
Index: /production/sydney_2006/project.py
===================================================================
--- /production/sydney_2006/project.py	(revision 3189)
+++ /production/sydney_2006/project.py	(revision 3190)
@@ -30,6 +30,6 @@
 nmaxviz = 6283000
 
-basename = 'slump_315res'
-basename2 = 'slump_1000res'
+basename = 'slump_duncan'
+basename2 = 'slump_1000res_1min'
 basename4 = 'slump_poly_ingo_test'
 
@@ -62,4 +62,5 @@
 
 gauge_filename = gaugedir + 'sydney_gauges.xya'
+buildings_filename = gaugedir + 'all_bld_ind.csv'
 #gauge_filename = gaugedir + 'sydney_gauges_test.xya'
 gauge_outname = gaugedir + 'gauges_max_output.xya'
@@ -137,10 +138,23 @@
 
 # clipping used for demo purposes to fit around poly3 and poly4 (Ingo's files)
-j0 = [318000, 6249000]
-j1 = [387000, 6249000]
-j2 = [387000, 6270000]
-j3 = [318000, 6270000]
-
-demopoly = [j0, j1, j2, j3]
+#j0 = [318000, 6249000]
+#j1 = [387000, 6249000]
+#j2 = [387000, 6270000]
+#j3 = [318000, 6270000]
+
+#demopoly = [j0, j1, j2, j3]
+
+# second demo poly which isn't rectangular
+j0 = [385000, 6280000]
+j1 = [360000, 6272500]
+j2 = [335000, 6272500]
+j3 = [330000, 6265000]
+j31 = [325000, 6260000]
+j4 = [316000, 6260000]
+j5 = [316000, 6247000]
+j6 = [350000, 6247000]
+j7 = [385000, 6238000]
+
+demopoly = [j0, j1, j2, j3, j31, j4, j5, j6, j7]
 
 # to put chunk back in
Index: /production/sydney_2006/run_sydney_smf.py
===================================================================
--- /production/sydney_2006/run_sydney_smf.py	(revision 3189)
+++ /production/sydney_2006/run_sydney_smf.py	(revision 3190)
@@ -27,4 +27,5 @@
 from pyvolution.quantity import Quantity
 from Numeric import allclose
+from utilities.polygon import inside_polygon
 
 # Application specific imports
@@ -70,17 +71,56 @@
 
 # original issue to Benfield
-#interior_res = 5000
-#interior_regions = [[project.harbour_polygon_2, interior_res],
-#                    [project.botanybay_polygon_2, interior_res]]
+interior_res = 5000
+interior_regions = [[project.harbour_polygon_2, interior_res],
+                    [project.botanybay_polygon_2, interior_res]]
 
 # used for finer mesh
-interior_res1 = 5000
-interior_res2 = 315
-interior_regions = [[project.newpoly1, interior_res1],
-                    [project.south1, interior_res1],
-                    [project.finepolymanly, interior_res2],
-                    [project.finepolyquay, interior_res2]]
-
-print 'number of interior regions', len(interior_regions)
+#interior_res1 = 5000
+#interior_res2 = 315
+#interior_regions = [[project.newpoly1, interior_res1],
+#                    [project.south1, interior_res1],
+#                    [project.finepolymanly, interior_res2],
+#                    [project.finepolyquay, interior_res2]]
+
+# used for coastal polygons constructed by GIS guys
+def get_polygon_from_file(filename):
+    """ Function to read in output from GIS determined polygon
+    """
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+    
+    polygon = []
+    for line in lines[1:]:
+       fields = line.split(',')
+       x = float(fields[1])
+       y = float(fields[2])
+       polygon.append([x, y])
+        
+    return polygon
+
+num_polygons = 9
+fileext = '.csv'
+filename = project.polygonptsfile
+
+#interior_res = 1000
+#interior_regions = []
+#bounding_polygon = project.diffpolygonall#project.demopoly
+#count = 0
+#for p in range(1, num_polygons+1):
+#    thefilename = filename + str(p) + fileext
+#    print 'reading in polygon points', thefilename
+#    interior_polygon = get_polygon_from_file(thefilename)
+#    interior_regions.append([interior_polygon, interior_res])
+#    n = len(interior_polygon)
+#    # check interior polygon falls in bounding polygon
+#    if len(inside_polygon(interior_polygon, bounding_polygon,
+#                   closed = True, verbose = False)) <> len(interior_polygon):
+#        print 'WARNING: interior polygon %d is outside bounding polygon' %(p)
+#        count += 1
+#    # check for duplicate points in interior polygon
+    
+print 'number of interior polygons: ', len(interior_regions)
+#if count == 0: print 'interior polygons OK'
 
 #FIXME: Fix caching of this one once the mesh_interface is ready
@@ -89,9 +129,16 @@
 # original + refined region
 _ = cache(create_mesh_from_regions,
-          project.diffpolygonall,
+          #project.demopoly,
+          project.diffpolygonall2,
+          #{'boundary_tags': {'bottom': [0], 
+          #                   'right1': [1], 'right0': [2],
+          #                   'right2': [3], 'top': [4], 'left1': [5],
+          #                   'left2': [6], 'left3': [7]},
           {'boundary_tags': {'bottom': [0], 
-                             'right1': [1], 'right0': [2],
-                             'right2': [3], 'top': [4], 'left1': [5],
+                             'bottom1': [1], 'right': [2],
+                             'top1': [3], 'top': [4], 'left1': [5],
                              'left2': [6], 'left3': [7]},
+          #{'boundary_tags': {'bottom': [0], 'right': [1],
+          #                   'top': [2], 'left': [3]},
            'maximum_triangle_area': 100000,
            'filename': meshname,           
@@ -126,9 +173,9 @@
                                radius=3330,
                                dphi=0.23,
-                               x0=project.slump_origin[0], 
-                               y0=project.slump_origin[1], 
+                               x0=project.slump_origin2[0], 
+                               y0=project.slump_origin2[1], 
                                alpha=0.0, 
-                               domain=domain)
-
+                               domain=domain,
+                               verbose=True)
 
 #-------------------------------------------------------------------------------                                 
@@ -138,5 +185,5 @@
 # apply slump after 30 mins, initialise to water level of tide = 0
 domain.set_quantity('stage', 0.0)
-domain.set_quantity('friction', 0.03) 
+domain.set_quantity('friction', 0.04) 
 domain.set_quantity('elevation',
                     filename = project.combineddemname + '.pts',
@@ -158,12 +205,16 @@
 #                      'right2': Br, 'top': Br, 'left1': Br,
 #                      'left2': Br, 'left3': Br} )
-
+# for new tests 4 April 2006
+#domain.set_boundary( {'bottom': Br, 'bottom1': Br, 'right': Br,
+#                      'top1': Br, 'top': Br, 'left1': Br,
+#                      'left2': Br, 'left3': Br} )
+                             
 # enforce Dirichlet BC - from 30/03/06 Benfield visit
-domain.set_boundary( {'bottom': Bd, 'right1': Bd, 'right0': Bd,
-                      'right2': Bd, 'top': Bd, 'left1': Bd,
+domain.set_boundary( {'bottom': Bd, 'bottom1': Bd, 'right': Bd,
+                      'top1': Bd, 'top': Bd, 'left1': Bd,
                       'left2': Bd, 'left3': Bd} )
 
 # increasingly finer interior regions
-#domain.set_boundary( {'bottom': Br, 'right': Br, 'left': Br, 'top': Br} )
+#domain.set_boundary( {'bottom': Bd, 'right': Bd, 'left': Bd, 'top': Bd} )
 
 
@@ -177,19 +228,67 @@
 fid = open(thisfile, 'w')
 
-for t in domain.evolve(yieldstep = 120, finaltime = 18000):
+# save every 10 secs leading up to slump initiation
+for t in domain.evolve(yieldstep = 10, finaltime = 60): # 6 steps
     domain.write_time()
     domain.write_boundary_statistics(tags = 'bottom')
-
-    # calculate integral
-    thisstagestep = domain.get_quantity('stage') 
-    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
-    fid.write(s)
-
-    # add slump after 30 mins
-    if allclose(t, 30*60):
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+    # add slump after 30 secs
+    if allclose(t, 30):
         slump = Quantity(domain)
         slump.set_values(tsunami_source)
         domain.set_quantity('stage', slump + thisstagestep)
-    
+        #test_stage = domain.get_quantity('stage')
+        #test_elevation = domain.get_quantity('elevation')
+        #test_depth = test_stage - test_elevation
+        #test_max = max(test_depth.get_values())
+        #print 'testing', test_max
+
+import sys; sys.exit()
+
+# save every two minutes leading up to interesting period
+for t in domain.evolve(yieldstep = 120, finaltime = 660, # steps
+                       skip_initial_step = True): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom')
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+
+
+# save every thirty secs during interesting period
+for t in domain.evolve(yieldstep = 60, finaltime = 5000, # steps
+                       skip_initial_step = True):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+
+
+import sys; sys.exit()
+# save every two mins for next 5000 secs
+for t in domain.evolve(yieldstep = 120, finaltime = 10000, # about 42 steps
+                       skip_initial_step = True):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+    
+# save every half hour to end of simulation    
+for t in domain.evolve(yieldstep = 1800, finaltime = 10*60*60, # 14 steps
+                       skip_initial_step = True):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage'
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
     
 print 'That took %.2f seconds' %(time.time()-t0)
Index: /production/sydney_2006/run_sydney_smf_polyingo.py
===================================================================
--- /production/sydney_2006/run_sydney_smf_polyingo.py	(revision 3189)
+++ /production/sydney_2006/run_sydney_smf_polyingo.py	(revision 3190)
@@ -20,5 +20,5 @@
 
 # Related major packages
-from pyvolution.shallow_water import Domain, Reflective_boundary, Time_boundary
+from pyvolution.shallow_water import Domain, Reflective_boundary, Time_boundary, Dirichlet_boundary
 from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
 from pyvolution.combine_pts import combine_rectangular_points_files 
@@ -26,4 +26,6 @@
 from pyvolution.quantity import Quantity
 from geospatial_data import *
+from utilities.polygon import inside_polygon
+from Numeric import allclose
 
 # Application specific imports
@@ -81,7 +83,6 @@
 meshname = project.meshname4+'.msh'
 
-interior_res1 = 25000
+interior_res1 = 1000
 interior_res2 = 1315
-print 'interior resolution', interior_res1, interior_res2
 
 # Read in files containing polygon points (generated by GIS)
@@ -104,10 +105,11 @@
     return polygon
 
-num_polygons = 3
+num_polygons = 4
 fileext = '.csv'
 filename = project.polygonptsfile
 
 interior_regions = []
-
+bounding_polygon = project.diffpolygonall
+count = 0
 for p in range(3, num_polygons+1):
     thefilename = filename + str(p) + fileext
@@ -115,8 +117,14 @@
     interior_polygon = get_polygon_from_file(thefilename)
     interior_regions.append([interior_polygon, interior_res1])
-
+    n = len(interior_polygon)
+    # check interior polygon falls in bounding polygon
+    if len(inside_polygon(interior_polygon, bounding_polygon,
+                   closed = True, verbose = False)) <> len(interior_polygon):
+        print 'WARNING: interior polygon %d is outside bounding polygon' %(p)
+        count += 1
+    # check for duplicate points in interior polygon
+    
 print 'number of interior polygons: ', len(interior_regions)
-
-#print 'Number of interior polygons: ', len(interior_regions)
+if count == 0: print 'interior polygons OK'
 
 # original
@@ -132,6 +140,8 @@
 #                    [project.finepolyquay, interior_res2]]
 
+
 #FIXME: Fix caching of this one once the mesh_interface is ready
 from caching import cache
+
 
 #_ = cache(create_mesh_from_regions,
@@ -141,11 +151,12 @@
 #                             'right2': [3], 'top': [4], 'left1': [5],
 #                             'left2': [6], 'left3': [7]},
-#           'maximum_triangle_area': 100000,
+#           'maximum_triangle_area': 150000,
 #           'filename': meshname,           
 #           'interior_regions': interior_regions},
 #          verbose = True)
 
+# for demo construction
 _ = cache(create_mesh_from_regions,
-          project.diffpolygonall_test,
+          project.demopoly,
           {'boundary_tags': {'bottom': [0], 'right': [1],
                              'top': [2], 'left': [3]},
@@ -171,6 +182,4 @@
 domain.set_datadir(project.outputdir)
 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
-
-
 
 
@@ -190,6 +199,4 @@
                                domain=domain)
 
-#print 'testing', tsunami_source.get_integral()
-
 #-------------------------------------------------------------------------------                                 
 # Setup initial conditions
@@ -199,8 +206,7 @@
 domain.set_quantity('stage', 0) #tsunami_source)
 domain.set_quantity('friction', 0.0) # supplied by Benfield, initial value 0.03
-domain.set_quantity('elevation', alpha = 10.0,#G, alpha = 550.0, 
-                    #filename = project.combineddemname + '.pts',
+domain.set_quantity('elevation', 
                     filename = project.coarsedemname + '.pts',
-                    use_cache = False,
+                    use_cache = True,
                     verbose = True)
 
@@ -231,9 +237,8 @@
 
 Br = Reflective_boundary(domain)
-
-#domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
-#                      'right2': Br, 'top': Br, 'left1': Br,
-#                      'left2': Br, 'left3': Br} )
-domain.set_boundary( {'bottom': Br, 'right': Br, 'left': Br, 'top': Br} )
+Bd = Dirichlet_boundary([0, 0, 0])
+
+# for demo
+domain.set_boundary( {'bottom': Bd, 'right': Bd, 'left': Bd, 'top': Bd} )
 #domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
 #                      'right2': Br, 'top': Br, 'left1': Br,
@@ -250,6 +255,5 @@
 fid = open(thisfile, 'w')
     
-# original
-for t in domain.evolve(yieldstep = 1, finaltime = 10):
+for t in domain.evolve(yieldstep = 10, finaltime = 600):
     domain.write_time()
     domain.write_boundary_statistics(tags = 'bottom')
@@ -257,4 +261,48 @@
     s = '%.2f, %.2f\n' %(t, stagestep.get_integral())
     fid.write(s)
+    if allclose(t, 30):
+        slump = Quantity(domain)
+        slump.set_values(tsunami_source)
+        domain.set_quantity('stage', slump + stagestep)
+
+# save every two minutes leading up to interesting period
+for t in domain.evolve(yieldstep = 120, finaltime = 660, # steps
+                       skip_initial_step = True): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom')
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+    
+# save every thirty secs during interesting period
+for t in domain.evolve(yieldstep = 30, finaltime = 3500, # steps
+                       skip_initial_step = True):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+
+# save every two mins for next 5000 secs
+for t in domain.evolve(yieldstep = 120, finaltime = 10000, # about 42 steps
+                       skip_initial_step = True):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+    
+# save every half hour to end of simulation    
+for t in domain.evolve(yieldstep = 1800, finaltime = 10*60*60, # 14 steps
+                       skip_initial_step = True):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage'
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
 
 print 'That took %.2f seconds' %(time.time()-t0)
Index: /production/sydney_2006/run_timeseries.py
===================================================================
--- /production/sydney_2006/run_timeseries.py	(revision 3189)
+++ /production/sydney_2006/run_timeseries.py	(revision 3190)
@@ -14,5 +14,5 @@
 
 # False to generate latex output, True otherwise
-title_on = False
+title_on = True
 
 # generate figures - see sww2timeseries for documentation
Index: /production/sydney_2006/smf_volume_conservation.py
===================================================================
--- /production/sydney_2006/smf_volume_conservation.py	(revision 3189)
+++ /production/sydney_2006/smf_volume_conservation.py	(revision 3190)
@@ -7,6 +7,6 @@
 from Numeric import ones
 from pylab import *
-import scipy
-from scipy.special import erf
+#import scipy
+#from scipy.special import erf
 
 ion()
@@ -39,4 +39,5 @@
 eta2 = func(X,0,1)
 
+#import sys; sys.exit()
 #eta = func(X,Y,kappad)
 #im1 = imshow(eta, cmap=cm.jet, interpolation='nearest')
