import os
from math import sqrt, pi
from shallow_water_1d import *
from Numeric import allclose, array, zeros, ones, Float, take, sqrt
from config import g, epsilon

def analytical_sol(C,t):
    
    #t  = 0.0     # time (s)
    g  = 9.81    # gravity (m/s^2)
    h1 = 10.0    # depth upstream (m)
    h0 = 0.0     # depth downstream (m)
    L = 2000.0   # length of stream/domain (m)
    n = len(C)    # number of cells

    u = zeros(n,Float)
    h = zeros(n,Float)
    x = C-L/2.0
    #x = zeros(n,Float)
    #for i in range(n):
    #    x[i] = C[i]-1000.0

    for i in range(n):
        # Calculate Analytical Solution at time t > 0
        u3 = 2.0/3.0*(sqrt(g*h1)+x[i]/t) 
        h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t)) 

        if ( x[i] <= -t*sqrt(g*h1) ):
            u[i] = 0.0 
            h[i] = h1 
        elif ( x[i] <= 2.0*t*sqrt(g*h1) ):
            u[i] = u3 
            h[i] = h3 
        else:
            u[i] = 0.0 
            h[i] = h0 

    return h , u*h

def newLinePlot(title='Simple Plot'):
    import Gnuplot
    gg = Gnuplot.Gnuplot(persist=0)
    gg.terminal(postscript)
    gg.title(title)
    gg('set data style linespoints') 
    gg.xlabel('x')
    gg.ylabel('y')
    return gg

def linePlot(gg,x1,y1,x2,y2):
    import Gnuplot
    plot1 = Gnuplot.PlotItems.Data(x1.flat,y1.flat,with="linespoints")
    plot2 = Gnuplot.PlotItems.Data(x2.flat,y2.flat, with="lines 3")
    g.plot(plot1,plot2)

#debug = False

print "TEST 1D-SOLUTION III -- DRY BED"

L = 2000.0     # Length of channel (m)
N = 100        # Number of compuational cells
cell_len = L/N # Origin = 0.0

points = zeros(N+1,Float)
for i in range(N+1):
    points[i] = i*cell_len

domain = Domain(points)

def stage(x):
    y = zeros(len(x),Float)
    for i in range(len(x)):
        if x[i]<=1000.0:
            y[i] = 10.0
        else:
            y[i] = 0.0
    return y


import time

finaltime = 30.0
yieldstep = 1
L = 2000.0     # Length of channel (m)
number_of_cells = [100]#,200,500,1000,2000,5000,10000,20000]
h_error = zeros(len(number_of_cells),Float)
uh_error = zeros(len(number_of_cells),Float)
k = 0
for i in range(len(number_of_cells)):
    N = int(number_of_cells[i])
    print "Evaluating domain with %d cells" %N
    cell_len = L/N # Origin = 0.0
    points = zeros(N+1,Float)
    for j in range(N+1):
        points[j] = j*cell_len
        
    domain = Domain(points)
    
    domain.set_quantity('stage', stage)
    domain.set_boundary({'exterior': Reflective_boundary(domain)})
    domain.default_order = 2
    domain.default_time_order = 1
    domain.cfl = 1.0
    domain.limiter = "vanleer"

    t0 = time.time()
    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
        pass
    N = float(N)
    StageC = domain.quantities['stage'].centroid_values
    XmomC = domain.quantities['xmomentum'].centroid_values
    C = domain.centroids
    h, uh = analytical_sol(C,domain.time)
    h_error[k] = 1.0/(N)*sum(abs(h-StageC))
    uh_error[k] = 1.0/(N)*sum(abs(uh-XmomC))
    print "h_error %.10f" %(h_error[k])
    print "uh_error %.10f"% (uh_error[k])
    k = k+1
    print 'That took %.2f seconds' %(time.time()-t0)
    X = domain.vertices.flat
    StageQ = domain.quantities['stage'].vertex_values.flat
    XmomQ = domain.quantities['xmomentum'].vertex_values.flat
    h, uh = analytical_sol(X,domain.time)
    #from pylab import * # this command stuffs up setting of new domain
    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
    print 'test2'
    #hold(False)
    #subplot(211)
    #plot1 = plot(X,h,X,StageQ)
    hold(False)
    plot1 = subplot(211)
    plot(X,h,X,StageQ)
    plot1.set_ylim([-1,11])
    #plot1 = plot(X,StageQ,'o',X,h)
    #set(markersize=1)
    #title('Free Surface Elevation of a Dry Dam-Break')
    xlabel('Position')
    ylabel('Stage')
    legend(('Analytical Solution', 'Numerical Solution'),
           'upper right', shadow=True)
    plot2 = subplot(212)
    plot(X,uh,X,XmomQ)
    plot2.set_ylim([-5,35])
    #title('Xmomentum Profile of a Dry Dam-Break')
    xlabel('Position')
    ylabel('Xmomentum')
    #legend(('Analytical Solution', 'Numerical Solution'),
    #       'upper right', shadow=True)
    file = "dry_bed_"
    file += str(number_of_cells[i])
    file += ".eps"
    #savefig(file)
    show()
    
print "Error in height", h_error
print "Error in xmom", uh_error