Index: /anuga_work/production/carnarvon/ph2_compare.py
===================================================================
--- /anuga_work/production/carnarvon/ph2_compare.py	(revision 5994)
+++ /anuga_work/production/carnarvon/ph2_compare.py	(revision 5994)
@@ -0,0 +1,139 @@
+"""
+Read in event time series and determine max stage for Ph2 comparison
+Compare with Green's function (no focussing)
+Compare with ANUGA outputs (full model, 250m no polys, 250m all polys)
+Leharne Fountain and Jane Sexton, 2008
+"""
+
+import project
+from pylab import plot, xlabel, ylabel, savefig, ion, close, axis, title, legend, grid, figure
+from os import sep
+
+###################################################
+# Relevant definitions
+###################################################
+
+def get_max_boundary_data(filename):
+    from anuga.utilities.numerical_tools import mean
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+    stage = []
+    for line in lines[1:]:
+        fields = line.split(',')
+        stage.append(float(fields[3]))
+    return mean(stage)
+
+def get_stage_data(filename,depth,no_models):
+    from anuga.utilities.numerical_tools import mean
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+    
+    mean_stages = zeros((len(depth),no_models), Float)
+    max_stages = zeros((len(depth),no_models), Float)
+    for i in range(no_models):
+        stage5 = []
+        stage10 = []
+        stage20 = []
+        stage50 = []
+        for line in lines[1:]:
+            fields = line.split(',')
+            x = float(fields[0])
+            # based on csv file output from ArcGIS: depth, x,y, no poly, all poly, and orig
+            y = float(fields[i+3]) 
+            if x == -depth[0]:
+                stage5.append(y)
+            if x == -depth[1]:
+                stage10.append(y)
+            if x == -depth[2]:
+                stage20.append(y)
+            if x == -depth[3]:
+                stage50.append(y)
+            
+        mean_stages[:,i] = [mean(stage5), mean(stage10), mean(stage20), mean(stage50)]
+        max_stages[:,i] = [max(stage5), max(stage10), max(stage20), max(stage50)]
+    return mean_stages, max_stages
+
+###################################################
+# Determine max stage at boundary points
+###################################################
+
+event_number = 27283
+boundary_mean = get_max_boundary_data(project.boundaries_dir+str(event_number)+sep+'max_sts_stage.csv')
+
+###################################################
+# Read in max data from all models
+###################################################
+
+from Numeric import zeros, Float
+directory = project.home+project.state+sep+project.scenario+sep+'map_work'
+depth = [5.0,10.,20.,50.]
+no_models = 3
+filename = directory + sep + 'Carnarvon_ph2_compare.csv'
+mean_stages, max_stages = get_stage_data(filename, depth, no_models)
+
+###################################################
+# Compare with Green's function and plot
+###################################################
+
+from anuga.abstract_2d_finite_volumes.util import greens_law
+from Numeric import arange
+d1 = 100.
+d2 = arange(d1,1,-0.1)
+h1 = boundary_mean
+green = []
+for d in d2:
+    h2 = greens_law(d1,d,h1)
+    green.append(h2)
+
+ion()
+figure(1)
+plot(depth,mean_stages[:,2],'>g',d2,green,'-g')
+xlabel('depth (m)') 
+ylabel('stage (m)')
+title('ANUGA outputs (average stage) versus Green\'s approximation \n \
+for event 27283 at Carnarvon')
+legend(['original','Green\'s law'])
+#axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
+grid(True)
+figname = 'ph2compare_Carnarvon_mean_ORIG_' + str(event_number) + '_mean'
+savefig(figname)
+
+figure(2)
+plot(depth,mean_stages[:,1],'+r',depth,mean_stages[:,2],'>g',d2,green,'-g')
+xlabel('depth (m)') 
+ylabel('stage (m)')
+title('ANUGA outputs (average stage) versus Green\'s approximation \n \
+for event 27283 at Carnarvon')
+legend(['250m poly','original','Green\'s law'])
+#axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
+grid(True)
+figname = 'ph2compare_Carnarvon_mean_250AP_' + str(event_number) + '_mean'
+savefig(figname)
+
+figure(3)
+plot(depth,mean_stages[:,0],'ob',depth,mean_stages[:,1],'+r',depth,mean_stages[:,2],'>g',d2,green,'-g')
+xlabel('depth (m)') 
+ylabel('stage (m)')
+title('ANUGA outputs (average stage) versus Green\'s approximation \n \
+for event 27283 at Carnarvon')
+legend(['250m no poly','250m poly','original','Green\'s law'])
+#axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
+grid(True)
+figname = 'ph2compare_Carnarvon_mean_ALL_' + str(event_number) + '_mean'
+savefig(figname)
+
+figure(4)
+plot(depth,max_stages[:,0],'ob',depth,max_stages[:,1],'+r',depth,max_stages[:,2],'>g',d2,green,'-g')
+xlabel('depth (m)') 
+ylabel('stage (m)')
+title('ANUGA outputs (max stage) versus Green\'s approximation \n \
+for event 27283 at Carnarvon')
+legend(['250m no poly','250m poly','original','Green\'s law'])
+#axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
+grid(True)
+figname = 'ph2compare_Carnarvon_max_ALL_' + str(event_number) + '_mean'
+savefig(figname)
+close('all')
+
