#! /usr/bin/env python
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
#
# Copyright © 2017 putanowr <putanowr@foo>
#
# Distributed under terms of the MIT license.

"""
Newton-Cotes numerical integration formulas
"""

class NewtonCotes:
  def __init__(self, degree):
    self.degree = degree
  def nodes(self, a, b):
    n = self.degree+1
    dx = (b-a)/n
    return [a+i*dx for i in range(n)]

  def weights(self, a, b):
    return ((b-a)*x for x in self.ref_weights)

class TrapezoidQuadrature(NewtonCotes):
    ref_weights = (0.5, 0.5)
    def __init__(self):
        super().__init__(1)

class SimpsonQuadrature(NewtonCotes):
    ref_weights = (x/6.0 for x in (1, 4, 1))
    def __init__(self):
        super().__init__(2)

if __name__ == '__main__':
  qr = TrapezoidQuadrature()
  a = -3;
  b = 5;
  if abs(sum(qr.weights(a,b)) - (b-a)) < 1.e-8:
    print("Tesing weights: OK")
  else:
    print("Testing weights: FAILED")
  print("TrapezoidQuadrature degree: ", qr.degree)

  qr = SimpsonQuadrature()
  if abs(sum(qr.weights(a,b)) - (b-a)) < 1.e-8:
    print("Tesing weights: OK")
  else:
    print("Testing weights: FAILED")
  print("SimpsonQuadrature degree: ", qr.degree)