source: development/momentum_sink/PCS.py @ 2415

"""Read in sww file, interpolate at specified locations and
    Cross section at specified T

"""
from os import sep
import Numeric
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 *
#from compare_sww import gauge_locations

    
#swwfile = project.newoutputname + '.sww'
#swwfile = project.outputname
swwfile_B = project.outputdir + sep  + 'Buildings_3790.sww'
swwfile_F = project.outputdir + sep  + 'Buildings_3662.sww'
#swwfile_F = project.outputdir + sep  + 'friction_3135.sww'

gauge_depth = Numeric.arrayrange(0, 700, 10)
gauge_breadth = 100
gauge_locations = []

for GD in gauge_depth:
    gauge_location = [GD,gauge_breadth]
    gauge_locations.append(gauge_location)


#Read model output
quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']

f_B = file_function(swwfile_B,
                  quantities = quantities,
                  interpolation_points = gauge_locations,
                  verbose = True,
                  use_cache = True)

f_F = file_function(swwfile_F,
                  quantities = quantities,
                  interpolation_points = gauge_locations,
                  verbose = True,
                  use_cache = True)


T = Numeric.arrayrange(0,500,20)
#T = [ 20, 50, 100, 150, 200 ]

for t in T:
    model_time = []
    stages_B = []
    stages_F = []
    elevations = []
    momenta = []
    velocity = []
    for i, GL in enumerate(gauge_depth):
        wB = f_B(t, point_id = i)[0]
        zB = f_B(t, point_id = i)[1]
        uhB = f_B(t, point_id = i)[2]
        vhB = f_B(t, point_id = i)[3]

        wF = f_F(t, point_id = i)[0]        
        uhF = f_B(t, point_id = i)[2]
        vhF = f_B(t, point_id = i)[3]
        #m = sqrt(uh*uh + vh*vh)   #Absolute momentum
        velB = sqrt(uhB*uhB + vhB*vhB) / (wB + 1.e-30) #Absolute velocity
        velF = sqrt(uhF*uhF + vhF*vhF) / (wF + 1.e-30) #Absolute velocity                   
        #print velB, "buiding velocity"
        #print velF, "friction velocity"

        stages_B.append(wB)
        stages_F.append(wF)
        #elevations.append(zB)  #Should be constant over time
        #momenta.append(mB)
        #velocityB.append(velB)
        #velocityF.append(velF)
        
    #print len(stages), "stages"
    diff = abs(Numeric.array(stages_B) - Numeric.array(stages_F))
    #print diff, "difference"


    ion()
    hold(False)
    plot(gauge_depth, stages_B, '-b',gauge_depth,stages_F,'-r')
    title('Time_%d_ vs water depth' %t)
    xlabel('gauge length(m)')
    ylabel('water depths (m)')
    savefig('Comp_Time%d_B_vs_F' %t)
    #plot(gauge_depth, diff, 'g')
    #title('differences')
    # raw_input('Next')
print "finished"