"""Calculate integral.
"""
import sys

def main():
  [fun, a, b] = parse_command_line(sys.argv)
  mesh = read_mesh("mesh.dat")
  a = mesh[0][0]
  b = mesh[-1][0] 
  integral = integrate(mesh, fun)

  print("Integral of %s in [%g %g] is %g" %(fun, a, b, integral))

def parse_command_line(argv):
  """Parse command line"""
  argc = len(sys.argv)
  if len(sys.argv) > 1:
    fun = sys.argv[1]
  else:
    print("Using default function f(x) = x^2")
    fun = "x**2"

  if len(sys.argv) > 3:
    a = float(sys.argv[2])
    b = float(sys.argv[3])
  else:
    print("Using default domian [0,1]")
    a = 0.0
    b = 1.0
 
  return [fun, a, b]

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

def integrate(mesh, fun):
  """Integrate function fun using numerical integration on given mesh"""
  mesh_fun = make_mesh_fun(mesh, fun)
  nelem = len(mesh)-1
  integral = 0.0
  for i in range(nelem):
    integral = integrate_element(i, mesh, mesh_fun)
  return integral

def make_mesh_fun(mesh, fun):
  xcoord = [ node[0] for node in mesh ]
  discrete_fun = [ eval(fun) for x in xcoord ] 
  return discrete_fun 

def integrate_element(i, mesh, mesh_fun):
  x1 = mesh[i][0]
  x2 = mesh[i+1][0]
  f1 = mesh_fun[i]
  f2 = mesh_fun[i+1]
  h = x2 - x1
  integral = h * (f1+f2) / 2.0
  return integral

if __name__ == '__main__':
  main()