#! /usr/bin/env python
"""
 Gef.py is a program  which calculate the Gaussian error function E(x)
 by intgrating exp(-t**2) from 0 to x.   The program calculates E(x)
 for  values of x from 0 to 3 in steps of 0.1.  The inegration use the 
 trapezodial rule.  Final results are displayed as a graph of E(x) vs x
 along with a table of values.

 Paul Eugenio
 PHZ4151C
 13 Feb 2019
"""
from __future__ import division, print_function

import numpy as np
import matplotlib.pyplot as plt

# use latex rendering in lables 
plt.rc('text', usetex=True)

def int_by_Trap(fun, lowerLim, upperLim, N):
    """
    Integration by trapezodial rule using N steps
    This function integrates the function fun(x) 
    from lower limit to upper limit.
    """
    w = (upperLim-lowerLim)/N
    E = 0.5 * w * ( fun( lowerLim ) + fun( lowerLim + N*w ) )
    for k in range(N):
        E += w * fun( lowerLim + k*w )
    return E


#
# main

def errorFunction(t):
    """ integrand function exp(-x^2)"""
    return np.exp(-t*t)

xMin = 0.0        # lower plotting limit
xMax = 3.0        # upper plotting limit
nPoints = 30      # number of points to plot

N = 100           # number of integration slices

# plotting points

xPoints = np.linspace(xMin, xMax, nPoints+1)
yPoints = []#np.zeros(nPoints+1, float)

# now calculate ypoints by intergration of f(t) from 0 to x
print("x\tErrorFunction(x)")

for x in xPoints:
    E = int_by_Trap(errorFunction, xMin, x, N)
    yPoints += [E]
    print("{0:5.3f}\t{1:5.3f}".format(x, E))


# make plot and dress up the figure
plt.plot(xPoints, yPoints)
plt.title("Gaussian error function")
plt.xlabel("x")
plt.ylabel(r"$\int_{0}^{x} e^{-t^2} dt$")
plt.show()