# -*- coding: utf-8 -*-
"""
Created on Mon Apr 26 11:45:02 2021
This does the integral/linear algebra version of finding the values of the amplitudes
of fhat.  Note it first does the recursive discrete approximation, with the integral 
method at the end.  The integral method is much shorter and simpler
(Uses only 5 intermediate points in a linear 'fhat' approximation to 'f')
@author: neil
"""

import numpy as np
import matplotlib.pyplot as plt

def points(x,x0,x1,a0,a1):     # function to return the line or a point within an element
    # x is a point or list of points, x0,x1, are the edges of the element, a0,a1 are the
    # values of the 'f' function at the element edges
    temp = (x - x0)/(x1 - x0)    # note that x can be a list of numbers
    return a0*(1-temp) + a1*temp # can return a list, if x is a list

xi = [0.5,1.5,4,5.5,7]                   # list of element boundaries
xpoints = np.linspace(xi[0],xi[-1],100)  # for plotting, make a list of points
f = np.sin(xpoints) + 0.1*xpoints + 2    # this is the 'f' function we are approximating

fig = plt.figure()           # plot the 'f' function in red, leave the axis to plot 'fhat'
ax  = fig.add_subplot(111)  
#ax.plot(xpoints,f,"r-*")  #
ax.set_ylabel("f(x)-red, discrete fhat(x)-green, integral fhat-blue")
ax.set_title("element by element fhat by the integral method")

ao = np.sin(xi[0])+0.1*xi[0]+2
n=5
for i in range(len(xi)-1):              # do each element in sequence
    x = np.linspace(xi[i],xi[i+1],n)    # make n points inside the element
    f = np.sin(x) + 0.1*x + 2           # produce the f curve inside the element
    s = 0                               # s will be our error sum over the lement
    t = 0                               # t the denominator sum (used to avoid a divide by 0)
    for j in range(0,n):
        s += (f[j]*(xi[i+1]-xi[i]) - ao*(xi[i+1]-x[j]))
        t += x[j]-xi[i]
    a1 = s/(t)                          # the division of the two sums is the next 'a'
    fhat = points(x,xi[i],xi[i+1],ao,a1)
    ax.plot(x,fhat,"g-+")               # plot fhat that was returned by points
    ax.plot(x,f,"r-*")                  # plot the f(x) values in this element
    ax.plot([xi[i],xi[i]],[1,3.5],"k:") # add the element edges
    ao = a1
        
ax.plot([xi[-1],xi[-1]],[1,3.5],"k:") 

# now do the integral method, using the fact that the integral of the fhat function
# inside an element is just: a*length of element/2, as developed in class
A = np.zeros([len(xi),len(xi)])  # A and C matrix/vectors for linear algebra
C = np.zeros(len(xi))
C[0] = np.sin(xi[0]) + 0.1*xi[0] + 2  # this sets the a0 amplitude to match the f(x)
A[0,0] = 1
b = np.diff(xi)/2                  # 'b' is the integral of each 'phi's in the element

for i in range(1,len(xi)):         # do each element in sequence

    # integral of 'f(x)' from x1 to x2, (integral done in class)
    x1 = xi[i-1]
    x2 = xi[i]
    If = -np.cos(x2)+np.cos(x1)+0.05*(x2*x2-x1*x1)+2*(x2-x1)
    
    A[i,i]   = b[i-1]              # the A matrix is integrals of phi
    A[i,i-1] = b[i-1]
    C[i] = If                      # the C matrix is integrals of f(x)
    
fhati = np.linalg.solve(A,C)    
ax.plot(xi,fhati,"b-d")