## def animatesww2d_alt(sww_filename, movie_filename, range):
##     """Plot cross section of model output
##     """

##     # Read SWW file	
##     from Scientific.IO.NetCDF import NetCDFFile
##     fid = NetCDFFile(sww_filename, 'r')
    
##     x = fid.variables['x']
##     y = fid.variables['y']
##     volumes = fid.variables['volumes']
    
##     elevation = fid.variables['elevation']
##     time = fid.variables['time']
##     stage = fid.variables['stage']
##     xmomentum = fid.variables['xmomentum']
##     ymomentum = fid.variables['ymomentum']
    
def animatesww2d(swwfile):
    """Read in sww file and plot cross section of model output
    """

    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']

    try:
        fid = open(swwfile)
    except Exception, e:
        msg = 'File "%s" could not be opened: Error="%s"'\
              %(swwfile, e)
        raise msg

    print 'swwfile', swwfile

    # interpolation points are for y = 0
    # number of increments in x
    m = 1000 
    max_x = 100000.
    x = 0
    points = []
    x_points = []
    for i in range(m):
        x += max_x/m
        points.append([x,0.])
        x_points.append(x)
            
    time_thinning = 1
    use_cache = True
    verbose = True

    from anuga.abstract_2d_finite_volumes.util import file_function
    
    f = file_function(swwfile,
                      quantities = sww_quantity,
                      interpolation_points = points,
                      time_thinning = time_thinning,
                      verbose = verbose,
                      use_cache = use_cache)

    n = len(f.get_time()) # number of time steps

    from Numeric import zeros, Float
    from math import sqrt
    stage = zeros((n,m), Float)
    elevation = zeros((n,m), Float)
    xmom = zeros((n,m), Float)
    ymom = zeros((n,m), Float)
    momenta = zeros((n,m), Float)
    speed = zeros((n,m), Float)
    max_stages = zeros((n,m), Float)
    max_stage = None
    min_stages = zeros((n,m), Float)
    min_stage = 100.
    fid = open('stage.csv', 'w')
    for k, location in enumerate(x_points):
        for i, t in enumerate(f.get_time()):
            w = f(t, point_id = k)[0]
            z = f(t, point_id = k)[1]
            uh = f(t, point_id = k)[2]
            vh = f(t, point_id = k)[3]
            depth = w-z      
            m = sqrt(uh*uh + vh*vh)
            if depth < 0.001:
                vel = 0.0
            else:
                vel = m / (depth + 1.e-6/depth)

            stage[i,k] = w
            s = '%.2f,' %(w)
            fid.write(s)
            elevation[i,k] = z
            xmom[i,k] = uh
            ymom[i,k] = vh
            momenta[i,k] = m
            speed[i,k] = vel

            if w > max_stage: max_stage = w
            if w < min_stage: min_stage = w
            max_stages[i,k] = max_stage
            min_stages[i,k] = min_stage
        s = '\n'
        fid.write(s)

    min_stage = min(min(min_stages))
    max_stage = max(max(max_stages))
    stage_axis = [0, max_x, min_stage-1, max_stage+1]
    from pylab import plot, xlabel, ylabel, savefig, close, hold, axis, title,show
    hold(False)
    figs = []
    j = 1
    #max_vec = zeros(100, Float)

    #if sys.platform == 'win32':
    try:
        for i in range(0,n,20):
            #max_vec[i] = max(stage[i,:])
            plot(x_points,stage[i,:]) #x_points,max_vec,'+') 
            xlabel('x (m)')
            ylabel('stage (m)')
            axis(stage_axis)
            name = 'time_%i' %i
            savefig(name)
            title(name)
            figs.append(name)
            show()
    except:
        pass
    

    close('all')

if __name__ == '__main__':
    import sys
    swwfile = sys.argv[1]
    
    animatesww2d(swwfile)
    
#swwfile = join(getenv('INUNDATIONHOME'),'data','validation','convergence_study','5000','myexample2.sww')
#animatesww2d(swwfile)

# once png files have been created, they can then be dragged into
# Windows Movie Maker and exported to a mpeg.