"""
Ex A.14 in A primer on...

a)Formulate the difference eqs:
s[0] = 0
a[0] = x
s[i] = s[i-1] + a[i-1]
a[i] = -(x**2/((2*i+1)*2*i)) * a[i-1]

Update formula without using arrays:
s = s + a
a = -(x**2/((2*i+1)*2*i)) * a
"""
import numpy as np
import matplotlib.pyplot as plt

#b)
"""
Compute sequence s_0,...s_{n+1}
s_i = S(x,i-1)
Return s_{n+1} = S(x,n)
"""

def sin_Taylor(x,n):
    a = x
    s = 0

    for i in range(1,n+2):
        #s[i] = s[i-1] + a[i-1]
        #a[i] = -(x**2/((2*i+1)*2*i)) * a[i-1]
        s = s + a
        a = -(x**2/((2*i+1)*2*i)) * a

    return s, abs(a)

#simple test to see if the function seems to work:
x = np.pi/2
for n in [1,2,10]:
    print(f'n = {n}, approx = {sin_Taylor(x,n)[0]}')


#c) more rigorous test with a test function
"""
n = 2 by hand:
j = 0 -> x
j = 1 -> -x**3/3!
j = 2 -> x**5/5!
"""

def test_sin_Taylor():
    from math import factorial
    x = 0.5
    tol = 1e-10
    expected = x - x**3/factorial(3)+x**5/factorial(5)
    computed = sin_Taylor(x,2)[0]
    success = abs(expected-computed) < tol
    assert success

test_sin_Taylor()

#d)
x_max = 4*np.pi
x = np.linspace(0,x_max,1001)

for n in range(10):
    y = sin_Taylor(x,n)[0]
    plt.plot(x,y,label=f'n={n}')

plt.plot(x,np.sin(x),label='Exact')
plt.axis([0,x_max,-2,2])
plt.legend()
plt.show()

"""
Terminal> python sin_Taylor_series_diffeq.py
n = 1, approx = 0.9248322292886504
n = 2, approx = 1.0045248555348174
n = 10, approx = 1.0000000000000002
"""