Index: /misc/tools/event_selection/wave_energy/wave_energy.py
===================================================================
--- /misc/tools/event_selection/wave_energy/wave_energy.py	(revision 7554)
+++ /misc/tools/event_selection/wave_energy/wave_energy.py	(revision 7555)
@@ -19,4 +19,5 @@
 import glob
 import numpy
+import run_up
 
 # Change path here!!!
@@ -42,4 +43,5 @@
 max_energy_list = []
 max_energy_inst_list = []
+run_up_list = []
 counter = 0
 
@@ -69,4 +71,6 @@
     max_energy = 0
     max_energy_inst = 0
+    max_stage = 0
+    temp_max_stage = 0
     
     ###############################################################
@@ -85,4 +89,5 @@
         x_vel = x_momentum[i]/(100+stage[i])
         y_vel = y_momentum[i]/(100+stage[i])
+        velocity = numpy.sqrt(x_vel*x_vel,y_vel*y_vel)
         # Kinetic energy KE = 1/2 roh h u^2
         KE = (1./2.)*(numpy.power(x_vel,2)+numpy.power(y_vel,2))*(1000)*(100+stage[i])
@@ -105,7 +110,10 @@
             if temp_max_energy > max_energy:
                 max_energy = temp_max_energy
-                
+            if temp_max_stage > max_stage:
+                max_stage = temp_max_stage
+                max_time = time[i] - start_time   
             energy_sum = energy
             start_time = time[i]
+            max_stage_temp = 0
             
         elif stage[i] < 0 and sign != -1:
@@ -114,5 +122,7 @@
             if temp_max_energy > max_energy:
                 max_energy = temp_max_energy
-                
+            if temp_max_stage > max_stage:
+                max_stage = temp_max_stage
+                max_time = time[i] - start_time                 
             energy_sum = energy
             start_time = time[i]
@@ -120,9 +130,13 @@
         elif stage[i] == 0 and sign != 0:
             sign = 0
+            if temp_max_stage > max_stage:
+                max_stage = temp_max_stage
+                max_time = time[i] - start_time 
             energy_sum = energy
             start_time = time[i]
-            
+                      
         else:
             energy_sum += energy
+            temp_max_stage = max(temp_max_stage,stage[i])
             
     # Add maximum values to list        
@@ -130,5 +144,19 @@
     max_energy_list.append( max_energy)
 
+    # Account for both halves of the wave
+    max_time = max_time*2
+    print 'max stage',max_stage
+    print 'max time',max_time
+    #####################################################################
+    # call run-up module
+    #####################################################################
+    Run_up = run_up.Madsen(100,0.01,max_stage,max_time)
+    print 'runup', Run_up
+    run_up_list.append(Run_up)
+
+    
     counter += 1
+    
+
 
     
@@ -137,4 +165,8 @@
 total_max_inst_energy = 0
 
+
+
+
+
 # Print results for each file and calculate maximum value across all events
 for filename in filepathlist:
@@ -142,5 +174,6 @@
     print 'max crest-integrated energy:', max_energy_list[counter],'J.s/m^2'
     print 'max instantatneous energy:', max_energy_inst_list[counter],'J/m^2'
-
+    print 'Run_up:', run_up_list[counter], 'm'
+    
     if max_energy_list[counter] > total_max_energy:
         total_max_energy = max_energy_list[counter]
@@ -156,2 +189,4 @@
 print '\n File with maximum instantaneous energy is:', filepathlist[max_energy_inst_index]
 print 'Energy is', max_energy_inst_list[max_energy_inst_index], 'J/m^2'
