source: development/momentum_sink/PCS.py @ 2578

"""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 create_buildings import  WidtH
# from run_friction import fric
from pylab import *

fric = 'Dbl'    
#swwfile = project.newoutputname + '.sww'
#swwfile = project.outputname

#swwfile_B = project.outputdir + sep  + 'Buildings_3790.sww'
swwfile_B = project.outputdir + sep  + 'buildings_Str_Off_BR=15_5_2287.sww'
#swwfile_B = project.outputdir + sep  + 'friction_n=10_3135.sww'
swwfile_F = project.outputdir + sep  + 'friction_n=Dbl_33_0.sww'

# gauge_depth = Numeric.arrayrange(51, 1.5*WidtH, 25) # Random special
gauge_depth = Numeric.arrayrange(0, 5*WidtH, 25)
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', '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,1000,10)
#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):

        # mines data from specified time values in sww file
        wB = f_B(t, point_id = i)[0]
        uhB = f_B(t, point_id = i)[1]
        vhB = f_B(t, point_id = i)[2]

        wF = f_F(t, point_id = i)[0]        
        uhF = f_B(t, point_id = i)[1]
        vhF = f_B(t, point_id = i)[2]
        #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(True)
    subplot(211)
    
    plot(gauge_depth, stages_B, '-b',gauge_depth,stages_F,'-r')
    title('Stage@Time_%d_ vs Gauge length (m)[red: fric, blue: buildings] ' %t)
    xlabel('gauge length(m)')
    ylabel('water depths (m)')
    #savefig('Comp_Time%d_B_vs_F' %t)
    hold(False)
#savefig('Final_comp_n=%s' %fric)
subplot(212)
subplots_adjust(hspace = 0.4)
plot(gauge_depth, diff, 'g')
title('Final Differences in water levels')
xlabel('gauge length(m)')
ylabel('water differences (m)')
savefig('Final_Water_diff_n=%s' %fric)
print "finished"