# -*- 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

import numpy as np
import matplotlib.pyplot as plt

# ************************* header ***********************************
plen    = 4*np.pi         # span of solution plot in radians, 2 cycles
npoints = 71              # 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,'-r')      # plot cosine intgral (note the -1 to shift the plot)
#h2 = text(4,.8,'cosine','Color','red');

# now to integrate numerically using a simple 'trapezoid' scheme to estimate the area under the curve
#and then the areas to get the integral up to that point in x
quad  = -np.ones_like(x)            # create a place to put the integrals and (somewhat cute) add the 
                                    # constant of integration, a 1 (by making it the first of the sum)

for n in np.arange(0,len(x)-1):     # integrate by adding the area from our start point (0)
                                    # to the end point (x) of little
                                    # trapezoids of width dx
    dx = x[n+1] - x[n]              # calc dx each time allows for uneven x spacing
    quad[n+1] = quad[n] + (np.sin(x[n]) + np.sin(x[n+1]))*0.5*dx

ax.plot(x,quad,'-k*')
ax.plot(x,-np.cos(x),'r')
plt.show()