+
+
Index: /misc/tools/event_selection/wave_energy/wave_energy_and_momentum.py
===================================================================
--- /misc/tools/event_selection/wave_energy/wave_energy_and_momentum.py	(revision 7555)
+++ /misc/tools/event_selection/wave_energy/wave_energy_and_momentum.py	(revision 7555)
@@ -0,0 +1,177 @@
+"""
+Program to calcualte the energy density in offshore tsunami time series.
+Energy is integrated over the time of the each wave crest. 
+Also calculates time integrated momentum quantities from the boundary time series.
+Used to investiage different methods of estimating offshore tsunami hazard,
+rather than just wave height.
+
+Energy calculations from Johnson 1997. 'A modern introduction ot the mathematical
+theory of water waves' p23.
+
+Creator: Jonathan Griffin, Geoscience Australia
+Created: 12 October 2009
+"""
+
+import os
+import sys
+from os.path import join
+import csv
+import glob
+import numpy
+
+filepath = r'/nas/gemd/georisk_models/inundation/data/new_south_wales/batemans_bay_tsunami_scenario_2009/anuga/boundaries'
+index = 2872
+filebasename = 'sts_gauge_' + str(index) +'.csv'
+filename = join(filepath, filebasename)
+
+filepathlist = []
+filenamelist = []
+
+
+# Find all relevant gauges files and put into list
+for root, dirs, files in os.walk(filepath):
+
+    for filenames in files:
+        if filenames == filebasename:
+            filepathname = join(root, filenames)
+            filepathlist.append(filepathname)
+            filenamelist.append(filenames)
+            
+# Initialise lists for storing maximum values    
+max_energy_list = []
+max_momentum_list = []
+max_energy_inst_list = []
+counter = 0
+
+###########################################################
+# Loop over gauges files
+###########################################################    
+for filename in filepathlist:
+    time = []
+    stage = []
+    x_momentum = []
+    y_momentum = []
+    momentum = []
+    print 'reading file ', filename
+    csvreader = csv.reader(open(filename))
+    header = csvreader.next()
+    for line in csvreader:  
+        time.append(float(line[0]))
+        stage.append(float(line[1]))
+        x_momentum.append(float(line[2]))
+        y_momentum.append(float(line[3]))
+        # Calculate momentum norm
+        momentum.append(numpy.sqrt(numpy.power(float(line[2]),2)+numpy.power(float(line[3]),2)))
+
+    time_diff = time[1]-time[0]
+
+    sign = 0
+    max_momentum = 0
+    momentum_sum = 0
+    energy_sum = 0
+    max_energy = 0
+    max_energy_inst = 0
+    
+    ###############################################################
+    # Loop through time series and add values (integrate over time)
+    # Sum set to zero when stage height changes sign
+    ###############################################################
+    for i in range(len(time)):
+
+        ###############################################################
+        # Energy calculations - note that quantities from gauge file are
+        # already depth integrated.
+        ###############################################################
+        # Calculate velocities from 'momentum' terms (Note not true momentum,
+        # defined as u*h).
+        x_vel = x_momentum[i]/(100+stage[i])
+        y_vel = y_momentum[i]/(100+stage[i])
+        # Kinetic energy
+        KE = (1./2.)*(numpy.power(x_vel,2)+numpy.power(y_vel,2))*(1000)
+        # Potential energy
+        PE = 9.81*(stage[i])*1000
+        energy = KE+PE
+        print 'KE', KE
+        print 'PE', PE
+        # Maximum instantaneous energy density
+        if energy > max_energy_inst:
+            max_energy_inst = energy
+
+        ###############################################################
+        # If sign of wave changes (goes from trough to crest or vice-versa)
+        # multiply energy sum by timestep and reset sum value to zero
+        ###############################################################
+        
+        if stage[i] > 0 and sign != 1:
+            sign = 1
+            temp_max_momentum = momentum_sum*time_diff
+            temp_max_energy = energy_sum*time_diff
+            if temp_max_momentum > max_momentum:
+                max_momentum = temp_max_momentum
+            if temp_max_energy > max_energy:
+                max_energy = temp_max_energy
+                
+            energy_sum = energy
+            momentum_sum = momentum[i]
+            start_time = time[i]
+            
+        elif stage[i] < 0 and sign != -1:
+            sign = -1
+            temp_max_momentum = momentum_sum*time_diff
+            temp_max_energy = energy_sum*time_diff
+            if temp_max_momentum > max_momentum:
+                max_momentum = temp_max_momentum
+            if temp_max_energy > max_energy:
+                max_energy = temp_max_energy
+                
+            energy_sum = energy
+            momentum_sum = momentum[i]
+            start_time = time[i]
+            
+        elif stage[i] == 0 and sign != 0:
+            sign = 0
+            momentum_sum = momentum[i]
+            energy_sum = energy
+            start_time = time[i]
+            
+        else:
+            momentum_sum += momentum[i]
+            energy_sum += energy
+            
+    # Add maximum values to list        
+    max_energy_inst_list.append(max_energy_inst)
+    max_energy_list.append( max_energy)
+    max_momentum_list.append(max_momentum)
+
+    counter += 1
+
+    
+counter = 0
+total_max_energy = 0
+total_max_momentum = 0
+total_max_inst_energy = 0
+
+# Print results for each file and calculate maximum value across all events
+for filename in filepathlist:
+    print filename
+    print 'max energy', max_energy_list[counter]
+    print 'max instantatneous energy', max_energy_inst_list[counter]
+    print 'max momentum', max_momentum_list[counter]
+
+    if max_energy_list[counter] > total_max_energy:
+        total_max_energy = max_energy_list[counter]
+        max_energy_index = counter
+    if max_energy_inst_list[counter] > total_max_inst_energy:
+        total_max_inst_energy = max_energy_inst_list[counter]
+        max_energy_inst_index = counter
+    if max_momentum_list[counter] > total_max_momentum:
+        total_max_momentum = max_momentum_list[counter]
+        max_momentum_index = counter
+    counter+=1
+
+print '\n File with maximum energy integrated over a wave crest is: ', filepathlist[max_energy_index]
+print 'Energy is ', max_energy_list[max_energy_index], ' J.s/m^2'
+print '\n File with maximum instantaneous energy is: ', filepathlist[max_energy_inst_index]
+print 'Energy is ', max_energy_inst_list[max_energy_inst_index], ' J/m^2'
+print '\n File with maximum momentum is: ', filepathlist[max_momentum_index]
+print 'Momentum is ', max_energy_list[max_momentum_index], 'kg.m/s'
