source: anuga_work/development/momentum_sink/scripts/PCS.py @ 6985

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

"""
from os import sep
import Numeric
import project
from anuga.pyvolution.util import file_function
from create_buildings import  WidtH
# from run_friction import fric
from pylab import *

fric = 'mid_limit_0'    # FILE LABEL 
#swwfile = project.newoutputname + '.sww'
#swwfile = project.outputname

#swwfile_B = project.outputdir + sep  + 'Buildings_3790.sww'
swwfile_B = project.outputdir + sep  + 'buildings_Str_D=5.8_976.sww'
#swwfile_F = project.outputdir + sep  + 'buildings_Str_half_Off_D=6.89_1116.sww'
#swwfile_B = project.outputdir + sep  + 'friction_n=6_3135.sww'
swwfile_F = project.outputdir + sep  + 'friction_n=0.9_3135.sww'

# gauge_depth = Numeric.arrayrange(51, 1.5*WidtH, 25) # Random special

#gauge_depth = Numeric.arrayrange(0, 5*WidtH, 5)
gauge_depth = [49.600000000000001, 74.6, 99.6, 124.6, 149.6, 174.6, 199.6, 224.6, 240.00999999999999]

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 = False)

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


T = Numeric.arrayrange(700,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)

print stages_F, " stages"
print gauge_depth, " gauge depths"
diff_area_F = trapz(gauge_depth,stages_F)
diff_area_B = trapz(gauge_depth,stages_B)
diff_area = abs(diff_area_F) - abs(diff_area_B)
print diff_area, " <<< diff area*********"
#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"