In [9]:
%matplotlib inline
import os
import sys
import platform
import matplotlib
import pandas as pd
import numpy as np
import scipy as sp
import statsmodels.api as sm
import matplotlib.pyplot as plt
import matplotlib.dates as dates
import matplotlib.ticker as tick
from matplotlib.backends.backend_pdf import PdfPages
from datetime import datetime,timedelta
#from pylab import rcParams

#rcParams['figure.figsize'] = 10, 10

In [10]:
print("Operating System " + platform.system() + " " + platform.release())
print("Python Version " + str(sys.version))
print("Pandas Version " + str(pd.__version__))
print("Numpy Version " + str(np.__version__))
print("Scipy Version " + str(sp.__version__))
print("Matplotlib Version " + str(matplotlib.__version__))


Operating System Windows 7
Python Version 2.7.8 (default, Jun 30 2014, 16:03:49) [MSC v.1500 32 bit (Intel)]
Pandas Version 0.17.1
Numpy Version 1.10.4
Scipy Version 0.16.1
Matplotlib Version 1.5.1

In [96]:
def Theis(dist, T, S, pump, t):
    if t <= 0:
        return 0
    else:
        u = (np.power(dist,2)*S)/(4*T*t)
        qpart = (pump/(4.0*np.pi*T))
        lnpart = (-0.5772-np.log(u)+u-(np.power(u,2)/(2.0*np.math.factorial(2)))+(np.power(u,3)/(3*np.math.factorial(3)))-(np.power(u,4)/(4*np.math.factorial(4))))
        return qpart*lnpart

def pumprec(dist,T,S,pump,t,t1,t2):
    return Theis(dist,T,S,pump,t-t1) + Theis(dist,T,S,-1*pump,t-t2)

In [100]:
T = 1000
S = 0.001
dist = 100

rates = [10000,14000,12000,2000,500]
starts = [0,170,250,200,400]
ends = [100,300,500,600,450]

timesteps = np.arange(0,ends[-1]*1.5,1)  

h,g,i = [],[],[]


def TheisPlot(T,S,dist,col):
    h = []
    for t in timesteps:
        pint = 0
        for i in range(len(rates)):
            pint1 = pumprec(dist,T,S,rates[i],t,starts[i],ends[i])
            pint = pint + pint1
        h.append(pint)
    plt.plot(timesteps, h, color=col)

TheisPlot(1000,0.001,100, 'red')    
TheisPlot(1000,0.001,500, 'blue')
TheisPlot(1500,0.001,500, 'green')

plt.scatter(t2, Theis(dist,T,S,-1*pump,t2-dur)+Theis(dist,T,S,pump,t2), color='red')
#plt.ylim(0,1)


Out[100]:
<matplotlib.collections.PathCollection at 0x1af58810>

In [8]:
print(len(timesteps))


1000