# -*- coding: utf-8 -*-
"""
@author: neil
"""
# program to plot a sine curve and its derivative (cosine), and integral (-cosine)
# the cosine is calculated simply as the slope of the sine curve by
# taking (sine(x+dx)-sine(x))/((x+dx)-x)  as dy/dx
# integration is done numerically by the trapezoidal sum rule

# the main reason for this program is to illustrate simple python coding, introduce 'numpy'
# and show a standard way of plotting on the screen using 'matplotlib'

import numpy as np
import matplotlib.pyplot as plt

# start of program, note that comments are GOOD!, also note that numbers are hard coded in this program which
# is poor programing, all numbers should be defined at the beginning of a program in a header section

# start with a simple sine wave plot 
x = np.linspace(0,7,20)               #linspace is a very useful generator, look it up! (google 'numpy linspace')
w = np.sin(x)                         #create a list of x points and the matching sin(x) points, in w

#use the plotting routine we have already developed. After it is created, you can plot onto the axes (ax) at any point in the program
fig = plt.figure()                    #create a place to plot
ax  = fig.add_subplot(111)            #create axes.  you can create multiple plots on the same figure with subplot(2,3,1) example
ax.plot(x,w,'r-o')                    #simple plot, you can add lots of options, common is a string like 'r-o' for a red line with markers
ax.set_title("Sine wave, Derivative and Integral, Finite Difference and 'canned'")
ax.set_xlabel("This is an x label, Radians")
ax.set_ylabel("red-sin, blue-cos, black-FD sin, green-FD integral, magenta-SPIPY")

#Now calculate the derivative using a first order FINITE DIFFERENCE approximation
diff = np.zeros(len(w)-1)             #many ways to make an array, using np.zeros(size) makes a zero filled 1D array
for i in np.arange(0,len(w)-1):       #'for' loop, remember it always has a ':' to denote the following is part of the loop
    diff[i] = w[i+1] - w[i]           #calc the difference between all the adjacent values of the sine curve

delx = x[1]-x[0]	                  #since x is evenly spaced we can define a del x
diff = diff/delx                      #divide the difference by 'dx' and we have an estimate of 'dw/dx'
                                      #it is legal to assign 'diff' after modifying right back into variable 'diff'

# plot our estimate of the derivative (or the cosine) curve
ax.plot(x[:-1]+delx/2,diff,'k-+')     # note the x positions should be in the middle of the steps, so we add delx/2
                                      # also the x[:-1]; is pythonic 'slicing' for 'part of array x, starting at 0 excluding the last
# and for comparison plot the canned numpy version of cosine
ax.plot(x,np.cos(x),'b-d')            # plot the cosine, or the analytic derivative									  
plt.show()                            # Spyder will often let you skip this line, but it is required for a standalone program

# use trapezoid rule to calculate the incremental areas. This is just adding the area under the curve.. in vertical slices
# accumulate the sum of the areas of delx*(w[i+1]+w[i])*.5
iareas = np.zeros(len(w)-1)
suma   = np.zeros(len(w))             # note the areas of all the slices (iareas) is 1 less than length of x
for i in np.arange(0,len(w)-1):   
    iareas[i] = w[i+1]+w[i]           # these are the areas of each little trapezoid under the curve
    suma[i+1] = suma[i] + iareas[i]*0.5*delx  # sum all the areas from 0 to x[i]
                                      # although the number of areas are 1 less than length x, the end points of the slices goes from 0 to len(w)   
ax.plot(x,suma-1,'g-*')               # The definite integral of sine does not include the undetermined constant, so we subtract 1 to match
                                      # the analitic solution [that d(sin)/dx = -cosine]
                                      
#finally we plot the numpy (actualy scipy) integration, using the quad method
from scipy.integrate import quad      #you can import modules anywhere, but it has to be before you use them
                                      #google 'scipy integrate quad' to get documentation
									  #quad() requires a function name as an 'arguement', remember np.sin()is a function but np.sin is the name

integ   = np.zeros(len(w))            # note the areas of all the slices (iareas) is 1 less than length of x
for i in np.arange(1,len(w)):
    integ[i] = quad(np.sin,0,x[i])[0] # quad needs the function to integrate and the limits, it returns the integral plus the error term
                                      # note, passing functions is more advanced than our present state, we look at it next week
ax.plot(x,integ-1,'m-*')              # The definite integral of sine does not include the undetermined constant, so we subtract 1 to match
                                      # the analitic solution [that is d(sin)/dx = -cosine]