Changeset 2406


Ignore:
Timestamp:
Feb 14, 2006, 4:49:34 PM (18 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:
1 edited

Legend:

Unmodified
Added
Removed
  • production/sydney_2006/get_timeseries.py

    r2326 r2406  
    55import project
    66from pyvolution.util import file_function
    7 #from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
    87from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
    98from pylab import *
    10 
    11 
    12    
    13 #swwfile = project.newoutputname + '.sww'
    14 swwfile = project.outputname + '.sww'
     9from matplotlib.ticker import MultipleLocator, FormatStrFormatter
     10
     11swwfile = project.outputname2 + '.sww'
    1512
    1613#Time interval to plot
     
    3633       
    3734    #Return gauges and raw data for subsequent storage
    38     #return gauges, linesfs
    3935    return gauges, lines, gaugelocation
    4036
     
    5248
    5349               
    54 from math import sqrt
     50from math import sqrt, atan, degrees
     51from Numeric import ones
    5552N = len(gauges)
    5653for k, g in enumerate(gauges):
     
    6360    momenta = []
    6461    velocity = []
     62    xmom = []
     63    ymom = []
     64    bearings = []
    6565
    6666    max_depth = 0
    6767    max_momentum = 0
    6868    max_velocity = 0
     69
     70    due_east = 90.0*ones([len(f.T)])
     71    due_west = 270.0*ones([len(f.T)])
     72    maxT = max(f.T)
     73    tstep = maxT/(len(f.T)-1)
     74
    6975    for i, t in enumerate(f.T): # T is a list of times
    7076        #if tmin < t < tmax:
     
    7783        m = sqrt(uh*uh + vh*vh)   #Absolute momentum
    7884        vel = sqrt(uh*uh + vh*vh) / (w-z + 1.e-30) #Absolute velocity
    79         print vel
    80         #dep = w-z
    81         #vel = sqrt(uh*uh + vh*vh) / dep #Absolute velocity
    82        
     85        angle = degrees(atan(vh/(uh+1.e-15)))
     86        print 'angle', angle
     87        if (0 < angle < 90.0):
     88            if vh > 0:
     89                bearing = 90.0 - abs(angle)
     90            if vh < 0:
     91                bearing = 270.0 - abs(angle)
     92        if (-90 < angle < 0):
     93            if vh < 0:
     94                bearing = 90.0 - abs(angle)
     95            if vh > 0:
     96                bearing = 270.0 - abs(angle)
     97        if angle == 0:
     98            bearing = 0.0
     99               
    83100        model_time.append(t)       
    84101        stages.append(w)
     
    86103        momenta.append(m)
    87104        velocity.append(vel)
     105        xmom.append(uh)
     106        ymom.append(vh)
     107        bearings.append(bearing)
    88108
    89109        if w-z > max_depth:
     
    95115            print 'max speed', max_velocity
    96116
    97    
    98 
    99117    #Plot only those gauges that have been inundated by more than a threshold
    100118    #if max_depth < 0.2:
     
    124142    raw_input('Next')
    125143
    126 
     144    #X Momentum plot
     145    ion()
     146    hold(False)
     147    plot(model_time, xmom, '-r')
     148    title(name)
     149
     150    title('%s (x momentum)' %name)
     151    xlabel('time [s]')
     152    ylabel('uh [m^2/s]')   
     153    savefig('Gauge_%d_xmomentum' %k)
     154    #savefig('Gauge_%s_momentum' %myloc)
     155   
     156    raw_input('Next')
     157
     158    #y Momentum plot
     159    ion()
     160    hold(False)
     161    plot(model_time, ymom, '-r')
     162    title(name)
     163
     164    title('%s (y momentum)' %name)
     165    xlabel('time [s]')
     166    ylabel('vh [m^2/s]')   
     167    savefig('Gauge_%d_ymomentum' %k)
     168    #savefig('Gauge_%s_momentum' %myloc)
     169   
     170    raw_input('Next')
     171   
    127172    #Momentum plot
    128173    ion()
     
    139184    raw_input('Next')
    140185
     186    #Bearing plot
     187    ion()
     188    hold(False)
     189    ax = plot(model_time, bearings, '-b', model_time, due_west, '-.b',
     190         model_time, due_east, '-.b')
     191    title(name)
     192    ax = axis([0, maxT, 0, 360])
     193    text(maxT+tstep, 90, 'East')
     194    text(maxT+tstep, 270, 'West')
     195    majorLocator = MultipleLocator(45)
     196    #print 'major', majorLocator[1]
     197    #ax.yaxis.set_major_locator(majorLocator) #'yticklabels', range(30,390,30))
     198    # set(labels,color='g',rotation=45)
     199
     200    title('%s (bearing)' %name)
     201    xlabel('time [s]')
     202    ylabel(' atan(vh/uh) [degrees from North]')   
     203    savefig('Gauge_%d_bearing' %k)
     204    #savefig('Gauge_%s_bearing' %myloc)
     205   
     206    raw_input('Next')
     207
    141208    #Speed plot
    142209    ion()
Note: See TracChangeset for help on using the changeset viewer.