#-------------------------------------------------------------------------------
# Name:        waterloop-peilverandering.py
# Purpose:     Rise of groundwater table after surface water level rise
#              Ernst(1958) gives three formula's:
#               1. River without entry resitance
#               2. Compensation for evaporation increase (no entry resistance)
#               3. River with entry resitance
#              This script contains formula 1
# Source:      L.F. Ernst, 1959. Verhoging van grondwaterstanden en vermindering
#              van afvoer door opstuwing van beken. Instituut voor Cultuurtechniek
#              en Watrerhuishouding, Wageningen. Mededeling 8.
# Created:     20-02-2012
# Copyright:   (c) Grondwaterformules
# Licence:     Python 2.6 with SciPy
#-------------------------------------------------------------------------------
#!/usr/bin/env python
from scipy import *
from scipy import sqrt as sqrt
from scipy.special import erfc as erfc
from pylab import *

def Ernst1958_1(x,t,kD,S):
    beta = sqrt(kD/S)
    return erfc(x/(2*beta*sqrt(t)))

def plotdh(x,t,dh):

    # Define colors
    orange = '#FF9900'
    blue   = '#2F64B2'

    # Create plots
    graph = subplot(111)

    # Set ticklines
    xticklines = graph.get_xticklines()
    yticklines = graph.get_yticklines()

    for line in xticklines+yticklines:
        line.set_linewidth(3)

    for line in xticklines:
        line.set_visible(False)

    # Set gridlines
    gridlines = graph.get_xgridlines()
    gridlines.extend( graph.get_ygridlines() )
    for line in gridlines:
        line.set_linestyle('-')
    grid(True)

    # Set ticklabels
    xticklabels = graph.get_xticklabels()
    yticklabels = graph.get_yticklabels()

    for label in xticklabels+yticklabels:
        label.set_color('k')
        label.set_fontsize('medium')
        label.set_visible(True)

    # Create annotation
    ylabel('verandering grondwaterstand (%h0)')
    xlabel('afstand tot de watergang (m)')

    # Plot lowring in time
    for ti in range(len(t)):
        dhi = dh[ti,:]*100
        graph.plot(x, dhi, color=blue, lw=1.0)

    # Filename
    filename = "waterloop-peilverandering.png"

    # Save figure to file
    fig = figure(1)
    mydpi = 100.0
    fig.set_dpi(mydpi)
    fig.set_size_inches(6.0,4.0)
    fig.savefig(filename,dpi=mydpi, facecolor='w', edgecolor='w')

    # Show the plot on screen
    show()

if __name__ == '__main__':

    x = linspace(0,700,100)
    t = array([1,4,9,16,25,36,49,64,81,100,121])
    kD = 700.0
    S = 0.25
    [X,T] = meshgrid(x,t)
    hfrac=Ernst1958_1(X,T,kD,S)
    plotdh(x,t,hfrac)