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

"""
Calculate integral of scalar function in 1D
"""
import sys
import os
import itertools
import newtoncotes

class Integrand:
  """Represents integrand as 1D scalar function
  """
  def __init__(self, expression):
    self.expression = expression

  def evaluate(self, x):
    return eval(self.expression)

  def __call__(self, x):
    return self.evaluate(x)

class Mesh:
  """Represents on dimensional mesh
  """
  def __init__(self):
    self.nodes = list()

  def load(self, filename):
    """Read from file a list of nodes"""
    with open(filename, 'r') as f:
      for l in f:
        self.nodes.append(tuple(float(x) for x in l.split()))

  def nelem(self):
    return len(self.nodes)-1

class MeshIntegrator:
  """Calculate integral over a domain discretised by a mesh
  """

  def integrate(self, mesh, integrand):
    """Integrate function fun using numerical integration on given mesh
    """
    quadrature = newtoncotes.TrapezoidQuadrature();
    integral = 0.0
    for i in range(mesh.nelem()):
      integral += self.integrate_element(i, mesh, integrand, quadrature)
    return integral
 
  def make_mesh_fun(self, mesh, fun):
    xcoord = [ node[0] for node in mesh.nodes ]
    discrete_fun = [ fun(x) for x in xcoord ] 
    return discrete_fun 
 
  def integrate_element(self, i, mesh, integrand, quadrature):
    x1 = mesh.nodes[i][0]
    x2 = mesh.nodes[i+1][0]
    qr = quadrature
    integral = 0.0
    for (x, w) in itertools.zip_longest(qr.nodes(x1, x2), qr.weights(x1,x2)):
      integral += w * integrand(x);
    return integral

def parse_command_line(argv):
  """Parse command line"""
  argc = len(sys.argv)
  if len(sys.argv) < 2:
    print("Mesh file name must be given")
    sys.exit(22)
  meshfile = sys.argv[1]
  if len(sys.argv) > 2:
    fun = sys.argv[2]
  else:
    print("Using default function f(x) = x^2")
    fun = "x**2"
 
  return [meshfile, fun]
 
def main():
  [meshfile, fun] = parse_command_line(sys.argv)
  mesh = Mesh()
  mesh.load(meshfile)
  integrand = Integrand(fun)
  integrator = MeshIntegrator()
  integral = integrator.integrate(mesh, integrand)
  print("Integral of %s is %g" %(fun, integral))
 
if __name__ == '__main__':
  main()