# -*- coding: utf-8 -*-
"""
Created on Thu Jan 25 11:29:06 2021
@author: neil
"""

# lab week 1, analytic vs simple numerical integral, 
# NOTE indefinite integrals have little numerical meaning since there is
# always an unknown constant.  We learn that -cos(x) is the integral of
# sin(x), but actually the indefinite integral is -cos(x) + C
# in general, only definite integrals are numerically evaluated

# this program compares the numerical integral of sin(x) with the analytic
# integral (cos(x)) by plotting the numerical integral from 0 to x of sin(x),
# also plot -cos(x) + cos(0), which should be the same as the integral
# This version uses the numpy (scipy) integration method 'quad'
# scipy.integrate.quad(func, start,stop), returns [integral,error estimate]
# it is almost identical to our previous integration program

from scipy.integrate import quad  #scipy is a large package of numerical routines
import numpy as np                #that are designed to work well with numpy arrays
import matplotlib.pyplot as plt

# ************************* header ***********************************
plen    = 4*np.pi      # span of solution plot in radians, 2 cycles
npoints = 81           # number of points

# ********************** prepare an empty plot figure ****************
fig = plt.figure()
ax  = fig.add_subplot(111)

# ********************** main program ********************************
x = np.linspace(0,plen,npoints)      # make a list of radians for plotting sine and cos
ax.plot(x,np.sin(x),'-b')            # plot sine
ax.set_title('sine and -cos(x) from 0 to x [in radians]')
ax.set_ylabel('trig function value, and integral value')
ax.set_xlabel('radians')
ax.plot(x,np.cos(0)-np.cos(x)-1,'-ro')# plot cosine intgral (note the -1 to shift the plot)

# now to integrate numerically using the scipy.integrate.quad method
q = np.zeros([len(x)])
q[0] = -1                            # constant of integration (required to center around zero)
for n in np.arange(0,len(x)-1):      # integrate by adding the area from our start point (0)
    q[n+1] = q[n] + quad(np.sin,x[n],x[n+1])[0]
    # note the [0] above, quad returns a list, we want only the first, [0], which is the integral
ax.plot(x,q,'-k*');                  # plot the integral 
plt.show()                           # don't forget to show the plot on the screen!