Changeset 2406

Timestamp:

Feb 14, 2006, 4:49:34 PM (19 years ago)

Author:

sexton

Message:

Incorporating velocity bearing to determine direction of wave at maximum speed (still to fix up plotting, in particular yaxis ticks)

File:

• production/sydney_2006/get_timeseries.py (modified) (9 diffs)

production/sydney_2006/get_timeseries.py

@@ -5 +5 @@
 import project
 from pyvolution.util import file_function
-#from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
 from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
 from pylab import *
-
-
-    
-#swwfile = project.newoutputname + '.sww'
-swwfile = project.outputname + '.sww'
+from matplotlib.ticker import MultipleLocator, FormatStrFormatter
+
+swwfile = project.outputname2 + '.sww'
 
 #Time interval to plot
@@ -36 +33 @@
         
     #Return gauges and raw data for subsequent storage
-    #return gauges, linesfs
     return gauges, lines, gaugelocation
 
@@ -52 +48 @@
 
                
-from math import sqrt
+from math import sqrt, atan, degrees
+from Numeric import ones
 N = len(gauges)
 for k, g in enumerate(gauges):
@@ -63 +60 @@
     momenta = []
     velocity = []
+    xmom = []
+    ymom = []
+    bearings = []
 
     max_depth = 0
     max_momentum = 0
     max_velocity = 0
+
+    due_east = 90.0*ones([len(f.T)])
+    due_west = 270.0*ones([len(f.T)])
+    maxT = max(f.T)
+    tstep = maxT/(len(f.T)-1)
+
     for i, t in enumerate(f.T): # T is a list of times
         #if tmin < t < tmax:
@@ -77 +83 @@
         m = sqrt(uh*uh + vh*vh)   #Absolute momentum
         vel = sqrt(uh*uh + vh*vh) / (w-z + 1.e-30) #Absolute velocity
-        print vel
-        #dep = w-z
-        #vel = sqrt(uh*uh + vh*vh) / dep #Absolute velocity
-        
+        angle = degrees(atan(vh/(uh+1.e-15)))
+        print 'angle', angle
+        if (0 < angle < 90.0):
+            if vh > 0:
+                bearing = 90.0 - abs(angle)
+            if vh < 0:
+                bearing = 270.0 - abs(angle)
+        if (-90 < angle < 0):
+            if vh < 0:
+                bearing = 90.0 - abs(angle)
+            if vh > 0:
+                bearing = 270.0 - abs(angle)
+        if angle == 0:
+            bearing = 0.0
+                
         model_time.append(t)        
         stages.append(w)
@@ -86 +103 @@
         momenta.append(m)
         velocity.append(vel)
+        xmom.append(uh)
+        ymom.append(vh)
+        bearings.append(bearing)
 
         if w-z > max_depth:
@@ -95 +115 @@
             print 'max speed', max_velocity
 
-    
-
     #Plot only those gauges that have been inundated by more than a threshold
     #if max_depth < 0.2:
@@ -124 +142 @@
     raw_input('Next')
 
-
+    #X Momentum plot
+    ion()
+    hold(False)
+    plot(model_time, xmom, '-r')
+    title(name)
+
+    title('%s (x momentum)' %name)
+    xlabel('time [s]')
+    ylabel('uh [m^2/s]')    
+    savefig('Gauge_%d_xmomentum' %k)
+    #savefig('Gauge_%s_momentum' %myloc)
+    
+    raw_input('Next')
+
+    #y Momentum plot
+    ion()
+    hold(False)
+    plot(model_time, ymom, '-r')
+    title(name)
+
+    title('%s (y momentum)' %name)
+    xlabel('time [s]')
+    ylabel('vh [m^2/s]')    
+    savefig('Gauge_%d_ymomentum' %k)
+    #savefig('Gauge_%s_momentum' %myloc)
+    
+    raw_input('Next')
+    
     #Momentum plot
     ion()
@@ -139 +184 @@
     raw_input('Next')
 
+    #Bearing plot
+    ion()
+    hold(False)
+    ax = plot(model_time, bearings, '-b', model_time, due_west, '-.b',
+         model_time, due_east, '-.b')
+    title(name)
+    ax = axis([0, maxT, 0, 360])
+    text(maxT+tstep, 90, 'East')
+    text(maxT+tstep, 270, 'West')
+    majorLocator = MultipleLocator(45)
+    #print 'major', majorLocator[1]
+    #ax.yaxis.set_major_locator(majorLocator) #'yticklabels', range(30,390,30))
+    # set(labels,color='g',rotation=45)
+
+    title('%s (bearing)' %name)
+    xlabel('time [s]')
+    ylabel(' atan(vh/uh) [degrees from North]')    
+    savefig('Gauge_%d_bearing' %k)
+    #savefig('Gauge_%s_bearing' %myloc)
+    
+    raw_input('Next')
+
     #Speed plot
     ion()