#
# Functie voor het berekenen van de imaginaire error funtie erfi(x)
# op basis van een machtsreeks
#
# Bron: Simpson, Clement en Yeomans, 2003.
# "Analytical Model for Computing Residence Times Near a Pumping Well."
# Ground Water Vol 41 No.3 May-June 2003 pages 351 - 354
#
# Thomas de Meij, november 2011
#

from pylab import *


def fac(n):
# fac(n) berekent n! (n faculteit) met een recursiefunctie

    if n==1 or n==0:
        return 1
    else:
        return n*fac(n-1)

def erfi(x):
# erfi(x) berekent de imaginaire errer functie erfi(x) met een machtsreeks

    if x == 0:
        return 0
    else:
        fout = 1.0
        crit = 1/1000.0
        i = 1
        while fout>crit:
            y = x**(2*i-1)/((2*i-1)*fac(i-1))
            if i == 1:
                ySum = y
            else:
                fout = y/ySum
                ySum = ySum+y
            i = i+1
        return 2/math.sqrt(pi)*ySum

def gErfi(xmin,xmax):
# gErfi maakt een grafiek van erfi(x) tussen xmin en xmax

    x = linspace(xmin,xmax,100)
    f = zeros(len(x))
    for i in range(len(x)):
        f[i] = erfi(x[i])
    figure()
    hold(True)
    xlabel('x')
    ylabel('erfi(x)')
    grid()

    plot(x,f,'b',lw=1)
    show()

if __name__ == '__main__':

    # Wolfram's Mathematica berekent de volgende waarden voor x,Erfi(x):
    Wolfram = [[-2.0,-18.5648],[-1.5, -4.58473],[-1.0, -1.65043],
    [-0.5, -0.614952],[0.0, 0.0],[0.5, 0.614952],[1.0, 1.65043],
    [1.5,4.58473],[2.0, 18.5648]]

    # Print tabel met Wolfram's waarden en uitkomst erfi(x)
    W = Wolfram
    print "           x         Erfi         Erfi"
    print "                  Wolfram  Machtsreeks"
    for i in range(len(W)):
        print '%12.1f'%W[i][0],'%12.5f'%W[i][1],'%12.5f'%erfi(W[i][0])

    # Maak een grafiek van Erfi(x) tussen x = -2.5 en x = 2.5
    gErfi(-2.5,2.5)