#-------------------------------------------------------------------------------
# Name:        waterloop-met-weerstand.py
# Purpose:     Rise of groundwater table after surface water level rise,
#              water course with entry resitance w (d/m)
# Source:      L.F. Ernst, 1962. Grondwaterstromingen in de verzadigde zone
#              en hun berekening bij aanwezigheid van horizontale evenwijdige
#              open leidingen. Proefschrift, ook verschenen als Verslagen van
#              landbouwkundige onderzoekingen nr 67.15
#              Formula 388 at page 153
# Created:     10-03-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 import exp as exp
from scipy.special import erfc as erfc
from scipy.special import erf as erf
from pylab import *

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

def Ernst1962(x,t,kD,S,w):
    ksi = x/(2*kD*w)
    tau = t/(4*S*kD*w**2)
    return erfc(ksi/(2*sqrt(tau)))-exp(tau+ksi)*erfc(ksi/(2*sqrt(tau))+sqrt(tau))

def plotdh(x,t,dh1,dh2):

    # Define colors
    orange = '#FF9900'
    blue   = '#2F64B2'
    purple = '#000066'
    greyblue = '#EAF3FE'
    mycolors = [orange,blue,purple,greyblue]

    black = '#000000'
    darkgrey = '#333333'
    grey = '#58585A'
    lightgrey = '#808080'
    greys = [black,darkgrey,grey,lightgrey]

    # 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)):
        dh1i = dh1[ti][:]*100
        dh2i = dh2[ti][:]*100
        graph.plot(x, dh1i, color=blue, lw=2.0)
        graph.plot(x, dh2i, color=orange, lw=2.0)

    # Filename
    filename = "waterloop-met-weerstand.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__':

    kD = 700.0
    S = 0.25
    w = 0.15

    # Make a plot
    x = linspace(0,700,100)
    t = array([1,4,16,64])
    [X,T]  = meshgrid(x,t)
    hfrac1 = Ernst1958_1(X,T,kD,S)
    hfrac2 = Ernst1962(X,T,kD,S,w)
    plotdh(x,t,hfrac1,hfrac2)