############################################################################
## CALCULUS - BIOINFORMATICS, 2018-19 T3
## Authors: Toni Guillamon, Alex Haro
############################################################################

##################################################################################################################
### Integrating tables of points, not necessary to know the underlying function f(x)                           ###
### https://docs.scipy.org/doc/scipy/reference/integrate.html#module-scipy.integrate                           ###
##################################################################################################################

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import*

### SOME PRELIMINARIES

### We work with double precision

###from decimal import *
###getcontext().prec = 16
###Decimal(1) / Decimal(7)

### Different ways to define functions
a1=4.0014
b1=0.20549
c1=0.018858

def fsin(x):
  y=np.exp(-c1*x)*(a1+b1*x)
  return y

### A) TRAPEZOIDS' AND SIMPSON'S RULES

## A.1) Trapezoids' rule from Scipy integration package
x = np.linspace(0, 120, num=17, endpoint=True)# Uniformly spaced, but including endpoint; it is n=16
y = fsin(x);
sT1=np.trapz (y,x)

## A.2) A self-made implementation of Trapezoids' rule.

#Warning: Remind to check that you have defined the function f3 that you want.
n=16 #number of intervals
a=0 #endpoints
b=120
h=(b-a)/n #integration step
s=fsin(a)+fsin(b) #we evaluate at the endpoints
for j in range(1,n): #loop to compute the inner components of the trapezoids' rule
  xj=a+h*j
  s=s+2*fsin(xj)

sT2=s*h/2

## A.3) Simpson's rule from Scipy integration package

x = np.linspace(0, 120, num=16+1, endpoint=True) # Uniformly spaced, but including endpoint; it is n=16
y = fsin(x)
sS1=simps(y,x) ### compare the agreement with Slide 8 from Lecture 9.

## A.4) A self-made implementation of Simpson's rule.

#Warning: Remind to check that you have defined the function f3 that you want.
n=16 #number of intervals (must be EVEN)
a=0 #endpoints
b=120
h=(b-a)/n #integration step
s=fsin(a) #we evaluate at the endpoints
for j in range(1,int(n/2)+1): #loop to compute the inner components of the trapezoids' rule
  x1=a+h*(2*j-1)
  x2=a+h*2*j
  s=s+4*fsin(x1)+2*fsin(x2)

s=s-fsin(b) #we subtract an extra f(b) computed inside the loop
sS2=s*h/3

### B) ROMBERG'S METHOD (extrapolation)

#Default values: r=romberg(fsin, 0, 2, tol=1.48e-08, rtol=1.48e-08, show=False, divmax=10, vec_func=False)
#r=romberg(fsin, 0, 2)
#r

### For more details, try the function "RombergScheme" from file "Lab080-2019-RombergMethod.py" to retrieve the table "Extrapolation" given in "trapezoids.ods".


# In[24]:

print(sT1,sT2)
print(sS1,sS2)


# In[36]:

for k in range(1,10):
    n=8
    x = np.linspace(0, 120, num=n*2**k+1, endpoint=True)# Uniformly spaced, but including endpoint; it is n=16
    y = fsin(x);
    sT1=np.trapz (y,x)
    sS1=simps (y,x)
    print(n*2**k,sT1,sS1)


# In[ ]: