from pylab import *

def TweeKanalen_fr_metN (x, k, N, L, hL, hR):
    '''Stroming tussen twee kanalen in een freatisch pakket met neerslag
       Formules uit Bear (1979) p.180'''
    def fphi(x,k,N,L,hL,hR):
        phisquare = hL**2-x/L*(hL**2-hR**2)+N/k*(L-x)*x
        return sqrt(phisquare)
        
    phi = fphi(x,k,N,L,hL,hR) 
    Q = k/(2*L)*(hL**2-hR**2) - N*(L/2-x)
    Xws = L/2-k/(2*N*L)*(hL**2-hR**2)
    phiMax = fphi(Xws,k,N,L,hL,hR)
    return [phi,Q,Xws,phiMax]

def gPhi (x,phi,Xws,hL,hR,L,phiMax):
    '''Grafiek van freatische grondwaterstand'''
    figure()
    hold(True)
    
    # limieten assen instellen ivm axvline()
    ymin = floor(0.9*min([hL,hR]))
    ymax = ceil(phiMax)
    axis([0,L,ymin,ymax])
    axvline(x=Xws, linewidth=1, ls='--', color='r')
    
    xlabel('afstand tot de linker watergang (in meters)')
    ylabel('freatische grondwaterstand (in meters)')
    grid() 
    plot(x,phi,'b',lw=2)
    show()

def gQ (x,Q,Xws):
    '''Grafiek van horizontaal debiet in watervoerend pakket'''
    figure()
    hold(True)
    xlabel('afstand tot de linker watergang (in meters)')
    ylabel('debiet in de watervoerende laag (in m2/dag)')
    grid()
    axhline(y=0, linewidth=1, color='b')
    axvline(x=Xws, linewidth=1, ls='--', color='r')
    plot(x,Q,'b',lw=2)
    show()

if __name__ == '__main__':

    k = 3.0
    N = 0.0008
    L = 500.0
    hL = 11.00
    hR = 10.00
    x = linspace(0,L,100)
    [phi,Q,Xws,phiMax] = TweeKanalen_fr_metN(x, k, N, L, hL, hR)

    gPhi(x,phi,Xws,hL,hR,L,phiMax)
    gQ(x,Q, Xws)