# simpson.py

import math

def f(x):
  return math.sin(x)


def simpson(f, a, b, n=500):
  """
  Return the approximation of the integral of f
  from a to b using simpson's rule with n intervals.
  """
  if n % 2 != 0:
    print(f"n={n} must be even integer!")
    n = n+1 
  print("Range:", a, "to", b, "with", n, "points")

  h = (b - a)/n
  sum1 = 0
  for i in range(1, n//2 + 1):
    sum1 += f(a + (2*i-1)*h)
  sum2 = 0
  for i in range(1, n//2):
    sum2 += f(a + 2*i*h)

  return (b-a)/(3*n)*(f(a) + f(b) + 4*sum1 + 2*sum2)



# ------------- main -----------------

if __name__ == "__main__":
  print("Integrate f(x) = sin(x)")
  print(simpson(f, 0, math.pi, 20))  # = 2
  print(simpson(f, 0, math.pi, 50))
  print(simpson(f, 0, math.pi, 1000))