source: anuga_work/development/Hinwood_2008/prepare_time_boundary.py @ 5350

"""

Script for running a breaking wave simulation of Jon Hinwoods wave tank.
Note: this is based on the frinction_ua_flume_2006 structure.


Duncan Gray, GA - 2007 



"""


#----------------------------------------------------------------------------
# Import necessary modules
#----------------------------------------------------------------------------

from Scientific.IO.NetCDF import NetCDFFile
from Numeric import array, zeros, Float
  
from anuga.utilities.numerical_tools import ensure_numeric

from interp import interp

def prepare_time_boundary(filename):
    """Convert benchmark 2 time series to NetCDF tms file.
    This is a 'throw-away' code taylor made for files like
    'Benchmark_2_input.txt' from the LWRU2 benchmark
    """



    print 'Creating', filename

    # Read the ascii (.csv) version of this file
    fid = open(filename[:-4] + '.csv')

    # Skip first line
    line = fid.readline()

    # Read remaining lines
    lines = fid.readlines()
    fid.close()


    N = len(lines)
    T = zeros(N, Float)  #Time
    Q = zeros(N, Float)  #Stage
    X = zeros(N, Float)  #XMomentum
    Y = zeros(N, Float)  #YMomentum

    for i, line in enumerate(lines):
        fields = line.split(',')

        T[i] = float(fields[0])
        Q[i] = depth = float(fields[1])
        X[i] = float(fields[2]) * depth
        try:
            Y[i] = float(fields[3]) * depth
        except:
            pass


    # Create tms NetCDF file

    fid = NetCDFFile(filename, 'w')
    fid.institution = 'Geoscience Australia'
    fid.description = 'Input wave for Benchmark 2'
    fid.starttime = 0.0
    fid.createDimension('number_of_timesteps', len(T))
    fid.createVariable('time', Float, ('number_of_timesteps',))
    fid.variables['time'][:] = T

    fid.createVariable('stage', Float, ('number_of_timesteps',))
    fid.variables['stage'][:] = Q[:]

    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
    fid.variables['xmomentum'][:] =  X[:]

    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
    fid.variables['ymomentum'][:] =  Y[:]

    fid.close()


def combine_velocity_depth(velocity_file, depth_file, out_file):
    """
    
    Convert the rawish velocity and depth values, which have values at
    different times to a csv file, with values at the same time, with
    SI units.

    The format for the velocity file is;
    [time, sec], [x-velocity +ve is towards the wave generator, m/sec],
    [y-velocity], [z-velocity]

    The format for the pressure file is
    [time, sec], [mm above SWL for sensor A], many other sensors...
    """

    # Read velocity file
    vfid = open(velocity_file)
    lines = vfid.readlines()
    vfid.close()

    
    n_velocity = len(lines)
    vtimes = zeros(n_velocity, Float)  #Time
    velocities = zeros(n_velocity, Float)  #

    for i, line in enumerate(lines):
        fields = line.split(',')

        vtimes[i] = float(fields[0])
        velocities[i] = float(fields[1])

    # Read the depth file
    dfid = open(depth_file)
    lines = dfid.readlines()
    dfid.close()


    n_depth = len(lines)
    dtimes = zeros(n_depth, Float)  #Time
    depths = zeros(n_depth, Float)  #
    
    for i, line in enumerate(lines):
        fields = line.split(',')

        dtimes[i] = float(fields[0])
        depths[i] = float(fields[1])

    depths_at_vtimes = interp(dtimes, depths, vtimes, missing=1e+20)
    depths_at_vtimes = ensure_numeric(depths_at_vtimes)
    depths_at_vtimes = depths_at_vtimes/1000.00 # convert from mm to m
    velocities = ensure_numeric(velocities)
    velocities = velocities * -1.0 # Swap axis around
    
    fid = open(out_file,'w')

    assert len(depths_at_vtimes) == len(vtimes)
    
    #for vtime, depth_at_vtime, velocity in map(vtimes, depths_at_vtimes,
    #                                           velocities):
    for i in xrange(len(vtimes)):
        fid.write(str(vtimes[i]) + ',' + str(depths_at_vtimes[i]) \
                  + ',' + str(velocities[i])+'\n')

    fid.close()
    
    
#-------------------------------------------------------------------
if __name__ == "__main__":
    combine_velocity_depth('T2R7velfilt.csv','T2R7pressfilt.csv', 'cyeah.csv')